Connection timeout when downloading full-depth LLC4320 regional subset

[Dickerman-coder]

Can anyone help me?

I am not very familiar with Python and the xmitgcm library. While I clearly understand the physical aspects of the 3D data I need for my research, I am struggling with the technical execution of correctly slicing the complex grid and downloading the full-depth (90 layers) LLC4320 data without disconnecting (I keep encountering connection timeouts).

I would greatly appreciate it if you could provide any direct guidance, specific code modifications to safely loop through all depths, or alternative methods to obtain this small 3D regional subset.

[Dickerman-coder]

#!/usr/bin/env python

-- coding: utf-8 --

"""
Download a very small LLC4320 region through direct NASA URLs.

This is a standalone script: it does not import any local project helper
modules. Runtime dependencies still need to be installed in the Python
environment, including xmitgcm, xarray, zarr, fsspec, dask, requests and numpy.

Example:
python 下载LLC4320小区域_直连版.py --load-first --check-values

Default download:
variable : Theta
layer : k=0
time : first output
source face label : 3
xarray face index : 2
region : i=0:50, j=1000:1050
output : LLC4320_small_region_direct.nc
"""

from future import annotations

import argparse
import re
import time
from datetime import datetime
from pathlib import Path
from typing import Iterable

DIRECT_DATA_BASE_URL = "https://data.nas.nasa.gov/ecco/llc_4320/compressed"
DIRECT_GRID_BASE_URL = "https://data.nas.nasa.gov/ecco/llc_4320/grid"
MASK_ZARR_URL = "https://storage.googleapis.com/pangeo-ecco/llc/masks/llc_4320_masks.zarr/"

LLC4320_ITER_START = 10368
LLC4320_ITER_STEP = 144
LLC4320_ITER_STOP = 1495153

def llc4320_data_url(varname: str = "Theta", iternum: int = LLC4320_ITER_START) -> str:
iter_text = f"{iternum:010d}"
return f"{DIRECT_DATA_BASE_URL}/{iter_text}/{varname}.{iter_text}.data.shrunk"

def create_direct_llc4320_model():
"""Create an xmitgcm LLC4320 reader that uses the direct NASA URL layout."""
from fsspec.implementations.http import HTTPFileSystem
from xmitgcm.llcreader import stores
from xmitgcm.llcreader.known_models import LLC4320Model

class DirectNASALLC4320Model(LLC4320Model):
    def __init__(self):
        fs = HTTPFileSystem()
        store = stores.NestedStore(
            fs,
            base_path=DIRECT_DATA_BASE_URL,
            grid_path=DIRECT_GRID_BASE_URL,
            mask_path=MASK_ZARR_URL,
            shrunk=True,
            join_char="/",
        )
        super().__init__(store)

return DirectNASALLC4320Model()

def patch_xarray_update_for_xmitgcm() -> None:
"""
Make xarray.Dataset.update return self for older xmitgcm code paths.

Some newer xarray versions changed Dataset.update into an in-place method
returning None, while the installed xmitgcm still assigns its return value.
"""
import functools

import xarray as xr

original_update = xr.Dataset.update
if getattr(original_update, "_llc4320_returns_self", False):
    return

test_ds = xr.Dataset()
if original_update(test_ds, {}) is not None:
    return

@functools.wraps(original_update)
def update_return_self(self, other=None, **kwargs):
    original_update(self, other, **kwargs)
    return self

update_return_self._llc4320_returns_self = True
xr.Dataset.update = update_return_self

def patch_zarr_array_astype_for_xmitgcm() -> None:
"""
Add a lazy astype wrapper for zarr 3 arrays used by older xmitgcm.

xmitgcm expects zarr.Array.astype to exist. zarr 3 removed it, so this
restores the small lazy behavior xmitgcm needs for mask arrays.
"""
import numpy as np
import zarr

if hasattr(zarr.Array, "astype"):
    return

class LazyAstypeArray:
    def __init__(self, array, dtype):
        self.array = array
        self.dtype = np.dtype(dtype)

    def __getitem__(self, key):
        return self.array[key].astype(self.dtype, copy=False)

    def __bool__(self):
        return True

    @property
    def shape(self):
        return self.array.shape

    @property
    def ndim(self):
        return self.array.ndim

    @property
    def chunks(self):
        return self.array.chunks

def astype(self, dtype, *args, **kwargs):
    return LazyAstypeArray(self, dtype)

zarr.Array.astype = astype

def open_llc4320_dataset(
varnames: Iterable[str] = ("Theta",),
k_levels: Iterable[int] = (0,),
*,
read_grid: bool = False,
grid_vars_to_coords: bool = False,
iter_start: int | None = None,
iter_stop: int | None = None,
iter_step: int | None = None,
iters: Iterable[int] | None = None,
):
patch_xarray_update_for_xmitgcm()
patch_zarr_array_astype_for_xmitgcm()
model = create_direct_llc4320_model()
return model.get_dataset(
varnames=list(varnames),
k_levels=list(k_levels),
read_grid=read_grid,
grid_vars_to_coords=grid_vars_to_coords,
iter_start=iter_start,
iter_stop=iter_stop,
iter_step=iter_step,
iters=None if iters is None else list(iters),
)

def _parse_csv_ints(value: str) -> list[int]:
return [int(item.strip()) for item in value.split(",") if item.strip()]

def apply_xmitgcm_compat_patches() -> None:
"""Apply process-local compatibility patches needed by this environment."""
patch_xarray_update_for_xmitgcm()
patch_zarr_array_astype_for_xmitgcm()

def _slice_existing_dims(
ds, *, face_index, i_start, i_stop, j_start, j_stop, time_index, time_count
):
wanted_slices = {
"face": slice(face_index, face_index + 1),
"i": slice(i_start, i_stop),
"i_g": slice(i_start, i_stop),
"j": slice(j_start, j_stop),
"j_g": slice(j_start, j_stop),
"time": slice(time_index, time_index + time_count),
}
region_slice = {
dim: indexer for dim, indexer in wanted_slices.items() if dim in ds.dims
}
return ds.isel(**region_slice)

def _source_cache_path(cache_dir: Path, varname: str, iternum: int) -> Path:
iter_text = f"{iternum:010d}"
return cache_dir / "compressed" / iter_text / f"{varname}.{iter_text}.data.shrunk"

def _remote_size(session, url: str) -> int | None:
try:
response = session.head(url, allow_redirects=True, timeout=(15, 60))
if response.status_code >= 400:
return None
content_length = response.headers.get("Content-Length")
return int(content_length) if content_length else None
except Exception:
return None

def _content_range_total(value: str | None) -> int | None:
if not value:
return None
match = re.search(r"/(\d+)$", value)
return int(match.group(1)) if match else None

def _download_source_file_resumable(
url: str,
destination: Path,
*,
retries: int,
retry_sleep: int,
) -> None:
import requests

destination.parent.mkdir(parents=True, exist_ok=True)
partial = destination.with_name(destination.name + ".part")

with requests.Session() as session:
    total_size = _remote_size(session, url)
    if destination.exists():
        local_size = destination.stat().st_size
        if total_size is None or local_size == total_size:
            print(f"Source cache exists: {destination} ({local_size / 1024 / 1024:.1f} MB)")
            return
        print(f"Cached source has wrong size; redownloading: {destination}")
        destination.unlink()

    last_error = None
    for attempt in range(1, retries + 1):
        downloaded = partial.stat().st_size if partial.exists() else 0
        headers = {"Range": f"bytes={downloaded}-"} if downloaded else {}
        mode = "ab" if downloaded else "wb"
        print(
            f"Downloading source, attempt {attempt}/{retries}, "
            f"already {downloaded / 1024 / 1024:.1f} MB"
        )
        try:
            with session.get(
                url,
                headers=headers,
                stream=True,
                timeout=(30, 120),
            ) as response:
                if downloaded and response.status_code == 200:
                    print("Server ignored Range header; restarting source download.")
                    downloaded = 0
                    mode = "wb"
                response.raise_for_status()

                range_total = _content_range_total(response.headers.get("Content-Range"))
                if range_total is not None:
                    total_size = range_total
                elif total_size is None and response.headers.get("Content-Length"):
                    total_size = downloaded + int(response.headers["Content-Length"])

                last_report = time.time()
                with open(partial, mode) as handle:
                    for chunk in response.iter_content(chunk_size=1024 * 1024):
                        if not chunk:
                            continue
                        handle.write(chunk)
                        downloaded += len(chunk)
                        now = time.time()
                        if now - last_report >= 5:
                            if total_size:
                                percent = downloaded / total_size * 100
                                print(
                                    f"  {downloaded / 1024 / 1024:.1f}/"
                                    f"{total_size / 1024 / 1024:.1f} MB "
                                    f"({percent:.1f}%)"
                                )
                            else:
                                print(f"  {downloaded / 1024 / 1024:.1f} MB")
                            last_report = now

            final_size = partial.stat().st_size
            if total_size is not None and final_size < total_size:
                raise IOError(
                    f"incomplete source download: {final_size} < {total_size} bytes"
                )
            partial.replace(destination)
            print(f"Source cached: {destination} ({final_size / 1024 / 1024:.1f} MB)")
            return
        except Exception as exc:
            last_error = exc
            print(f"Source download failed on attempt {attempt}: {type(exc).__name__}: {exc}")
            if attempt < retries:
                print(f"Retrying source download after {retry_sleep} seconds...")
                time.sleep(retry_sleep)
    raise last_error

def _cache_source_files(
*,
cache_dir: Path,
varnames: list[str],
iternums: list[int],
retries: int,
retry_sleep: int,
) -> None:
for varname in varnames:
for iternum in iternums:
url = llc4320_data_url(varname, iternum)
destination = _source_cache_path(cache_dir, varname, iternum)
print(f"Source URL: {url}")
_download_source_file_resumable(
url,
destination,
retries=retries,
retry_sleep=retry_sleep,
)

def open_cached_llc4320_dataset(
*,
cache_dir: Path,
varnames: list[str],
k_levels: list[int],
iternums: list[int],
):
from fsspec.implementations.http import HTTPFileSystem
from fsspec.implementations.local import LocalFileSystem
from xmitgcm.llcreader import stores
from xmitgcm.llcreader.known_models import LLC4320Model

class CachedDirectLLC4320Model(LLC4320Model):
    def __init__(self):
        local_fs = LocalFileSystem()
        http_fs = HTTPFileSystem()
        store = stores.NestedStore(
            local_fs,
            base_path=(cache_dir / "compressed").resolve().as_posix(),
            grid_path=DIRECT_GRID_BASE_URL,
            mask_fs=http_fs,
            mask_path=MASK_ZARR_URL,
            shrunk=True,
            join_char="/",
        )
        super().__init__(store)

model = CachedDirectLLC4320Model()
return model.get_dataset(
    varnames=varnames,
    k_levels=k_levels,
    read_grid=False,
    grid_vars_to_coords=False,
    iters=iternums,
)

def _print_variable_stats(ds, label: str) -> None:
import numpy as np

data_vars = list(ds.data_vars)
if not data_vars:
    return
varname = data_vars[0]
values = ds[varname]
print()
print(f"{varname} stats ({label}):")
print(f"  shape = {tuple(values.shape)}")
print(f"  dims  = {values.dims}")
data = values.values
finite = np.isfinite(data)
print(f"  finite/total = {int(finite.sum())}/{data.size}")
print(f"  nan count    = {int(np.isnan(data).sum())}")
if finite.any():
    print(f"  min/max      = {float(np.nanmin(data)):.6g}/{float(np.nanmax(data)):.6g}")
    print(f"  mean         = {float(np.nanmean(data)):.6g}")
else:
    print("  all values are NaN/non-finite")

def _load_with_retries(ds, *, retries: int, retry_sleep: int, progress: bool):
from dask.diagnostics import ProgressBar

last_error = None
for attempt in range(1, retries + 1):
    print(f"Loading remote data, attempt {attempt}/{retries}...")
    try:
        if progress:
            with ProgressBar():
                return ds.load()
        return ds.load()
    except Exception as exc:
        last_error = exc
        print(f"Load failed on attempt {attempt}: {type(exc).__name__}: {exc}")
        if attempt < retries:
            print(f"Retrying after {retry_sleep} seconds...")
            time.sleep(retry_sleep)
raise last_error

def _write_with_retries(
ds,
output: Path,
*,
netcdf_engine: str,
retries: int,
retry_sleep: int,
progress: bool,
) -> None:
from dask.diagnostics import ProgressBar

last_error = None
for attempt in range(1, retries + 1):
    print(f"Writing NetCDF, attempt {attempt}/{retries}...")
    if output.exists():
        output.unlink()
    try:
        if progress:
            with ProgressBar():
                ds.to_netcdf(output, engine=netcdf_engine)
        else:
            ds.to_netcdf(output, engine=netcdf_engine)
        return
    except Exception as exc:
        last_error = exc
        if output.exists():
            output.unlink()
        print(f"Write failed on attempt {attempt}: {type(exc).__name__}: {exc}")
        if attempt < retries:
            print(f"Retrying after {retry_sleep} seconds...")
            time.sleep(retry_sleep)
raise last_error

def backup_existing_file(path: Path) -> None:
if not path.exists():
return
stamp = datetime.now().strftime("%Y%m%d%H%M%S")
backup = path.with_name(f"{path.stem}.bak_{stamp}{path.suffix}")
path.replace(backup)
print(f"Existing output moved to: {backup}")

def download_small_region(
*,
output: Path,
k_levels: list[int],
face_label: int,
i_start: int,
i_stop: int,
j_start: int,
j_stop: int,
time_index: int,
time_count: int,
iternum: int,
iter_step: int,
cache_source: bool,
cache_dir: Path,
check_values: bool,
load_first: bool,
netcdf_engine: str,
progress: bool,
retries: int,
retry_sleep: int,
max_horizontal_points: int,
dry_run: bool,
) -> None:
varnames = ["Theta"]
face_index = face_label - 1
i_count = i_stop - i_start
j_count = j_stop - j_start
if i_count <= 0 or j_count <= 0:
raise ValueError("Region must have positive size: require i_stop > i_start and j_stop > j_start.")
horizontal_points = i_count * j_count
if max_horizontal_points and horizontal_points > max_horizontal_points:
raise ValueError(
f"Requested horizontal region is {horizontal_points} points "
f"({i_count} x {j_count}), larger than --max-horizontal-points "
f"({max_horizontal_points}). Shrink i/j bounds or raise the limit explicitly."
)

print("Example direct data URLs:")
iternums = [iternum + n * iter_step for n in range(time_count)]
for varname in varnames:
    for current_iternum in iternums:
        print(f"  {llc4320_data_url(varname, current_iternum)}")

print()
print("Small-region request:")
print(f"  varnames = {varnames}")
print(f"  k_levels = {k_levels}")
print(f"  face label = {face_label}")
print(f"  face index = {face_index}")
print(f"  i/i_g    = {i_start}:{i_stop}")
print(f"  j/j_g    = {j_start}:{j_stop}")
print(f"  points   = {horizontal_points} ({i_count} x {j_count})")
print(f"  time     = {time_index}:{time_index + time_count}")
print(f"  iternums = {iternums}")
print(f"  output   = {output}")
print(f"  variable = {', '.join(varnames)} only")
print("  grid     = disabled (read_grid=False)")
print(f"  cache    = {cache_source}")
if cache_source:
    print(f"  cache dir = {cache_dir}")
    print("  note      = cache-source downloads the full Theta shrunk file before slicing")
print(f"  engine   = {netcdf_engine}")
print(f"  check    = {check_values}")
print(f"  load     = {load_first}")
print(f"  progress = {progress}")
print(f"  retries  = {retries}")

if dry_run:
    print("  download = no")
    return

apply_xmitgcm_compat_patches()
print("Applied compatibility patches for xarray/zarr + xmitgcm.")

if cache_source:
    _cache_source_files(
        cache_dir=cache_dir,
        varnames=varnames,
        iternums=iternums,
        retries=retries,
        retry_sleep=retry_sleep,
    )
    ds = open_cached_llc4320_dataset(
        cache_dir=cache_dir,
        varnames=varnames,
        k_levels=k_levels,
        iternums=iternums,
    )
else:
    ds = open_llc4320_dataset(
        varnames=varnames,
        k_levels=k_levels,
        read_grid=False,
        grid_vars_to_coords=False,
        iters=iternums,
    )
if ds is None:
    raise RuntimeError(
        "model.get_dataset(...) returned None. This should not happen with "
        "xmitgcm's LLC4320Model. Please check the imported xmitgcm version "
        "and whether local files named xmitgcm.py are shadowing the package."
    )
region = _slice_existing_dims(
    ds,
    face_index=face_index,
    i_start=i_start,
    i_stop=i_stop,
    j_start=j_start,
    j_stop=j_stop,
    time_index=time_index,
    time_count=time_count,
)
print()
print("Selected dataset:")
print(region)
region_to_write = region
if load_first or check_values:
    print()
    print("Loading remote data into memory first...")
    region_to_write = _load_with_retries(
        region, retries=retries, retry_sleep=retry_sleep, progress=progress
    )
    if check_values:
        _print_variable_stats(region_to_write, "loaded data")
print()
print("Writing NetCDF. This is the step that reads remote NASA files unless load=True...")
_backup_existing_file(output)
_write_with_retries(
    region_to_write,
    output,
    netcdf_engine=netcdf_engine,
    retries=retries,
    retry_sleep=retry_sleep,
    progress=progress,
)
print(f"saved: {output}")

def build_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(description="Download a small LLC4320 region.")
parser.add_argument("--output", type=Path, default=Path("LLC4320_small_region_direct.nc"))
parser.add_argument("--k-levels", default="0")
parser.add_argument(
"--face",
type=int,
default=3,
help="Face label used in the original script/figure. Xarray uses face-1.",
)
parser.add_argument("--i-start", type=int, default=0)
parser.add_argument("--i-stop", type=int, default=50)
parser.add_argument("--j-start", type=int, default=1000)
parser.add_argument("--j-stop", type=int, default=1050)
parser.add_argument("--time-index", type=int, default=0)
parser.add_argument("--time-count", type=int, default=1)
parser.add_argument("--iternum", type=int, default=LLC4320_ITER_START)
parser.add_argument("--iter-step", type=int, default=144)
parser.add_argument("--max-horizontal-points", type=int, default=2500)
parser.add_argument("--cache-source", action="store_true")
parser.add_argument("--cache-dir", type=Path, default=Path("LLC4320_source_cache"))
parser.add_argument("--check-values", action="store_true")
parser.add_argument("--load-first", action="store_true")
parser.add_argument("--netcdf-engine", default="scipy", choices=["scipy", "h5netcdf", "netcdf4"])
parser.add_argument("--no-progress", action="store_true")
parser.add_argument("--retries", type=int, default=5)
parser.add_argument("--retry-sleep", type=int, default=10)
parser.add_argument("--dry-run", action="store_true")
return parser

def main() -> int:
args = build_parser().parse_args()
download_small_region(
output=args.output,
k_levels=_parse_csv_ints(args.k_levels),
face_label=args.face,
i_start=args.i_start,
i_stop=args.i_stop,
j_start=args.j_start,
j_stop=args.j_stop,
time_index=args.time_index,
time_count=args.time_count,
iternum=args.iternum,
iter_step=args.iter_step,
cache_source=args.cache_source,
cache_dir=args.cache_dir,
check_values=args.check_values,
load_first=args.load_first,
netcdf_engine=args.netcdf_engine,
progress=not args.no_progress,
retries=args.retries,
retry_sleep=args.retry_sleep,
max_horizontal_points=args.max_horizontal_points,
dry_run=args.dry_run,
)
return 0

if name == "main":
raise SystemExit(main())

[cspencerjones]

I am struggling with the technical execution of correctly slicing the complex grid and downloading the full-depth (90 layers) LLC4320 data without disconnecting

As far as I know, this is not possible. The data is compressed in such a way that you cannot download a small spatial area in the horizontal. I believe that you can download one vertical layer at a time, but you must download the whole globe for each variable, and then throw away the part that you don't want.

The method ECCOPortalLLC4320Model may be helpful to you, but in the most recent xmitgcm versions the path to the data is commented out, because this data was unavailable for a long time (and much of it was corrupted). If you download the package and uncomment the paths:

xmitgcm/xmitgcm/llcreader/known_models.py

Lines 210 to 211 in 115d6b4

	# base_path = 'https://data.nas.nasa.gov/ecco/download_data.php?file=/eccodata/llc_4320/compressed'
	# grid_path = 'https://data.nas.nasa.gov/ecco/download_data.php?file=/eccodata/llc_4320/grid'

you may be able to use the method to download the data that is on the NAS portal (I have not tested this), but you still must download the full globe.