diff --git a/doc/changes_to_master.rst b/doc/changes_to_master.rst index a63f05e..91fe8f0 100644 --- a/doc/changes_to_master.rst +++ b/doc/changes_to_master.rst @@ -84,6 +84,9 @@ See `amrclaw diffs Changes to geoclaw ------------------ +- **Overview of changes related to topo and dtopo files** + See :ref:`topochanges`. + - **Topography preprocessing attributes.** :class:`~clawpack.geoclaw.topotools.Topography` now supports seven preprocessing attributes (``crop_extent``, ``coarsen``, ``buffer``, @@ -118,9 +121,9 @@ Changes to geoclaw now emit a ``DeprecationWarning``. Setting any preprocessing attribute before reading a type-1 file raises ``NotImplementedError``. To convert:: - t = Topography() - t.read('old.tt1', topo_type=1) # DeprecationWarning - t.write('new.tt2', topo_type=2) + topo = Topography() + topo.read('old.tt1', topo_type=1) # gives DeprecationWarning + 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), diff --git a/doc/community.rst b/doc/community.rst index d42c8f1..b6da0be 100644 --- a/doc/community.rst +++ b/doc/community.rst @@ -77,11 +77,15 @@ Available for streaming Upcoming ^^^^^^^^ -No events currently planned. +- None planned. Recent ^^^^^^^^ +- GeoClaw Users' Workshop on Zoom, March 19-20, 2026. + For more information and the registration link, see + `__ + - Tutorial on GeoClaw and D-Claw at the `CSDMS Annual Meeting `__ in May, 2024. diff --git a/doc/dtopo.rst b/doc/dtopo.rst index 7ad42d3..20538ef 100644 --- a/doc/dtopo.rst +++ b/doc/dtopo.rst @@ -2,6 +2,10 @@ .. _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 ***************************************************************** diff --git a/doc/geoclaw.rst b/doc/geoclaw.rst index e37368c..d3bc208 100644 --- a/doc/geoclaw.rst +++ b/doc/geoclaw.rst @@ -46,10 +46,12 @@ More will eventually appear in the :ref:`apps`. geoclaw_started geohints + topochanges topo grid_registration topo_order topodata_format + netcdf topotools dtopo kmltools_module diff --git a/doc/netcdf.rst b/doc/netcdf.rst index d82f35d..c44d5e5 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -1,10 +1,16 @@ .. _netcdf_input: +.. 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 (first released after v5.14.0). It has two sections: a user +``refactor-netcdf-support`` +[PR #701](https://github.com/clawpack/geoclaw/pull/701) +(merged but not yet released). It has two sections: a user guide for scientists who want to use NetCDF files as input, and a developer reference for those working on the implementation. diff --git a/doc/topo.rst b/doc/topo.rst index 22a15b8..a0e98f5 100644 --- a/doc/topo.rst +++ b/doc/topo.rst @@ -6,6 +6,10 @@ Topography data ***************************************************************** +.. 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:`topotools` - :ref:`grid_registration` diff --git a/doc/topo_order.rst b/doc/topo_order.rst index 311394f..990f9d2 100644 --- a/doc/topo_order.rst +++ b/doc/topo_order.rst @@ -6,6 +6,10 @@ Topography file ordering ***************************************************************** +.. 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. + .. warning :: This feature has bugs in v5.13.0 and v5.13.1 and should not be used in those versions. These bugs were fixed (we hope) in v5.14.0. See :ref:`release_5_14_0` for more discussion. diff --git a/doc/topochanges.rst b/doc/topochanges.rst new file mode 100644 index 0000000..f1d26f0 --- /dev/null +++ b/doc/topochanges.rst @@ -0,0 +1,155 @@ +.. _topochanges: + +***************************************************************** +Changes to topo and dtopo handling (planned for v5.15.0) +***************************************************************** + +Many changes are being implemented in the way topo and dtopo +files are handled, in both the Python tools and the Fortran code. + +.. warning :: Some of this work is still being debugged and code + on the master branch may not work as advertised. Interfaces + and parameter names are still evolving. Please test this code + but be aware of these warnings. + +.. seealso:: + - :ref:`topotools` + - :ref:`topodata_format` + - :ref:`netcdf_input` + - :ref:`topotools` + - :ref:`topo_order` + - :ref:`dtopo` + - :ref:`dtopotools` + +Major changes to topo file handling +=================================== + +Cropping and coarsening or minor adjustments to a topo DEM +---------------------------------------------------------- + +In the Fortran GeoClaw code, rather than specifying a list of +topo files in `setrun.py`, we now list a set of topo grids. +Each topography grid comes from a particular file, but has +additional parameters indicating if the DEM in the file should +be cropped or coarsened to create the topo grid (along with +other options, such as shifting the longitude by 360 degrees +or shifting the values vertically). + +This has several advantages: + +- A single large topo DEM can be dynamically cropped and coarsened + as needed, rather than potentially needing to make many smaller + files with different extents or resolutions. The same DEM can + be used multiple times in the `setrun.py` to specify different + topo grids with different extents and coarsening factors, if + desired. + +- If the same topo file is needed sometimes in latitude W (e.g. + in modeling a nearfield tsunami on the US West Coast) and sometimes + in longitude E (e.g. when modeling a farfield tsunami from Japan), + you no longer need two versions of the large file with only one + change in the header. + +See :ref:`setrun_topo_preprocessing` for a full attribute table with types +and defaults, and non-obvious behavior notes. + +NetCDF file handling +-------------------- + +Improvements have also been made in how netCDF files are used +in both Python and Fortran. See :ref:`netcdf_input`. + + + +Ordering of topo DEM priorities +------------------------------- + +When creating grid cell topo values, normally the topo grids specified +are used with the finest available grid having highest priority. +However, this can be over-ridden as is sometime necessary as described +in :ref:`topo_order`. Previously the ordering was figured out in the +Fortran code. Now the ordering is figured out in Python, in the process +of generating `topo.data` from `setrun.py`, and the order the grids +are listed in `topo.data` is exactly the order that will be used in +the Fortran code. This should make it easier for the user to confirm +the desired ordering is being used. + + +Major changes to dtopo file handling +==================================== + +NetCDF file handling +-------------------- + +Improvements have also been made in how netCDF files are used +in both Python and Fortran. See :ref:`netcdf_input`. + + + +Copied from :ref:`changes_to_master` +==================================== + + +- **Topography preprocessing attributes.** + :class:`~clawpack.geoclaw.topotools.Topography` now supports seven + preprocessing attributes (``crop_extent``, ``coarsen``, ``buffer``, + ``align``, ``x_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. + +- **CF-aware NetCDF reading.** + NetCDF topography files (``topo_type=4``) are now read via + :class:`~clawpack.geoclaw.netcdf_utils.TopoInspector`, which auto-detects + coordinate variable names using CF conventions (``standard_name``, ``axis``, + and common fallback names). Files with non-standard coordinate names and + non-standard dimension orders are handled automatically. + :meth:`~clawpack.geoclaw.topotools.Topography.read_header` also uses CF + detection for type-4 files, enabling a lazy-load pattern where coordinates + are available without loading the elevation array. + See :ref:`topotools` for an example. + +- **Python-owned priority ordering.** + Topography files in ``topo.data`` are now sorted entirely in Python by + :meth:`~clawpack.geoclaw.data.TopographyData._compute_priority_order` + before writing. Files are written coarsest-first (finest last), matching + the traditional GeoClaw listing order; the last file listed in ``topo.data`` + is assigned rank 1 (highest priority) by Fortran, with no Fortran-side + sorting. ``rundata.topo_data.override_order = True`` preserves + user-specified list order; when used, the finest (highest-resolution) file + should be listed last. See :ref:`topo_order`. + +- **topo_type=1 deprecated.** + Reading and writing ``topo_type=1`` (``x y z`` one-point-per-line ASCII) + now emit a ``DeprecationWarning``. Setting any preprocessing attribute + before reading a type-1 file raises ``NotImplementedError``. To convert:: + + topo = Topography() + topo.read('old.tt1', topo_type=1) # gives DeprecationWarning + 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), + recording all preprocessing attributes. See :ref:`topodata_format` for + the complete format specification. + +- **dtopo NetCDF (**\ ``dtopo_type=4``\ **).** + :class:`~clawpack.geoclaw.dtopotools.DTopography` reads and writes + CF-compliant NetCDF dtopo files. The optional ``time_reference`` attribute + selects a CF datetime time axis (``units = "seconds since "``) that a + plain ``xarray.open_dataset`` decodes to ``datetime64``; without it a bare + ``units = "seconds"`` (simulation-relative) axis is written. The reader + scales the time axis by its CF ``units``, fixing an earlier bug where a + ``"hours"``/``"minutes"`` axis was misread as seconds. See + :ref:`netcdf_input` and :ref:`dtopo_formats`. + +- **NetCDF write dtype control.** + ``Topography.write(topo_type=4, z_dtype=...)`` and + ``DTopography.write(dtopo_type=4, dz_dtype=...)`` override the on-disk + elevation/deformation dtype (default ``"float32"``; pass ``"float64"`` for + full double precision). + +- **Robust NetCDF file opening.** + A NetCDF backend engine is now selected explicitly, so a valid file opens + even when its name uses a non-standard extension (e.g. ``.dtt3``) that + xarray's extension-based engine guessing would not recognize. diff --git a/doc/topodata_format.rst b/doc/topodata_format.rst index e849578..2396180 100644 --- a/doc/topodata_format.rst +++ b/doc/topodata_format.rst @@ -1,5 +1,9 @@ .. _topodata_format: +.. 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. + topo.data File Format ===================== diff --git a/doc/topotools.rst b/doc/topotools.rst index 1bff5a1..22263af 100644 --- a/doc/topotools.rst +++ b/doc/topotools.rst @@ -4,144 +4,13 @@ Python tools for working with topo and dtopo -------------------------------------------- -.. seealso:: - - :ref:`topo` - - :ref:`topo_order` - - :ref:`topodata_format` - - -Preprocessing attributes -~~~~~~~~~~~~~~~~~~~~~~~~~ - -:class:`~clawpack.geoclaw.topotools.Topography` objects support seven -preprocessing attributes that are applied automatically by -:meth:`~clawpack.geoclaw.topotools.Topography.read` in this order: - -1. ``negate_z`` — flip sign of Z (independent of ``topo_type < 0``). -2. ``z_shift`` — add a constant to all non-missing Z values. -3. ``x_shift`` — add a constant to all x coordinates. -4. ``crop_extent``, ``buffer``, ``align``, ``coarsen`` — crop and subsample - via :meth:`~clawpack.geoclaw.topotools.Topography.crop`. - -Set them before calling ``read()``:: - - from clawpack.geoclaw.topotools import Topography - - t = Topography() - t.crop_extent = [-100., -60., 10., 50.] - t.coarsen = 2 - t.z_shift = 10.0 - t.read('bathymetry.nc', topo_type=4) - -See :ref:`setrun_topo_preprocessing` for a full attribute table with types -and defaults, and non-obvious behavior notes. - - -Lazy-load pattern for NetCDF (read_header) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -For NetCDF files (``topo_type=4``), -:meth:`~clawpack.geoclaw.topotools.Topography.read_header` uses -:class:`~clawpack.geoclaw.netcdf_utils.TopoInspector` to detect coordinate -variable names via CF conventions (``standard_name``, ``axis``, or common -fallback names such as ``lon``/``lat``, ``x``/``y``). After calling -``read_header()``, ``extent`` and ``delta`` are populated without loading the -elevation array. Accessing ``Z`` triggers a deferred ``read()``: - -.. code-block:: python - - t = Topography() - t.path = 'large_dem.nc' - t.topo_type = 4 - t.read_header() # fast: reads coordinate arrays only - print(t.extent) # available immediately - print(t.delta) # available immediately - z = t.Z # triggers full read on first access - -This pattern is also used internally by ``TopographyData._compute_priority_order`` -to determine file resolution without loading elevation data. - - -Vertical datum metadata -~~~~~~~~~~~~~~~~~~~~~~~~~ - -:class:`~clawpack.geoclaw.topotools.Topography` and -:class:`~clawpack.geoclaw.dtopotools.DTopography` carry an optional ``datum`` -attribute recording the vertical reference of the elevation/deformation data -(e.g. ``'MSL'`` or ``'NAVD88'``). It defaults to ``None``, is populated from a -NetCDF ``vertical_datum`` (or ``datum``) attribute on read, and is written -back out on NetCDF write (``topo_type=4`` / ``dtopo_type=4``). ASCII formats -have no datum field, so the value is not persisted for types 1/2/3. +.. 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. -The datum is **informational only** -- GeoClaw performs no vertical-datum -transformation and always uses the Z values as given. As a guard, both -:meth:`~clawpack.geoclaw.data.TopographyData.write` (when producing -``topo.data``) and :meth:`~clawpack.geoclaw.data.DTopoData.write` (``dtopo.data``) -issue a warning if the files they are given carry more than one distinct -datum, since mixing vertical references without converting between them is a -likely error. -.. deprecated:: - ``topo_type=1`` (three-column ``x y z`` ASCII, one point per line) is - deprecated. Reading emits a ``DeprecationWarning``; writing also emits a - ``DeprecationWarning``; setting any preprocessing attribute before reading - raises ``NotImplementedError``. - - To convert a type-1 file:: - - t = Topography() - t.read('old.tt1', topo_type=1) # DeprecationWarning - t.write('new.tt2', topo_type=2) # save as type 2 - - Genuinely unstructured (scattered) point data cannot be converted this - way. Either grid it externally (e.g. ``scipy.interpolate``, GMT) or use - :meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` (see - `Gridding unstructured (scattered) data`_ below). - - -Gridding unstructured (scattered) data -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -GeoClaw's solvers require topography on a regular, logically rectangular grid -(Cartesian x-y in meters, or lon-lat -- see ``coordinate_system`` in -:ref:`setrun_geoclaw`), but survey or sounding data often comes as scattered -``(x, y, z)`` points. A -:class:`~clawpack.geoclaw.topotools.Topography` constructed with -``unstructured=True`` holds such points in its ``x``, ``y``, and ``z`` -arrays, and -:meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` -interpolates them onto a regular grid, optionally filling gaps from one or -more coarser (structured or unstructured) "fill" topographies: - -.. code-block:: python - - from clawpack.geoclaw.topotools import Topography - - # Scattered survey points (x, y, z) - survey = Topography(unstructured=True) - survey.x = x_points - survey.y = y_points - survey.z = z_values - - # A coarser regional DEM used to fill gaps between survey points - regional = Topography(path='regional.tt3', topo_type=3) - regional.read() - - survey.interp_unstructured(regional, extent=[x1, x2, y1, y2], - proximity_radius=100.0) - survey.unstructured # now False -- it holds a regular grid - survey.write('combined.tt3', topo_type=3) - -The grid spacing is taken from the minimum spacing of the scattered points -(bounded below by ``delta_limit`` meters) unless set explicitly with -``delta``. A fill point is used only where it lies inside ``extent``, carries -valid data, and is more than ``proximity_radius`` meters from every scattered -point, so the higher-resolution survey data is preferred wherever it exists. -Missing fill values (``NaN`` in memory, or a numeric ``no_data_value``) are -dropped. See -:meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` for the -full set of parameters. +The following pages describe the Python tools in more detail: .. toctree:: @@ -151,3 +20,12 @@ full set of parameters. dtopotools_module geoclaw_util_module kmltools_module + + +.. seealso:: + - :ref:`topo` + - :ref:`dtopo` + - :ref:`netcdf_input` + + + diff --git a/doc/topotools_module.rst b/doc/topotools_module.rst index fbb9b74..822eca6 100644 --- a/doc/topotools_module.rst +++ b/doc/topotools_module.rst @@ -4,18 +4,157 @@ topotools module for working with topography data ================================================= +.. 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:`topo` - - :ref:`topotools` + - :ref:`netcdf_input` + + +Preprocessing attributes +~~~~~~~~~~~~~~~~~~~~~~~~~ + +:class:`~clawpack.geoclaw.topotools.Topography` objects support seven +preprocessing attributes that are applied automatically by +:meth:`~clawpack.geoclaw.topotools.Topography.read` in this order: + +1. ``negate_z`` — flip sign of Z (independent of ``topo_type < 0``). +2. ``z_shift`` — add a constant to all non-missing Z values. +3. ``x_shift`` — add a constant to all x coordinates. +4. ``crop_extent``, ``buffer``, ``align``, ``coarsen`` — crop and subsample + via :meth:`~clawpack.geoclaw.topotools.Topography.crop`. + +Set them before calling ``read()``:: + + from clawpack.geoclaw.topotools import Topography + + topo = Topography() + topo.crop_extent = [-100., -60., 10., 50.] + topo.coarsen = 2 + topo.z_shift = 10.0 + topo.read('bathymetry.nc', topo_type=4) + +See :ref:`setrun_topo_preprocessing` for a full attribute table with types +and defaults, and non-obvious behavior notes. + + +Lazy-load pattern for NetCDF (read_header) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For NetCDF files (``topo_type=4``), +:meth:`~clawpack.geoclaw.topotools.Topography.read_header` uses +:class:`~clawpack.geoclaw.netcdf_utils.TopoInspector` to detect coordinate +variable names via CF conventions (``standard_name``, ``axis``, or common +fallback names such as ``lon``/``lat``, ``x``/``y``). After calling +``read_header()``, ``extent`` and ``delta`` are populated without loading the +elevation array. Accessing ``Z`` triggers a deferred ``read()``: + +.. code-block:: python + + topo = Topography() + topo.path = 'large_dem.nc' + topo.topo_type = 4 + topo.read_header() # fast: reads coordinate arrays only + print(topo.extent) # available immediately + print(topo.delta) # available immediately + z = topo.Z # triggers full read on first access + +This pattern is also used internally by ``TopographyData._compute_priority_order`` +to determine file resolution without loading elevation data. + + +Vertical datum metadata +~~~~~~~~~~~~~~~~~~~~~~~~~ + +:class:`~clawpack.geoclaw.topotools.Topography` and +:class:`~clawpack.geoclaw.dtopotools.DTopography` carry an optional ``datum`` +attribute recording the vertical reference of the elevation/deformation data +(e.g. ``'MSL'`` or ``'NAVD88'``). It defaults to ``None``, is populated from a +NetCDF ``vertical_datum`` (or ``datum``) attribute on read, and is written +back out on NetCDF write (``topo_type=4`` / ``dtopo_type=4``). ASCII formats +have no datum field, so the value is not persisted for types 1/2/3. + +The datum is **informational only** -- GeoClaw performs no vertical-datum +transformation and always uses the Z values as given. As a guard, both +:meth:`~clawpack.geoclaw.data.TopographyData.write` (when producing +``topo.data``) and :meth:`~clawpack.geoclaw.data.DTopoData.write` (``dtopo.data``) +issue a warning if the files they are given carry more than one distinct +datum, since mixing vertical references without converting between them is a +likely error. + + +.. deprecated:: + ``topo_type=1`` (three-column ``x y z`` ASCII, one point per line) is + deprecated. Reading emits a ``DeprecationWarning``; writing also emits a + ``DeprecationWarning``; setting any preprocessing attribute before reading + raises ``NotImplementedError``. + + To convert a type-1 file:: + + topo = Topography() + topo.read('old.tt1', topo_type=1) # DeprecationWarning + topo.write('new.tt2', topo_type=2) # save as type 2 + + Genuinely unstructured (scattered) point data cannot be converted this + way. Either grid it externally (e.g. ``scipy.interpolate``, GMT) or use + :meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` (see + `Gridding unstructured (scattered) data`_ below). + + +Gridding unstructured (scattered) data +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +GeoClaw's solvers require topography on a regular, logically rectangular grid +(Cartesian x-y in meters, or lon-lat -- see ``coordinate_system`` in +:ref:`setrun_geoclaw`), but survey or sounding data often comes as scattered +``(x, y, z)`` points. A +:class:`~clawpack.geoclaw.topotools.Topography` constructed with +``unstructured=True`` holds such points in its ``x``, ``y``, and ``z`` +arrays, and +:meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` +interpolates them onto a regular grid, optionally filling gaps from one or +more coarser (structured or unstructured) "fill" topographies: + +.. code-block:: python + + from clawpack.geoclaw.topotools import Topography + + # Scattered survey points (x, y, z) + survey = Topography(unstructured=True) + survey.x = x_points + survey.y = y_points + survey.z = z_values + + # A coarser regional DEM used to fill gaps between survey points + regional = Topography(path='regional.tt3', topo_type=3) + regional.read() + + survey.interp_unstructured(regional, extent=[x1, x2, y1, y2], + proximity_radius=100.0) + survey.unstructured # now False -- it holds a regular grid + survey.write('combined.tt3', topo_type=3) + +The grid spacing is taken from the minimum spacing of the scattered points +(bounded below by ``delta_limit`` meters) unless set explicitly with +``delta``. A fill point is used only where it lies inside ``extent``, carries +valid data, and is more than ``proximity_radius`` meters from every scattered +point, so the higher-resolution survey data is preferred wherever it exists. +Missing fill values (``NaN`` in memory, or a numeric ``no_data_value``) are +dropped. See +:meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` for the +full set of parameters. The notebook `topotools_examples `__ -illustrates how to use some of the tools. +illustrates how to use some of the tools. **Needs to be updated.** + The file `$CLAW/geoclaw/tests/test_topotools.py` contains some tests of these tools. Looking at these test routines may also give some ideas on -how to use them. +how to use them. **Needs to be updated.** Documentation auto-generated from the module docstrings