Add fetch_remote_topo; deprecate read_netcdf; refactor etopotools

[mandli]

This PR adds a common function built on the newer topotools that handles remote fetching of files. This functionality replaces the internals of read_netcdf and etopotools, which remain as compatibility shims for now. Deprecation warnings have been added, but can be customized/changed depending on what we are willing to mark as deprecated.

Assisted-by: claude claude-opus-4-8

[mandli]

@rjleveque this was an attempt at removing some of the now duplicated functionality and unifying it. Another PR will also work towards unifying some of the divergent language.

[rjleveque]

This isn't working for me for some reason.

The url I specify gets a local file directory path pre-pended to it somehow, so it is looking for

/Users/rjl/git/CHTuser/sites/LagoonCreek/topo/https:/www.ngdc.noaa.gov/thredds/...

Here's the error:

In [9]: url = 'https://www.ngdc.noaa.gov/thredds/dodsC/regional/crescent_city_1
      ⋮ 3_mhw_2010.nc'
   ...: extent = [-124.12, -124.087, 41.578, 41.608]

In [10]: pdb
Automatic pdb calling has been turned ON

In [11]: topo13s = topotools.fetch_remote_topo(url, crop_extent=extent)
    ...: 
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/venv/geo4/lib/python3.13/site-packages/xarray/backends/file_manager.py:219, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
    218 try:
--> 219     file = self._cache[self._key]
    220 except KeyError:

File ~/venv/geo4/lib/python3.13/site-packages/xarray/backends/lru_cache.py:56, in LRUCache.__getitem__(self, key)
     55 with self._lock:
---> 56     value = self._cache[key]
     57     self._cache.move_to_end(key)

KeyError: [<class 'netCDF4._netCDF4.Dataset'>, ('/Users/rjl/git/CHTuser/sites/LagoonCreek/topo/https:/www.ngdc.noaa.gov/thredds/dodsC/regional/crescent_city_13_mhw_2010.nc',), 'r', (('clobber', True), ('diskless', False), ('format', 'NETCDF4'), ('persist', False)), '1174e964-0c49-4095-976a-192ffc5a183d']

Going up many times in pdb shows that the right thing was passed to xr.open_dataset, so I don't know what's wrong...

> /Users/rjl/git/clawpack/geoclaw/src/python/geoclaw/netcdf_utils.py(441)__init__()
    439         # and naively casting that to float (as a "seconds" value) reads back
    440         # the raw nanosecond tick count -- a billion-fold error.
--> 441         self.ds: xr.Dataset = xr.open_dataset(
    442             self.path, engine=engine, chunks=chunks, mask_and_scale=True,
    443             decode_timedelta=False,

ipdb> self.path
PosixPath('https:/www.ngdc.noaa.gov/thredds/dodsC/regional/crescent_city_13_mhw_2010.nc')
ipdb>

[mandli]

It looks like xarray was mangling the URL before it got to the netCDF stuff, I think this fixes that.

[mandli]

I should have added that this includes some fixes for naming conventions that were confusing.

The attribute crop_extent and crop_bounds already existed on main — what this PR changes is the argument spellings and the buffer type. Here are the tables with an accurate before/after column.

Region / crop vocabulary — before vs. after

Name	Role / frame	Before (main)	After (this PR)
crop_extent	requested crop rectangle · domain	attribute exists; not accepted as an argument	canonical attribute + argument on read, crop, interp_unstructured, fetch_remote_topo
extent	actual data bounds · domain	property with a setter; name also overloaded as a crop input on interp_unstructured	strictly the derived read-only result property
crop_bounds	requested crop · file coords	exists (NetCDF layer) → Fortran nc_crop_bounds	unchanged
crop()	crop operation	takes filter_region=	takes crop_extent= (see below)

Argument renames / deprecations

Call site	Before (main)	After (this PR)	Status
read(...)	filter_region=None	crop_extent=None	filter_region= kept as deprecated alias (warns; TypeError if both)
crop(...)	crop(filter_region=None, ...)	crop(crop_extent=None, ...)	filter_region= kept as deprecated alias
interp_unstructured(...)	extent=None	crop_extent=None	extent= kept as deprecated alias
fetch_remote_topo(...)	did not exist (new in this PR)	crop_extent=None	clean — never spelled filter_region
read_netcdf(...)	extent='all' (active function)	extent='all'	function itself now a deprecated shim over fetch_remote_topo
buffer (on read/crop/Topography)	annotated float, doc already said "int"	int + int()-coerced in crop()	type made honest; integer grid-point count, not meters (cf. interp_unstructured(buffer_length=))

Signatures confirming before → after:

# main
def crop(self, filter_region=None, coarsen=1, buffer=0, align=None):
def read(self, ..., filter_region=None, ...):
def interp_unstructured(self, fill_topo, extent=None, ...):

# this PR
def crop(self, crop_extent=None, coarsen=1, buffer=0, align=None, filter_region=_CROP_EXTENT_UNSET):
def read(self, ..., crop_extent=None, ..., filter_region=_CROP_EXTENT_UNSET):
def interp_unstructured(self, fill_topo, crop_extent=None, ..., extent=_CROP_EXTENT_UNSET):
def fetch_remote_topo(name_or_url, crop_extent=None, coarsen=1, buffer=0, align=None, ...):  # new

[rjleveque]

I'm looking at the latest version and it seems like the topotools.Topography.read() function is missing the argument align and instead of coarsen it has the argument stride.

When reading a netCDF file with topo_type==4 this seems to preserve alignment due to the comment I found in _netcdf_window_indices:

 - *stride* (int) - read stride; the low index is snapped **down** to a
   multiple of *stride* so the strided sub-window shares the same phase as
   striding the full array, keeping the sampled grid identical.

So e.g. if stride = [3,3] then the resulting array is aligned with the full array at the coarser resolution of 3*dx_original.

But when reading an ASCII file with topo_type==3 containing the same data and with essentially the same read command, with the same stride this is not true, and the grids may be misaligned by 1/3 or 2/3 of a coarse grid cell.

The crop function seems to work fine when coarsen and align are used, can we introduce these in the read function so netCDF and ASCII file reads behave the same?

See the attached notebook for an example:

TestTopotools.html

[mandli]

I think there's an incompatibility with the way stride is implemented in the ASCII and netCDF code which then leads to an issue with alignment. It looks like the netCDF is correct though and has just not been implemented fully in the ASCII. I am going to come up with a fix that ensures that these match and test it against your notebook example. Hopefully it won't be too much of a fix, just a little bit more coordinate arithmetic.

[rjleveque]

The way stride is implemented for netCDF, the resulting file is always aligned with the first x and y values in the original file, I think. But that may not always be what is wanted. E.g. some NCEI CUDEM tiles include a few cells to the left and/or right of the integer or quarter degree edges of the quarter degree square patches they cover. In this case when coarsening we might want to specify that the coarsened data should be aligned with the integer lat/lon rather than with a point that's an unknown distance outside.

So my preference would be to implement the align feature in netCDF rather than the stride in ASCII. The old topotools.read_netcdf implemented this, and perhaps could be carried over to read when topotype==4?

But happy to hear other views...

[mandli]

Agreed — that's exactly the direction I took. The align feature is now in read() for topo_type==4 (and for all topo types), carried over from the old read_netcdf. In fact the alignment code is now a single shared helper _crop_indices used by both crop() and the type-4 read, and it's the same argmin-fractional-offset approach the old read_netcdf used so netCDF and ASCII give identical grids by construction.

To your CUDEM point: the old stride behavior (snap to the file's first cell) is exactly what has been deprecated. Now:

• align=None (the default) starts the coarsened grid at the requested crop extent — matching ASCII/crop() and the old read_netcdf.
• align=[-124.25, 47.0] (or integer degrees) locks the coarsened grid to that lattice regardless of where the tile's first cell falls, so a CUDEM tile with a few cells outside the quarter-degree edge coarsens onto the quarter-degree lattice, not onto an unknown offset.

topo.read('tile.nc', topo_type=4, crop_extent=ce, coarsen=3,
          align=[-124.25, 47.0])   # integer/quarter-degree aligned

stride still works but warns and maps onto the scalar coarsen.

One thing to keep in mind with crop(), if the native grid isn't commensurate with the target lattice at the coarse spacing, alignment is the nearest index rather than exact.

I did use your notebook to verify against. Type 3 and type 4 now return identical x, y, Z for the same coarsen+align, across both N→S and S→N latitude storage orders, and the three shifted crops all give the same lattice.