Document new fetch_remote_topo and deprecation of read_netcdf and etopo1_download

doc/changes_to_master.rst

@@ -88,9 +88,9 @@ Changes to geoclaw
   See :ref:`topochanges`.
 
 - **Topography preprocessing attributes.**
-  :class:`~clawpack.geoclaw.topotools.Topography` now supports seven
+  :class:`~clawpack.geoclaw.topotools.Topography` now supports eight
   preprocessing attributes (``crop_extent``, ``coarsen``, ``buffer``,
-  ``align``, ``x_shift``, ``z_shift``, ``negate_z``) that are applied
+  ``align``, ``x_shift``, ``y_shift``, ``z_shift``, ``negate_z``) that are applied
   automatically when :meth:`~clawpack.geoclaw.topotools.Topography.read`
   loads a file.  See :ref:`setrun_topo_preprocessing` for the full table
   and :ref:`topotools` for usage examples and operation order.
@@ -106,6 +106,71 @@ Changes to geoclaw
   are available without loading the elevation array.
   See :ref:`topotools` for an example.
 
+- **Remote-DEM fetch helper;** ``read_netcdf`` **deprecated.**
+  :func:`~clawpack.geoclaw.topotools.fetch_remote_topo` resolves a nickname
+  (a key of ``topotools.remote_topo_urls``), URL, or local path and reads it
+  through the ``topo_type=4`` path, returning a
+  :class:`~clawpack.geoclaw.topotools.Topography` with the requested
+  ``crop_extent``/``coarsen``/``buffer``/``align`` applied.  The older
+  :func:`~clawpack.geoclaw.topotools.read_netcdf` is now a thin shim over it
+  that emits a ``DeprecationWarning``.  In ``etopotools``,
+  :func:`~clawpack.geoclaw.etopotools.fetch_etopo` is a convenience wrapper for
+  the etopo netCDF datasets, and
+  :func:`~clawpack.geoclaw.etopotools.etopo1_download` now always returns a
+  ``Topography``.  Because ``fetch_remote_topo`` uses the ``topo_type=4``
+  reader, the elevation variable must be in meters (or supply ``assume_units``
+  via ``nc_params``).  See :ref:`topo`.
+
+- **Unified coarsening and alignment;** ``stride`` **deprecated.**
+  :meth:`~clawpack.geoclaw.topotools.Topography.read` now accepts
+  ``crop_extent``, ``coarsen``, ``buffer`` and ``align`` directly as keyword
+  arguments; passing one is equivalent to setting the attribute of the same
+  name before calling ``read()``.  ASCII (``topo_type`` 2 or 3) and NetCDF
+  (``topo_type=4``) reads of the same data now return identical ``x``, ``y``
+  and ``Z`` for the same ``coarsen`` and ``align``: both drive off a single
+  shared index computation, so the type-4 lazy hyperslab (including files that
+  store latitude north-to-south) and
+  :meth:`~clawpack.geoclaw.topotools.Topography.crop` agree by construction.
+  Previously the NetCDF-only ``stride`` argument silently did nothing for ASCII
+  reads and used a different alignment convention, so the two paths could
+  disagree by a fraction of a coarsened cell.
+
+  ``stride`` is deprecated in favor of the scalar ``coarsen`` and emits a
+  ``DeprecationWarning``.  Per-axis striding is no longer supported and raises
+  ``ValueError``, as does passing both ``stride`` and a conflicting ``coarsen``.
+
+  .. warning::
+
+     This changes NetCDF behavior when no alignment is requested.  ``align=None``
+     (the default) starts the coarsened grid at the requested ``crop_extent``,
+     matching ASCII, ``crop()`` and the old ``read_netcdf``.  The old ``stride``
+     instead snapped the grid to the file's first cell.  A script that relied on
+     that implicit snap can now get a grid offset by up to ``coarsen - 1``
+     native cells; pass ``align`` to pin the lattice explicitly.
+
+  Passing ``align=[x, y]`` locks the coarsened grid to that lattice regardless
+  of where the file's first cell falls.  This matters for DEM tiles that carry a
+  few cells outside their nominal edge -- some NCEI CUDEM tiles extend slightly
+  beyond the quarter-degree square they cover -- where coarsening should follow
+  the quarter-degree lattice rather than an unknown offset::
+
+      topo = Topography()
+      topo.read('tile.nc', topo_type=4, crop_extent=[-124.2, -124.0, 47.1, 47.3],
+                coarsen=3, align=[-124.25, 47.0])
+
+  If the file's own grid does not contain the requested lattice point, alignment
+  selects the nearest native point instead of landing on the lattice exactly, so
+  the coarsened grid can sit up to half a native cell off.  It is still
+  deterministic and independent of ``crop_extent``: every crop of a given file
+  coarsens onto the same grid, which is the property that makes neighboring
+  tiles line up.  See :ref:`topotools` and :ref:`setrun_topo_preprocessing`.
+
+- **Remote NetCDF URL fix.**
+  :class:`~clawpack.geoclaw.netcdf_utils.NetCDFInspector` no longer passes
+  remote ``scheme://`` URLs through ``pathlib.Path``, which collapsed
+  ``https://`` to ``https:/`` and prepended the working directory, breaking
+  OPeNDAP/THREDDS reads.
+
 - **Python-owned priority ordering.**
   Topography files in ``topo.data`` are now sorted entirely in Python by
   :meth:`~clawpack.geoclaw.data.TopographyData._compute_priority_order`
@@ -126,7 +191,7 @@ Changes to geoclaw
       topo.write('new.tt2', topo_type=2)
 
 - **New** ``topo.data`` **format.**
-  Each per-file block in ``topo.data`` now contains 9 lines (up from 2),
+  Each per-file block in ``topo.data`` now contains 10 lines (up from 2),
   recording all preprocessing attributes.  See :ref:`topodata_format` for
   the complete format specification.
 

doc/dtopo.rst

@@ -2,14 +2,14 @@
 
 .. _dtopo:
 
-.. warning ::  Many changes are being implemented in the way topo and dtopo
-  files are handled, in both the Python tools and the Fortran code.
-  See :ref:`topochanges` for a summary.
-
 *****************************************************************
 Moving topography displacement files
 *****************************************************************
 
+.. warning ::  Many changes are being implemented in the way topo and dtopo
+  files are handled, in both the Python tools and the Fortran code.
+  See :ref:`topochanges` for a summary.
+
 .. seealso::
    - :ref:`dtopotools_module`
    - :ref:`topo`

doc/force_dry.rst

@@ -136,7 +136,7 @@ available from the NCEI thredds server:
 .. code:: ipython3
 
     path = 'https://www.ngdc.noaa.gov/thredds/dodsC/regional/puget_sound_13_mhw_2014.nc'
-    topo = topotools.read_netcdf(path, extent=extent)
+    topo = topotools.fetch_remote_topo(path, crop_extent=extent)
 
 Plot the topo we downloaded:
 

doc/marching_front.rst

@@ -257,7 +257,7 @@ available from the NCEI thredds server:
 .. code:: ipython3
 
     path = 'https://www.ngdc.noaa.gov/thredds/dodsC/regional/puget_sound_13_mhw_2014.nc'
-    topo = topotools.read_netcdf(path, extent=extent)
+    topo = topotools.fetch_remote_topo(path, crop_extent=extent)
 
 .. code:: ipython3
 

doc/netcdf.rst

@@ -1,12 +1,12 @@
 .. _netcdf_input:
 
+GeoClaw NetCDF Input System
+===========================
+
 .. warning ::  Many changes are being implemented in the way topo and dtopo
   files are handled, in both the Python tools and the Fortran code.
   See :ref:`topochanges` for a summary.
 
-GeoClaw NetCDF Input System
-===========================
-
 This document covers the NetCDF input pipeline introduced in the
 ``refactor-netcdf-support``
 [PR #701](https://github.com/clawpack/geoclaw/pull/701)
@@ -136,6 +136,29 @@ If your NetCDF file covers a larger area than your simulation domain
 Only the subset is read into memory at runtime. The full file is never
 loaded.
 
+.. note::
+
+   ``crop_bounds`` and ``crop_extent`` are two different cropping mechanisms
+   and it is worth knowing which one you want.
+
+   - ``crop_bounds`` (this section) is a *descriptor* crop, expressed in the
+     file's own coordinates.  It is recorded in ``topo.data`` and applied by
+     the Fortran code when it reads the NetCDF file at run time.  The file is
+     registered as-is; nothing is read or rewritten in Python.
+   - ``crop_extent`` (see :ref:`setrun_topo_preprocessing`) is a *preprocessing*
+     crop, expressed in domain coordinates.  It is applied by
+     :meth:`~clawpack.geoclaw.topotools.Topography.read` in Python, together
+     with the other preprocessing attributes ``coarsen``, ``buffer``,
+     ``align`` and the shifts.
+
+   Use ``crop_bounds`` when you simply want Fortran to read a window of a large
+   NetCDF file.  Use ``crop_extent`` when you want a
+   :class:`~clawpack.geoclaw.topotools.Topography` object in Python — to
+   coarsen it, shift it, plot it, or write it back out.  To fetch and crop a
+   remote DEM in one call, see
+   :func:`~clawpack.geoclaw.topotools.fetch_remote_topo` and
+   :ref:`noaa_thredds`.
+
 Checking CF compliance
 ^^^^^^^^^^^^^^^^^^^^^^
 

doc/setrun_geoclaw.rst

@@ -207,15 +207,16 @@ Set them on the object before appending to ``topofiles``:
    * - ``coarsen``
      - int
      - ``1``
-     - Stride-subsampling factor (1 = no coarsen).
+     - Subsampling factor (1 = no coarsening).
    * - ``buffer``
      - int
      - ``0``
      - Grid-point margin outside crop region.
    * - ``align``
      - tuple or None
      - ``None``
-     - ``(x, y)`` alignment for coarsened grids.
+     - ``(x, y)`` lattice to pin the coarsened grid to.  ``None`` starts
+       the coarsened grid at ``crop_extent``.
    * - ``x_shift``
      - float
      - ``0.0``
@@ -239,7 +240,9 @@ Set them on the object before appending to ``topofiles``:
 
    - ``crop_extent`` is a preprocessing input. It is distinct from
      ``Topography.extent``, which is a read-only property returning the
-     spatial bounds of the already-loaded data.
+     spatial bounds of the already-loaded data.  It is also distinct from the
+     NetCDF descriptor's ``crop_bounds``, which is applied by Fortran at run
+     time rather than by Python at read time; see :ref:`netcdf_input`.
 
    - ``buffer`` is in grid points, not degrees.  Float values are truncated
      via ``int()``, so ``buffer=0.5`` is equivalent to ``buffer=0``.
@@ -248,8 +251,13 @@ Set them on the object before appending to ``topofiles``:
      If both ``topo_type < 0`` **and** ``negate_z=True`` are active, two
      sign flips are applied and the result is the original Z (net identity).
 
-   - ``coarsen`` uses stride subsampling (every nth point), not cell
-     averaging.
+   - ``coarsen`` subsamples: it keeps every nth point, and does not average
+     over the coarsened cell.
+
+   - ``coarsen`` and ``align`` behave identically for ASCII and NetCDF files,
+     so the same settings applied to the same data in either format give the
+     same grid.  The older NetCDF-only ``stride`` argument to ``read()`` is
+     deprecated in favor of ``coarsen``; see :ref:`topotools_deprecated_args`.
 
 Example using several attributes::
 

doc/topo.rst

@@ -223,23 +223,27 @@ The `NOAA THREDDS server
 can be used to access a variety of topography data sets, including the etopo1
 global data set at 1 arcminute resolution and the etopo2 global data set at 2
 arcminute resolution.  These are available in netCDF format and can be read
-directly into GeoClaw.  As a convenience, you can use the `topotools.read_netcdf
-<topotools_module.html#clawpack.geoclaw.topotools.read_netcdf>`_ function.  Note
-that this also allows reading in only a subset of the data, both limiting the
-extent and the resolution, e.g. by sampling every other point (by setting
-`coarsen=2`). This is particularly useful if you only want a subset of a huge
+directly into GeoClaw.  As a convenience, you can use the
+`topotools.fetch_remote_topo
+<topotools_module.html#clawpack.geoclaw.topotools.fetch_remote_topo>`_ function
+(the older `topotools.read_netcdf` function is now deprecated in favor of it).
+Note that this also allows reading in only a subset of the data, both limiting
+the extent (via `crop_extent`) and the resolution, e.g. by sampling every
+other point (by setting `coarsen=2`). This is particularly useful if you only
+want a subset of a huge
 online netCDF file (e.g. coastal DEMs at 1/3 arcsecond resolution are typically
 several gigabytes).  See :ref:`netcdf_input` for more details on working with
 netCDF files.
 
 The dictionary `topotools.remote_topo_urls` contains some useful URLs for
 etopo and a few other NOAA THREDDS datasets. This allows reading etopo
-60 arc-second data, for example, via::
+30 arc-second data, for example, via::
 
     from clawpack.geoclaw import topotools
-    extent = [-135, -120, 38, 52]
-    topo = topotools.read_netcdf('etopo22_60s', extent=extent,
-                                 coarsen=1, verbose=True)
+    crop_extent = [-135, -120, 38, 52]
+    topo = topotools.fetch_remote_topo('etopo22_30sec',
+                                       crop_extent=crop_extent,
+                                       coarsen=1, verbose=True)
 
 A quick plot of the topography can then be created using::
 
@@ -253,14 +257,16 @@ and the topo can be saved as an ASCII raster topofile via, e.g.::
                grid_registration='llcenter', Z_format='%.1f')
 
 
-See `$CLAW/geoclaw/tests/test_etopo1.py` for one example, in which a very
-small patch from the global etopo1 database (which has 1 arcminute resolution)
-is downloaded at different resolutions.
+See the ETOPO integration tests in `$CLAW/geoclaw/tests/test_topotools.py`
+for an example, in which a very small patch from the global etopo1 database
+(which has 1 arcminute resolution) is downloaded at different resolutions.
 
-**Note:** Earlier versions of clawpack included `etopotools.py` providing a
-different way to download subsampled etopo1 topography.  That has been
-deprecated since the old way is no longer supported by NOAA and did not
-always do the subsampling properly.
+**Note:** The `etopotools.py` module provides `etopotools.fetch_etopo`, a thin
+convenience wrapper around `fetch_remote_topo` for the etopo netCDF datasets,
+along with a legacy `etopotools.etopo1_download` that fetches subsampled etopo1
+data from an older NOAA WCS endpoint.  That WCS endpoint is no longer reliably
+supported by NOAA, so reading the netCDF data directly via `fetch_remote_topo`
+/ `fetch_etopo` is preferred.
 
 **Note:** Data in the NOAA THREDDS server is referenced to NAVD88, not to MHW!
 

doc/topochanges.rst

@@ -91,9 +91,9 @@ Copied from :ref:`changes_to_master`
 
 
 - **Topography preprocessing attributes.**
-  :class:`~clawpack.geoclaw.topotools.Topography` now supports seven
+  :class:`~clawpack.geoclaw.topotools.Topography` now supports eight
   preprocessing attributes (``crop_extent``, ``coarsen``, ``buffer``,
-  ``align``, ``x_shift``, ``z_shift``, ``negate_z``) that are applied
+  ``align``, ``x_shift``, ``y_shift``, ``z_shift``, ``negate_z``) that are applied
   automatically when :meth:`~clawpack.geoclaw.topotools.Topography.read`
   loads a file.  See :ref:`setrun_topo_preprocessing` for the full table
   and :ref:`topotools` for usage examples and operation order.
@@ -109,6 +109,71 @@ Copied from :ref:`changes_to_master`
   are available without loading the elevation array.
   See :ref:`topotools` for an example.
 
+- **Remote-DEM fetch helper;** ``read_netcdf`` **deprecated.**
+  :func:`~clawpack.geoclaw.topotools.fetch_remote_topo` resolves a nickname
+  (a key of ``topotools.remote_topo_urls``), URL, or local path and reads it
+  through the ``topo_type=4`` path, returning a
+  :class:`~clawpack.geoclaw.topotools.Topography` with the requested
+  ``crop_extent``/``coarsen``/``buffer``/``align`` applied.  The older
+  :func:`~clawpack.geoclaw.topotools.read_netcdf` is now a thin shim over it
+  that emits a ``DeprecationWarning``.  In ``etopotools``,
+  :func:`~clawpack.geoclaw.etopotools.fetch_etopo` is a convenience wrapper for
+  the etopo netCDF datasets, and
+  :func:`~clawpack.geoclaw.etopotools.etopo1_download` now always returns a
+  ``Topography``.  Because ``fetch_remote_topo`` uses the ``topo_type=4``
+  reader, the elevation variable must be in meters (or supply ``assume_units``
+  via ``nc_params``).  See :ref:`topo`.
+
+- **Unified coarsening and alignment;** ``stride`` **deprecated.**
+  :meth:`~clawpack.geoclaw.topotools.Topography.read` now accepts
+  ``crop_extent``, ``coarsen``, ``buffer`` and ``align`` directly as keyword
+  arguments; passing one is equivalent to setting the attribute of the same
+  name before calling ``read()``.  ASCII (``topo_type`` 2 or 3) and NetCDF
+  (``topo_type=4``) reads of the same data now return identical ``x``, ``y``
+  and ``Z`` for the same ``coarsen`` and ``align``: both drive off a single
+  shared index computation, so the type-4 lazy hyperslab (including files that
+  store latitude north-to-south) and
+  :meth:`~clawpack.geoclaw.topotools.Topography.crop` agree by construction.
+  Previously the NetCDF-only ``stride`` argument silently did nothing for ASCII
+  reads and used a different alignment convention, so the two paths could
+  disagree by a fraction of a coarsened cell.
+
+  ``stride`` is deprecated in favor of the scalar ``coarsen`` and emits a
+  ``DeprecationWarning``.  Per-axis striding is no longer supported and raises
+  ``ValueError``, as does passing both ``stride`` and a conflicting ``coarsen``.
+
+  .. warning::
+
+     This changes NetCDF behavior when no alignment is requested.  ``align=None``
+     (the default) starts the coarsened grid at the requested ``crop_extent``,
+     matching ASCII, ``crop()`` and the old ``read_netcdf``.  The old ``stride``
+     instead snapped the grid to the file's first cell.  A script that relied on
+     that implicit snap can now get a grid offset by up to ``coarsen - 1``
+     native cells; pass ``align`` to pin the lattice explicitly.
+
+  Passing ``align=[x, y]`` locks the coarsened grid to that lattice regardless
+  of where the file's first cell falls.  This matters for DEM tiles that carry a
+  few cells outside their nominal edge -- some NCEI CUDEM tiles extend slightly
+  beyond the quarter-degree square they cover -- where coarsening should follow
+  the quarter-degree lattice rather than an unknown offset::
+
+      topo = Topography()
+      topo.read('tile.nc', topo_type=4, crop_extent=[-124.2, -124.0, 47.1, 47.3],
+                coarsen=3, align=[-124.25, 47.0])
+
+  If the file's own grid does not contain the requested lattice point, alignment
+  selects the nearest native point instead of landing on the lattice exactly, so
+  the coarsened grid can sit up to half a native cell off.  It is still
+  deterministic and independent of ``crop_extent``: every crop of a given file
+  coarsens onto the same grid, which is the property that makes neighboring
+  tiles line up.  See :ref:`topotools` and :ref:`setrun_topo_preprocessing`.
+
+- **Remote NetCDF URL fix.**
+  :class:`~clawpack.geoclaw.netcdf_utils.NetCDFInspector` no longer passes
+  remote ``scheme://`` URLs through ``pathlib.Path``, which collapsed
+  ``https://`` to ``https:/`` and prepended the working directory, breaking
+  OPeNDAP/THREDDS reads.
+
 - **Python-owned priority ordering.**
   Topography files in ``topo.data`` are now sorted entirely in Python by
   :meth:`~clawpack.geoclaw.data.TopographyData._compute_priority_order`
@@ -129,7 +194,7 @@ Copied from :ref:`changes_to_master`
       topo.write('new.tt2', topo_type=2)
 
 - **New** ``topo.data`` **format.**
-  Each per-file block in ``topo.data`` now contains 9 lines (up from 2),
+  Each per-file block in ``topo.data`` now contains 10 lines (up from 2),
   recording all preprocessing attributes.  See :ref:`topodata_format` for
   the complete format specification.