From 1ae9037a747fe7d4cdbdb3e59b1f1109fd5cc0cb Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Thu, 16 Apr 2026 09:24:46 -0400 Subject: [PATCH 01/16] Repurpose netcdf.rst for input only description Adds both user and developer level documentation. Need to still fix references to this document. --- doc/netcdf.rst | 411 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 403 insertions(+), 8 deletions(-) diff --git a/doc/netcdf.rst b/doc/netcdf.rst index fdcb9af..bfd16d6 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -1,12 +1,407 @@ -:orphan: +.. toctree:: + :maxdepth: 2 -.. _netcdf: + user_guide + developer_reference + descriptor_format + unit_contract + coordinate_normalization + technical_debt + adding_met_field + test_structure + next_steps -.. seealso:: :ref:`topo_netcdf`. +GeoClaw NetCDF Input System +=========================== -========================== -Using NetCDF output -========================== +This document covers the NetCDF input pipeline introduced in the +``refactor-netcdf-support`` PR. 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. -NetCDF output is not currently supported in Clawpack. This is not a suitable -format for AMR style data. +-------------- + +User Guide +---------- + +Topography from NetCDF +~~~~~~~~~~~~~~~~~~~~~~ + +GeoClaw can read topography/bathymetry directly from a NetCDF file +(``topo_type=4``). The file must contain a 2D elevation variable on a +regular latitude/longitude grid. Common sources include GEBCO, ETOPO, +NOAA coastal DEMs, and any file produced by standard GIS tools. + +What GeoClaw handles automatically +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- Longitude conventions: both [-180, 180] and [0, 360] are detected and + normalized at runtime. You do not need to preprocess the file. +- Latitude ordering: both S-to-N and N-to-S are handled correctly. +- Dimension ordering: ``(lat, lon)`` and ``(lon, lat)`` are both + supported. +- Fill values: ``_FillValue`` and ``missing_value`` attributes are + resolved automatically. If a fill value is found within your + simulation domain, GeoClaw will abort with an error -- this is + intentional, as missing bathymetry is a silent correctness hazard. + +What your file must provide +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- A 2D variable containing elevation in **meters** (positive up, + negative for ocean). If your file uses feet or another unit, use + ``CFNormalizer`` to convert before running (see below). +- 1D coordinate variables for latitude and longitude with recognizable + names (``lat``/``latitude``/``y`` and ``lon``/``longitude``/``x`` are + all detected). Curvilinear (2D coordinate) grids are not currently + supported for topography. + +Registering a NetCDF topo file in setrun.py +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code:: python + + from clawpack.geoclaw.netcdf_utils import TopoInterrogator + meta = TopoInterrogator('bathy.nc', var_name='z', + crop_bounds=(-100, -60, 15, 35)).interrogate_topo() + rundata.topo_data.topofiles.append([4, 'bathy.nc', meta]) + +That is the only change required in most cases. GeoClaw's Python layer +interrogates the file when you run ``setrun.py`` and writes the +necessary descriptor information into ``topo.data`` automatically. +Variable and coordinate names are auto-detected from CF attributes where +possible; pass ``var_name``, ``lon_name``, or ``lat_name`` explicitly to +``TopoInterrogator`` if auto-detection fails. + +Domain subsetting (crop) +^^^^^^^^^^^^^^^^^^^^^^^^ + +If your NetCDF file covers a larger area than your simulation domain +(common with global or regional datasets), you can crop at read time +without creating a smaller file: + +.. code:: python + + topo_data.topofiles.append([4, 'gebco_global.nc', + {'crop_bounds': [-100, -80, 20, 35]}]) + +Only the subset is read into memory at runtime. The full file is never +loaded. + +Checking CF compliance +^^^^^^^^^^^^^^^^^^^^^^ + +If you are unsure whether your file will be read correctly, the +``CFNormalizer`` utility can inspect and repair common issues: + +.. code:: python + + from clawpack.geoclaw.netcdf_utils import CFNormalizer + + cf = CFNormalizer('path/to/bathymetry.nc') + cf.report() # prints any issues found + cf.normalize('path/to/bathymetry_cf.nc') # writes a corrected copy + +``CFNormalizer`` adds missing ``standard_name``, ``axis``, and ``units`` +attributes, and resolves ``_FillValue``/``missing_value`` conflicts. It +does not resample or reproject. + +-------------- + +Storm surge met forcing from NetCDF +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For storm surge simulations using full gridded met forcing (wind and +pressure fields), GeoClaw can read directly from a NetCDF file. This +replaces the need to convert to OWI ASCII format. + +Supported source formats: + +- **ERA5** (ECMWF reanalysis): detected automatically from CF attributes +- **NWS13** (OWI NetCDF): detected automatically +- **Generic CF-compliant NetCDF**: requires variable name mapping if + names are non-standard + +Not yet supported: raw WRF output (requires preprocessing due to +curvilinear grid and string-encoded time axis). + +Required variables and units +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +GeoClaw expects the following in the NetCDF file **after any unit +conversion** (conversion happens automatically during Python +preprocessing): + +============================= ================================ +Variable Unit +============================= ================================ +Wind (u-component, eastward) m/s +Wind (v-component, northward) m/s +Surface pressure Pa +Time seconds from user-defined offset +============================= ================================ + +If your file uses hPa/mbar for pressure or knots for wind, +``MetInterrogator`` will convert automatically. + +Registering a NetCDF storm file +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In ``setrun.py``: + +.. code:: python + + surge_data.storm_specification_type = 'data' + surge_data.storm_file = 'isaac.storm' + +Then create the storm descriptor file, *e.g.* ``isaac.storm``, using the +Python API: + +.. code:: python + + from clawpack.geoclaw.surge.storm import Storm + import numpy as np + + storm = Storm() + storm.time_offset = np.datetime64('2012-08-29') + storm.file_format = 'netcdf' + storm.file_paths = ['path/to/forcing.nc'] + storm.write('isaac.storm', file_format='data') + +If your variable or dimension names are non-standard, provide a mapping: + +.. code:: python + + storm.write('isaac.storm', file_format='data', + dim_mapping={'t': 'valid_time'}, + var_mapping={'wind_u': 'u10', 'wind_v': 'v10', + 'pressure': 'msl'}) + +Time handling +^^^^^^^^^^^^^ + +GeoClaw works in seconds from a user-defined offset (typically landfall +or storm genesis). All CF time decoding -- including calendar handling +and unit conversion from hours/days -- is done in Python. The Fortran +runtime sees only seconds from offset. + +Set the offset when constructing the storm: + +.. code:: python + + storm.time_offset = np.datetime64('2012-08-29T00:00') # landfall time + +-------------- + +Developer Reference +------------------- + +Architecture overview +~~~~~~~~~~~~~~~~~~~~~ + +The system has a strict Python/Fortran split: + +**Python** handles: file interrogation, CF attribute parsing, coordinate +convention detection, fill value resolution, unit conversion, time +decoding, crop bound validation, and descriptor writing. + +**Fortran** handles: opening the NetCDF file at runtime using +information from the descriptor, index arithmetic for coordinate +normalization (no data copies), domain subsetting via +``start``/``count`` arguments to ``nf90_get_var``, and time-slice reads +for met forcing. + +Fortran assumes the unit contract from ``GEOCLAW_NETCDF_UNITS`` in +``units.py`` without checking. Python enforces it. + +Class hierarchy (``netcdf_utils.py``) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +:: + + NetCDFInterrogator + - open file (xarray, Dask-lazy) + - discover coordinate variables by name heuristics + CF standard_name + - detect lon convention, lat order, dim order + - resolve fill value (_FillValue wins over missing_value per CF spec) + - validate crop bounds against file extent + - output: NetCDFDescriptor dataclass + + TopoInterrogator(NetCDFInterrogator) + - detect fill values within crop region (hard error) + - verify and convert units to contract + - no multi-file coverage logic (Fortran handles compositing) + + MetInterrogator(NetCDFInterrogator) + - check wind_u, wind_v, pressure present and on same grid/time axis + - convert units to contract + - decode CF time to seconds from offset + - detect ensemble/member dimensions (hard error if non-singleton) + - partial domain coverage is allowed (Fortran fills edges) + + DescriptorWriter + - topo: writes key=value lines for inline inclusion in topo.data, using a + blank line to terminate the Fortran read loop + - met: writes &file_info namelist block + repeated &variable_info + blocks for *.storm file body + + CFNormalizer + - adds/fixes CF attributes without modifying data + - idempotent + +Descriptor format +~~~~~~~~~~~~~~~~~ + +Topo (lines in topo.data after topo_type) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +:: + + var_name = z + lon_name = longitude + lat_name = latitude + lon_convention = 180 + lat_order = S_to_N + dim_order = lat,lon + fill_value = -9999.0 + fill_action = abort + crop_bounds = -100.0 -80.0 20.0 35.0 + +``fill_action = abort`` is the only supported value for topography. +``crop_bounds`` is omitted if no crop is specified. + +Met forcing (body of \*.storm after format header) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +:: + + # format: netcdf + &file_info + source_file = /path/to/forcing.nc + lon_name = longitude + lat_name = latitude + time_name = valid_time + dim_order = time,lat,lon + lon_convention = 360 + lat_order = S_to_N + fill_value = -9999.0 + fill_action = warn + time_offset = 0.0 + crop_bounds = -100.0 -80.0 20.0 35.0 + / + &variable_info var_name=u10 geoclaw_role=wind_u / + &variable_info var_name=v10 geoclaw_role=wind_v / + &variable_info var_name=msl geoclaw_role=pressure / + +The ``&variable_info`` blocks use a manually parsed key=value format +(not true repeated Fortran namelists) for compiler portability. Fortran +reads them in a loop until EOF. Adding new roles (precipitation, +friction) requires no change to the Fortran parser. + +Unit contract +~~~~~~~~~~~~~ + +Defined in ``units.py``: + +.. code:: python + + GEOCLAW_NETCDF_UNITS = { + "topo": "m", + "wind_u": "m/s", + "wind_v": "m/s", + "pressure": "Pa", + "time": "s", + } + +All conversion happens in Python before the descriptor is written. +Fortran trusts the descriptor and never checks units. If you add a new +``geoclaw_role``, add its contract unit here first. + +Fortran coordinate normalization +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Fortran applies coordinate normalization via index arithmetic only -- no +arrays are duplicated. The descriptor provides ``lon_convention`` and +``lat_order``; the reader uses these to compute correct indices when +calling ``nf90_get_var``. For crop bounds, ``start`` and ``count`` are +computed from the coordinate arrays at runtime. + +All ``nf90_*`` calls go through a checked interface in +``topo_module.f90`` that prints a meaningful diagnostic and stops +cleanly on failure. This is important for SLURM batch jobs where a +silent Fortran ``STOP`` produces no useful output. + +Known technical debt +~~~~~~~~~~~~~~~~~~~~ + +- ``util.get_netcdf_names`` and the coordinate/variable discovery logic + in ``NetCDFInterrogator`` are parallel implementations. They should + be consolidated -- ``NetCDFInterrogator`` should be the single source + and ``util.get_netcdf_names`` should delegate to it. This is deferred + to the surge module refactor. +- ``Storm.write(file_format="data")`` does not yet call + ``DescriptorWriter`` directly for NetCDF storm entries. Currently the + tests write descriptors manually. The integration should happen as + part of the surge refactor. +- WRF raw output requires ``MetPreprocessor`` for string-time decoding + (``Times`` character array, not a numeric CF time variable) and + curvilinear grid handling (``XLAT``/``XLONG`` are 2D time-varying + arrays). A skip-marked test stub documents the gap in + ``test_storm.py``. + +Adding a new met forcing field (e.g. precipitation) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +1. Add the role and unit to ``GEOCLAW_NETCDF_UNITS`` in ``units.py`` +2. Add detection logic to ``MetInterrogator`` for the new variable +3. ``DescriptorWriter`` requires no change -- new ``&variable_info`` + blocks are written automatically for any role the interrogator + returns +4. Add the Fortran reader logic in the storm NetCDF module to consume + the new ``geoclaw_role`` value +5. Add unit tests in ``tests/netcdf/test_met_interrogator.py`` +6. Add a regression test in the appropriate storm surge example + +Test structure +~~~~~~~~~~~~~~ + +:: + + tests/netcdf/ unit tests for netcdf_utils.py + conftest.py in-memory NetCDF fixtures + test_base_interrogator.py + test_topo_interrogator.py + test_met_interrogator.py + test_descriptor_writer.py + test_cf_normalizer.py + + tests/test_storm.py extended for ERA5 and NWS13 variants + + examples/tsunami/bowl-slosh-netcdf/ coordinate variant regression tests + examples/tsunami/chile2010/ topotools write path + coord variants + examples/storm-surge/isaac/ ERA5 and NWS13 met forcing regression + +All test NetCDF files are generated in-memory or in ``tmp_path`` -- no +binary files are committed to the repository. + +-------------- + +Next steps +---------- + +- **Surge module refactor**: consolidate ``util.get_netcdf_names`` into + ``NetCDFInterrogator``, wire ``Storm.write`` to use + ``DescriptorWriter``, clean up parallel discovery paths in + ``storm.py`` +- **WRF support**: implement ``MetPreprocessor`` to handle string-time + axis and curvilinear grid; the skip-marked test stub in + ``test_storm.py`` documents the expected interface +- **Precipitation field**: reserved in ``GEOCLAW_NETCDF_UNITS``, add + ``MetInterrogator`` detection and Fortran consumer +- **Friction field**: same pattern as precipitation +- **CLI tool**: expose ``CFNormalizer`` as a command-line utility for + users who want to check or repair their NetCDF files before running + GeoClaw +- **Duplicate code audit**: other places in the GeoClaw codebase that + do ad-hoc NetCDF variable discovery should be identified and migrated + to ``NetCDFInterrogator`` From d28de7035dbe971e49b43b43824bfdc25a29da04 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Thu, 16 Apr 2026 17:07:51 -0400 Subject: [PATCH 02/16] Add docs for new netcdf input for GeoClaw --- doc/netcdf.rst | 33 +++++++++++++++++++++++++-------- doc/topo.rst | 48 +++++++++++++++++++++++------------------------- 2 files changed, 48 insertions(+), 33 deletions(-) diff --git a/doc/netcdf.rst b/doc/netcdf.rst index bfd16d6..7ef3bb6 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -11,12 +11,14 @@ test_structure next_steps +.. _netcdf_input: + GeoClaw NetCDF Input System =========================== This document covers the NetCDF input pipeline introduced in the -``refactor-netcdf-support`` PR. It has two sections: a user guide for -scientists who want to use NetCDF files as input, and a developer +``refactor-netcdf-support`` PR (VERSION ?5.15?). 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. -------------- @@ -59,6 +61,20 @@ What your file must provide Registering a NetCDF topo file in setrun.py ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +If your file meets the above requirements and either has the variable `z` or +only one 2D variable in the file, you can simply use the following in +``setrun.py``: + +.. code:: python + + rundata.topo_data.topofiles.append([4, 'bathy.nc']) + +This matches with what is expected for backwards compatibility. + +If your file meets the above requirements and has a non-standard variable name +or you want to specify crop bounds you can use the Python API to interrogate the +file and write a descriptor: + .. code:: python from clawpack.geoclaw.netcdf_utils import TopoInterrogator @@ -67,11 +83,12 @@ Registering a NetCDF topo file in setrun.py rundata.topo_data.topofiles.append([4, 'bathy.nc', meta]) That is the only change required in most cases. GeoClaw's Python layer -interrogates the file when you run ``setrun.py`` and writes the -necessary descriptor information into ``topo.data`` automatically. -Variable and coordinate names are auto-detected from CF attributes where -possible; pass ``var_name``, ``lon_name``, or ``lat_name`` explicitly to -``TopoInterrogator`` if auto-detection fails. +interrogates the file when you run ``setrun.py`` and writes the necessary +descriptor information into ``topo.data`` automatically. Variable and coordinate +names are auto-detected from CF attributes where possible; pass ``var_name``, +``lon_name``, or ``lat_name`` explicitly to ``TopoInterrogator`` if +auto-detection fails. This also should handle multiple different coordinate +layouts. Domain subsetting (crop) ^^^^^^^^^^^^^^^^^^^^^^^^ @@ -86,7 +103,7 @@ without creating a smaller file: {'crop_bounds': [-100, -80, 20, 35]}]) Only the subset is read into memory at runtime. The full file is never -loaded. +loaded. Note this is the same as using the Checking CF compliance ^^^^^^^^^^^^^^^^^^^^^^ diff --git a/doc/topo.rst b/doc/topo.rst index 42d4853..22a15b8 100644 --- a/doc/topo.rst +++ b/doc/topo.rst @@ -12,12 +12,13 @@ Topography data - :ref:`topo_order` - :ref:`tsunamidata` - :ref:`dtopo` + - :ref:`g_input` The :ref:`geoclaw` software for flow over topography requires at least one topo file to be input, see :ref:`setrun_geoclaw`. -Currently topo files are restricted to three possible formats as ASCII files, or -NetCDF files are also allowed. +Currently topo files are restricted to three possible formats as ASCII files, +and netCDF files. In the descriptions below it is assumed that the topo file gives the elevation of the topography (relative to some reference level) as a value of @@ -169,14 +170,8 @@ The recognized topotypes are: **topotype = 4** - This file type is not ASCII but rather in a NetCDF4 format supported by the - `CF MetaData conventions (v. 1.6) `_. Files - that conform to this standard can be read in by GeoClaw. The `topotools` - module also has support for reading and writing (including therefore - conversion) of these types of bathymetry files (see :ref:`topo_netcdf` - below). To use this functionality - you will need to add *-DNETCDF* to the *FFLAGS* variable either by the - command line or in the Makefile. + Read in netCDF formatted topography data. See :ref:`netcdf_input` for + NetCDF topography format information. The Fortran code will recognize headers for *topotype* 2 or 3 that have the labels first and then the parameter values. @@ -214,21 +209,24 @@ Several on-line databases are available for topograpy, see Some Python tools for working with topography files are available, see :ref:`topotools`. -.. _topo_netcdf: - -NetCDF format -^^^^^^^^^^^^^ - -Topofiles can be read in netCDF format, either from local `.nc` files or -from some online databases that provide netCDF servers, e.g. the -`NOAA THREDDS server `_. -Use the -`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 online netCDF file (e.g. coastal DEMs at 1/3 arcsecond -resolution are typically several gigabytes). +.. _noaa_thredds: + +NOAA THREDDS server +^^^^^^^^^^^^^^^^^^^^^ + +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 +`_ 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 +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 From 68acd0a45fe72d4b73d9934093aec06e33979ea7 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Thu, 16 Apr 2026 17:08:06 -0400 Subject: [PATCH 03/16] Remove mention of NetCDF output for AMR --- doc/output_styles.rst | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/doc/output_styles.rst b/doc/output_styles.rst index fefb5f8..ac6cae8 100644 --- a/doc/output_styles.rst +++ b/doc/output_styles.rst @@ -134,8 +134,6 @@ Each patch has a unique `grid_number` that usually isn't needed for visualization purposes. - - .. _output_binary: Raw binary output data formats @@ -163,15 +161,6 @@ Unlike the ASCII data files, the binary output files contain ghost cells as well as the interior cells (since a contiguous block of memory is dumped for each patch with a single `write` statement). - -.. _output_netcdf: - -NetCDF output data format ------------------------------- - -NetCDF output is not currently supported in Clawpack. This is not a suitable -format for AMR style data. - .. _output_aux: Output of aux arrays @@ -185,5 +174,3 @@ have the same format as the `q` data, as specifed by `output_format`. Set `output_aux_onlyonce` to `True` if the `aux` arrays do not change with time and you only want to output these arrays once. - - From 893bce5ab30a5a12138399b7c973640dd7c3dd29 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Fri, 17 Apr 2026 19:27:22 -0400 Subject: [PATCH 04/16] Add changes needed for met forcing Some of the less common cases are now handled better with auto-detection of variables and better support for time offsets. --- doc/netcdf.rst | 160 +++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 128 insertions(+), 32 deletions(-) diff --git a/doc/netcdf.rst b/doc/netcdf.rst index 7ef3bb6..3050123 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -171,43 +171,136 @@ In ``setrun.py``: surge_data.storm_specification_type = 'data' surge_data.storm_file = 'isaac.storm' -Then create the storm descriptor file, *e.g.* ``isaac.storm``, using the -Python API: +Then create the storm descriptor file, *e.g.* ``isaac.storm``. GeoClaw +uses a two-stage discovery process: standard variable names are found +automatically from CF ``axis`` / ``standard_name`` attributes and +built-in fallback lists; you only need to supply ``var_mapping`` for +roles whose names are non-standard. Coordinate names (``lon``, ``lat``, +``time``) are always discovered automatically and never need to be +specified. + +**Case 1 — all standard names (ERA5 and similar):** + +ERA5 variable names (``u10``, ``v10``, ``msl``) are in the built-in +fallback lists and are found automatically. No ``var_mapping`` is +required: .. code:: python - from clawpack.geoclaw.surge.storm import Storm import numpy as np + from clawpack.geoclaw.surge.storm import Storm storm = Storm() storm.time_offset = np.datetime64('2012-08-29') storm.file_format = 'netcdf' - storm.file_paths = ['path/to/forcing.nc'] + storm.file_paths = ['path/to/era5_forcing.nc'] storm.write('isaac.storm', file_format='data') -If your variable or dimension names are non-standard, provide a mapping: +**Case 2 — mixed: standard wind names, non-standard pressure:** + +A common pattern when a file has multiple pressure fields (e.g. mean +sea-level and surface pressure) or a non-standard pressure variable +name. Specify only the roles that cannot be discovered automatically; +the rest are still found from the fallback lists: .. code:: python storm.write('isaac.storm', file_format='data', - dim_mapping={'t': 'valid_time'}, - var_mapping={'wind_u': 'u10', 'wind_v': 'v10', - 'pressure': 'msl'}) + var_mapping={'pressure': 'prmsl'}) -Time handling -^^^^^^^^^^^^^ +Any name supplied in ``var_mapping`` is validated against the variables +actually present in the file before the descriptor is written, so a +typo raises an informative error immediately rather than producing +incorrect output silently. + +**Case 3 — all non-standard names (e.g. NWS13/OWI-NetCDF):** + +Supply all three roles explicitly when none of the variable names match +the built-in fallback lists: + +.. code:: python + + storm = Storm() + storm.time_offset = np.datetime64('2012-08-29') + storm.file_format = 'nws13' + storm.file_paths = ['path/to/nws13_forcing.nc'] + storm.write('isaac.storm', file_format='data', + var_mapping={'wind_u': 'uwnd', + 'wind_v': 'vwnd', + 'pressure': 'press'}) -GeoClaw works in seconds from a user-defined offset (typically landfall -or storm genesis). All CF time decoding -- including calendar handling -and unit conversion from hours/days -- is done in Python. The Fortran -runtime sees only seconds from offset. +**Advanced: pre-built MetInterrogator:** -Set the offset when constructing the storm: +If you need direct control over CF validation, unit checking, or +lon/lat convention detection before the descriptor is written, you can +construct a :class:`~clawpack.geoclaw.netcdf_utils.MetInterrogator` +explicitly and pass it via ``met_interrogator``. When a pre-built +interrogator is supplied, auto-discovery and ``var_mapping`` validation +are bypassed entirely: + +.. code:: python + + from clawpack.geoclaw.netcdf_utils import MetInterrogator + + mi = MetInterrogator('path/to/forcing.nc', + variable_map={'wind_u': 'u10', + 'wind_v': 'v10', + 'pressure': 'msl'}) + storm.write('isaac.storm', file_format='data', met_interrogator=mi) + +.. note:: + + ``dim_mapping`` is accepted by ``write()`` / ``write_data()`` for + backwards compatibility but has no effect -- coordinate names are + always discovered automatically via CF conventions. + +.. note:: + + ``storm.window`` (ramp width and application domain) and + ``MetInterrogator``'s ``crop_bounds`` (read-time spatial subset of + the NetCDF file) are independent. Setting one does not affect the + other. + + +Time handling +^^^^^^^^^^^^^^^ + +GeoClaw works in seconds relative to a user-defined reference time +(typically landfall or storm genesis). Set it when constructing the +storm: .. code:: python storm.time_offset = np.datetime64('2012-08-29T00:00') # landfall time +Python computes ``nc_time_offset`` — the elapsed seconds from +``storm.time_offset`` to the first record in the file — and writes it +to the descriptor. This value is independent of how the file encodes +time internally (Unix epoch, local epoch, hours-since-reference, etc.). + +Fortran converts raw time values using: + +.. code:: + + storm_time[i] = raw[i] − raw[0] + nc_time_offset + +Subtracting ``raw[0]`` converts any absolute encoding to +elapsed-since-first-record; adding ``nc_time_offset`` then anchors that +to ``storm.time_offset``. No unit conversion (hours→seconds, +days→seconds) is performed in Fortran: the descriptor stores +``nc_time_offset`` in seconds and the raw values are read as stored. +``write_data`` passes ``time_reference=storm.time_offset`` to +``MetInterrogator`` automatically; no additional user action is needed. + +Longitude convention +^^^^^^^^^^^^^^^^^^^^^^ + +Longitude convention ([0, 360] vs [-180, 180]) is detected and +normalized automatically by ``MetInterrogator``. ERA5 files that store +longitude in [0, 360] do not need to be preprocessed to [-180, 180] +before use. The detected convention is written to the descriptor and +Fortran applies index normalization at runtime without copying arrays. + -------------- Developer Reference @@ -219,14 +312,16 @@ Architecture overview The system has a strict Python/Fortran split: **Python** handles: file interrogation, CF attribute parsing, coordinate -convention detection, fill value resolution, unit conversion, time -decoding, crop bound validation, and descriptor writing. +convention detection, fill value resolution, unit conversion, +``nc_time_offset`` computation (elapsed seconds from ``storm.time_offset`` +to the first record), crop bound validation, and descriptor writing. **Fortran** handles: opening the NetCDF file at runtime using information from the descriptor, index arithmetic for coordinate normalization (no data copies), domain subsetting via -``start``/``count`` arguments to ``nf90_get_var``, and time-slice reads -for met forcing. +``start``/``count`` arguments to ``nf90_get_var``, time-slice reads +for met forcing, and converting raw time values to seconds from +``storm.time_offset`` via ``raw[i] − raw[0] + nc_time_offset``. Fortran assumes the unit contract from ``GEOCLAW_NETCDF_UNITS`` in ``units.py`` without checking. Python enforces it. @@ -351,21 +446,21 @@ silent Fortran ``STOP`` produces no useful output. Known technical debt ~~~~~~~~~~~~~~~~~~~~ -- ``util.get_netcdf_names`` and the coordinate/variable discovery logic - in ``NetCDFInterrogator`` are parallel implementations. They should - be consolidated -- ``NetCDFInterrogator`` should be the single source - and ``util.get_netcdf_names`` should delegate to it. This is deferred - to the surge module refactor. -- ``Storm.write(file_format="data")`` does not yet call - ``DescriptorWriter`` directly for NetCDF storm entries. Currently the - tests write descriptors manually. The integration should happen as - part of the surge refactor. +- ``util.get_netcdf_names`` (in ``util.py``) remains a parallel + implementation of variable name discovery alongside + ``NetCDFInterrogator``. ``Storm.write_data`` now uses + ``MetInterrogator`` for NetCDF met forcing (format 2), but falls back + to ``util.get_netcdf_names`` for variable name discovery when no + explicit ``var_mapping`` is provided. Full consolidation -- making + ``util.get_netcdf_names`` delegate to ``NetCDFInterrogator`` and + removing the duplicate -- remains a TODO. - WRF raw output requires ``MetPreprocessor`` for string-time decoding (``Times`` character array, not a numeric CF time variable) and curvilinear grid handling (``XLAT``/``XLONG`` are 2D time-varying arrays). A skip-marked test stub documents the gap in ``test_storm.py``. + Adding a new met forcing field (e.g. precipitation) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -406,10 +501,11 @@ binary files are committed to the repository. Next steps ---------- -- **Surge module refactor**: consolidate ``util.get_netcdf_names`` into - ``NetCDFInterrogator``, wire ``Storm.write`` to use - ``DescriptorWriter``, clean up parallel discovery paths in - ``storm.py`` +- **Consolidate util.get_netcdf_names**: ``Storm.write_data`` now + calls ``MetInterrogator`` and ``DescriptorWriter`` for NetCDF met + forcing; the remaining step is to make ``util.get_netcdf_names`` + delegate to ``NetCDFInterrogator`` and remove the duplicate discovery + logic in ``storm.py`` - **WRF support**: implement ``MetPreprocessor`` to handle string-time axis and curvilinear grid; the skip-marked test stub in ``test_storm.py`` documents the expected interface From 8bb4891133820799acd9195e1e0040f651b19a8f Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Thu, 7 May 2026 10:44:44 -0400 Subject: [PATCH 05/16] Update netcdf docs for topography var name change --- doc/netcdf.rst | 40 ++++++++++++++++++++++++++++------------ 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/doc/netcdf.rst b/doc/netcdf.rst index 3050123..3e989a3 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -61,34 +61,50 @@ What your file must provide Registering a NetCDF topo file in setrun.py ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -If your file meets the above requirements and either has the variable `z` or -only one 2D variable in the file, you can simply use the following in +If your file meets the above requirements, you can simply use the following in ``setrun.py``: .. code:: python - + rundata.topo_data.topofiles.append([4, 'bathy.nc']) This matches with what is expected for backwards compatibility. -If your file meets the above requirements and has a non-standard variable name -or you want to specify crop bounds you can use the Python API to interrogate the -file and write a descriptor: +If you want to specify crop bounds, or if GeoClaw cannot find the elevation +variable automatically, use the Python API to interrogate the file and write a +descriptor: .. code:: python from clawpack.geoclaw.netcdf_utils import TopoInterrogator - meta = TopoInterrogator('bathy.nc', var_name='z', + meta = TopoInterrogator('bathy.nc', crop_bounds=(-100, -60, 15, 35)).interrogate_topo() rundata.topo_data.topofiles.append([4, 'bathy.nc', meta]) That is the only change required in most cases. GeoClaw's Python layer interrogates the file when you run ``setrun.py`` and writes the necessary -descriptor information into ``topo.data`` automatically. Variable and coordinate -names are auto-detected from CF attributes where possible; pass ``var_name``, -``lon_name``, or ``lat_name`` explicitly to ``TopoInterrogator`` if -auto-detection fails. This also should handle multiple different coordinate -layouts. +descriptor information into ``topo.data`` automatically. + +The elevation variable is found automatically using a two-step search: + +1. **CF** ``standard_name`` **attribute** — variables whose ``standard_name`` + is one of ``surface_altitude``, ``height_above_mean_sea_level``, + ``height_above_reference_ellipsoid``, ``bedrock_altitude``, ``altitude``, + ``height``, or ``sea_floor_depth_below_geoid`` are matched first. +2. **Common variable names** — if no CF match is found, the variable name is + checked against a built-in list that includes ``z``, ``elevation``, + ``topo``, ``height``, ``altitude``, ``depth``, ``dem``, ``bathymetry``, + ``bathy``, ``Band1``, and several capitalisation variants. + +Pass ``var_name`` explicitly only when none of the names above match your file, +or when the file contains multiple variables that would both match and you need +to disambiguate:: + + meta = TopoInterrogator('bathy.nc', var_name='my_elevation', + crop_bounds=(-100, -60, 15, 35)).interrogate_topo() + +Coordinate names (``lon``/``latitude``/``x`` etc.) and dimension ordering are +always discovered automatically and never need to be specified. Domain subsetting (crop) ^^^^^^^^^^^^^^^^^^^^^^^^ From 917e494a274544f1250a45a068f287fed71936d9 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Thu, 7 May 2026 23:15:49 -0400 Subject: [PATCH 06/16] Add .gitignore Mostly just ignores build and doctree --- .gitignore | 65 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..dd68408 --- /dev/null +++ b/.gitignore @@ -0,0 +1,65 @@ +# ----------------------------------------------------------------------------- +# Sphinx build output +# ----------------------------------------------------------------------------- +_build1 + +# Generated API docs (if autogenerated) +api/generated/ + +# ----------------------------------------------------------------------------- +# Python +# ----------------------------------------------------------------------------- +__pycache__/ +*.py[cod] +*$py.class + +.pytest_cache/ +.mypy_cache/ +.ruff_cache/ +.coverage +coverage.xml +htmlcov/ + +# ----------------------------------------------------------------------------- +# Virtual environments +# ----------------------------------------------------------------------------- +.venv/ +venv/ +env/ + +# ----------------------------------------------------------------------------- +# Jupyter +# ----------------------------------------------------------------------------- +.ipynb_checkpoints/ + +# ----------------------------------------------------------------------------- +# Documentation tooling caches +# ----------------------------------------------------------------------------- +.doctrees/ + +# If using sphinx-gallery +auto_examples/ + +# If using jupyter-cache / myst-nb +.jupyter_cache/ + +# ----------------------------------------------------------------------------- +# OS / editor files +# ----------------------------------------------------------------------------- +.DS_Store +Thumbs.db + +.vscode/ +.idea/ + +# Vim +*.swp +*.swo + +# ----------------------------------------------------------------------------- +# Package managers / environments +# ----------------------------------------------------------------------------- +# uv +.uv/ +.pixi/ +conda-meta/ From 4788a81a1869210b67d51fba0ee7e3feefbeec9a Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Thu, 7 May 2026 23:17:56 -0400 Subject: [PATCH 07/16] Add docs for lon wrapping support in netcdf --- doc/netcdf.rst | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/doc/netcdf.rst b/doc/netcdf.rst index 3e989a3..bcea652 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -71,18 +71,19 @@ If your file meets the above requirements, you can simply use the following in This matches with what is expected for backwards compatibility. If you want to specify crop bounds, or if GeoClaw cannot find the elevation -variable automatically, use the Python API to interrogate the file and write a -descriptor: +variable automatically, use ``topo_entries()`` to interrogate the file: .. code:: python from clawpack.geoclaw.netcdf_utils import TopoInterrogator - meta = TopoInterrogator('bathy.nc', - crop_bounds=(-100, -60, 15, 35)).interrogate_topo() - rundata.topo_data.topofiles.append([4, 'bathy.nc', meta]) - -That is the only change required in most cases. GeoClaw's Python layer -interrogates the file when you run ``setrun.py`` and writes the necessary + with TopoInterrogator('bathy.nc', crop_bounds=(-100, -60, 15, 35)) as intr: + rundata.topo_data.topofiles.extend(intr.topo_entries()) + +``topo_entries()`` returns a list of ``[4, path, TopoMetadata]`` entries ready +to pass directly to ``topofiles``. In the common case the list has one entry, +but it may contain two when the crop region straddles the file's longitude cut +point (e.g. a global file cropped across the dateline) — ``extend`` handles +both cases without extra logic. GeoClaw's Python layer writes the necessary descriptor information into ``topo.data`` automatically. The elevation variable is found automatically using a two-step search: @@ -110,16 +111,18 @@ Domain subsetting (crop) ^^^^^^^^^^^^^^^^^^^^^^^^ If your NetCDF file covers a larger area than your simulation domain -(common with global or regional datasets), you can crop at read time -without creating a smaller file: +(common with global or regional datasets), pass ``crop_bounds`` to +``TopoInterrogator`` and use ``topo_entries()`` to register the file: .. code:: python - topo_data.topofiles.append([4, 'gebco_global.nc', - {'crop_bounds': [-100, -80, 20, 35]}]) + from clawpack.geoclaw.netcdf_utils import TopoInterrogator + with TopoInterrogator('gebco_global.nc', + crop_bounds=(-100, -80, 20, 35)) as intr: + rundata.topo_data.topofiles.extend(intr.topo_entries()) Only the subset is read into memory at runtime. The full file is never -loaded. Note this is the same as using the +loaded. Checking CF compliance ^^^^^^^^^^^^^^^^^^^^^^ @@ -388,7 +391,7 @@ Topo (lines in topo.data after topo_type) var_name = z lon_name = longitude lat_name = latitude - lon_convention = 180 + lon_offset = 0.0 lat_order = S_to_N dim_order = lat,lon fill_value = -9999.0 From 3c2a1470a94502aa4c093e79b6d13cb5c5abf7ee Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Fri, 22 May 2026 14:22:11 -0400 Subject: [PATCH 08/16] =?UTF-8?q?Rename=20Interrogator=E2=86=92Inspector?= =?UTF-8?q?=20in=20NetCDF=20docs;=20add=20API=20module=20page=20and=20stor?= =?UTF-8?q?m=5Ftime=5Fscale?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - netcdf.rst: rename all TopoInterrogator/MetInterrogator/NetCDFInterrogator references to TopoInspector/MetInspector/NetCDFInspector, rename met_interrogator parameter to met_inspector, rename test file stubs to match (test_base_inspector.py etc.) - netcdf_utils_module.rst: new Sphinx automodule page covering NetCDFInspector, TopoInspector, MetInspector, CFNormalizer, DescriptorWriter - geoclaw.rst: add netcdf and netcdf_utils_module to the toctree so both the user guide and API reference are reachable from the GeoClaw index - setrun_geoclaw.rst: document the new storm_time_scale parameter Co-Authored-By: Claude Sonnet 4.6 --- doc/geoclaw.rst | 2 + doc/netcdf.rst | 84 ++++++++++++++++++------------------- doc/netcdf_utils_module.rst | 39 +++++++++++++++++ doc/setrun_geoclaw.rst | 11 ++++- 4 files changed, 93 insertions(+), 43 deletions(-) create mode 100644 doc/netcdf_utils_module.rst diff --git a/doc/geoclaw.rst b/doc/geoclaw.rst index c3ee6d0..4dabffe 100644 --- a/doc/geoclaw.rst +++ b/doc/geoclaw.rst @@ -69,6 +69,8 @@ More will eventually appear in the :ref:`apps`. nearshore_interp tsunamidata surgedata + netcdf + netcdf_utils_module marching_front force_dry sphere_source diff --git a/doc/netcdf.rst b/doc/netcdf.rst index bcea652..173fd1c 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -71,13 +71,13 @@ If your file meets the above requirements, you can simply use the following in This matches with what is expected for backwards compatibility. If you want to specify crop bounds, or if GeoClaw cannot find the elevation -variable automatically, use ``topo_entries()`` to interrogate the file: +variable automatically, use ``topo_entries()`` to inspect the file: .. code:: python - from clawpack.geoclaw.netcdf_utils import TopoInterrogator - with TopoInterrogator('bathy.nc', crop_bounds=(-100, -60, 15, 35)) as intr: - rundata.topo_data.topofiles.extend(intr.topo_entries()) + from clawpack.geoclaw.netcdf_utils import TopoInspector + with TopoInspector('bathy.nc', crop_bounds=(-100, -60, 15, 35)) as insp: + rundata.topo_data.topofiles.extend(insp.topo_entries()) ``topo_entries()`` returns a list of ``[4, path, TopoMetadata]`` entries ready to pass directly to ``topofiles``. In the common case the list has one entry, @@ -101,8 +101,8 @@ Pass ``var_name`` explicitly only when none of the names above match your file, or when the file contains multiple variables that would both match and you need to disambiguate:: - meta = TopoInterrogator('bathy.nc', var_name='my_elevation', - crop_bounds=(-100, -60, 15, 35)).interrogate_topo() + meta = TopoInspector('bathy.nc', var_name='my_elevation', + crop_bounds=(-100, -60, 15, 35)).inspect_topo() Coordinate names (``lon``/``latitude``/``x`` etc.) and dimension ordering are always discovered automatically and never need to be specified. @@ -112,14 +112,14 @@ Domain subsetting (crop) If your NetCDF file covers a larger area than your simulation domain (common with global or regional datasets), pass ``crop_bounds`` to -``TopoInterrogator`` and use ``topo_entries()`` to register the file: +``TopoInspector`` and use ``topo_entries()`` to register the file: .. code:: python - from clawpack.geoclaw.netcdf_utils import TopoInterrogator - with TopoInterrogator('gebco_global.nc', - crop_bounds=(-100, -80, 20, 35)) as intr: - rundata.topo_data.topofiles.extend(intr.topo_entries()) + from clawpack.geoclaw.netcdf_utils import TopoInspector + with TopoInspector('gebco_global.nc', + crop_bounds=(-100, -80, 20, 35)) as insp: + rundata.topo_data.topofiles.extend(insp.topo_entries()) Only the subset is read into memory at runtime. The full file is never loaded. @@ -178,7 +178,7 @@ Time seconds from user-defined offset ============================= ================================ If your file uses hPa/mbar for pressure or knots for wind, -``MetInterrogator`` will convert automatically. +``MetInspector`` will convert automatically. Registering a NetCDF storm file ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -248,24 +248,24 @@ the built-in fallback lists: 'wind_v': 'vwnd', 'pressure': 'press'}) -**Advanced: pre-built MetInterrogator:** +**Advanced: pre-built MetInspector:** If you need direct control over CF validation, unit checking, or lon/lat convention detection before the descriptor is written, you can -construct a :class:`~clawpack.geoclaw.netcdf_utils.MetInterrogator` -explicitly and pass it via ``met_interrogator``. When a pre-built -interrogator is supplied, auto-discovery and ``var_mapping`` validation +construct a :class:`~clawpack.geoclaw.netcdf_utils.MetInspector` +explicitly and pass it via ``met_inspector``. When a pre-built +inspector is supplied, auto-discovery and ``var_mapping`` validation are bypassed entirely: .. code:: python - from clawpack.geoclaw.netcdf_utils import MetInterrogator + from clawpack.geoclaw.netcdf_utils import MetInspector - mi = MetInterrogator('path/to/forcing.nc', - variable_map={'wind_u': 'u10', - 'wind_v': 'v10', - 'pressure': 'msl'}) - storm.write('isaac.storm', file_format='data', met_interrogator=mi) + mi = MetInspector('path/to/forcing.nc', + variable_map={'wind_u': 'u10', + 'wind_v': 'v10', + 'pressure': 'msl'}) + storm.write('isaac.storm', file_format='data', met_inspector=mi) .. note:: @@ -276,7 +276,7 @@ are bypassed entirely: .. note:: ``storm.window`` (ramp width and application domain) and - ``MetInterrogator``'s ``crop_bounds`` (read-time spatial subset of + ``MetInspector``'s ``crop_bounds`` (read-time spatial subset of the NetCDF file) are independent. Setting one does not affect the other. @@ -309,13 +309,13 @@ to ``storm.time_offset``. No unit conversion (hours→seconds, days→seconds) is performed in Fortran: the descriptor stores ``nc_time_offset`` in seconds and the raw values are read as stored. ``write_data`` passes ``time_reference=storm.time_offset`` to -``MetInterrogator`` automatically; no additional user action is needed. +``MetInspector`` automatically; no additional user action is needed. Longitude convention ^^^^^^^^^^^^^^^^^^^^^^ Longitude convention ([0, 360] vs [-180, 180]) is detected and -normalized automatically by ``MetInterrogator``. ERA5 files that store +normalized automatically by ``MetInspector``. ERA5 files that store longitude in [0, 360] do not need to be preprocessed to [-180, 180] before use. The detected convention is written to the descriptor and Fortran applies index normalization at runtime without copying arrays. @@ -330,7 +330,7 @@ Architecture overview The system has a strict Python/Fortran split: -**Python** handles: file interrogation, CF attribute parsing, coordinate +**Python** handles: file inspection, CF attribute parsing, coordinate convention detection, fill value resolution, unit conversion, ``nc_time_offset`` computation (elapsed seconds from ``storm.time_offset`` to the first record), crop bound validation, and descriptor writing. @@ -350,7 +350,7 @@ Class hierarchy (``netcdf_utils.py``) :: - NetCDFInterrogator + NetCDFInspector - open file (xarray, Dask-lazy) - discover coordinate variables by name heuristics + CF standard_name - detect lon convention, lat order, dim order @@ -358,12 +358,12 @@ Class hierarchy (``netcdf_utils.py``) - validate crop bounds against file extent - output: NetCDFDescriptor dataclass - TopoInterrogator(NetCDFInterrogator) + TopoInspector(NetCDFInspector) - detect fill values within crop region (hard error) - verify and convert units to contract - no multi-file coverage logic (Fortran handles compositing) - MetInterrogator(NetCDFInterrogator) + MetInspector(NetCDFInspector) - check wind_u, wind_v, pressure present and on same grid/time axis - convert units to contract - decode CF time to seconds from offset @@ -467,11 +467,11 @@ Known technical debt - ``util.get_netcdf_names`` (in ``util.py``) remains a parallel implementation of variable name discovery alongside - ``NetCDFInterrogator``. ``Storm.write_data`` now uses - ``MetInterrogator`` for NetCDF met forcing (format 2), but falls back + ``NetCDFInspector``. ``Storm.write_data`` now uses + ``MetInspector`` for NetCDF met forcing (format 2), but falls back to ``util.get_netcdf_names`` for variable name discovery when no explicit ``var_mapping`` is provided. Full consolidation -- making - ``util.get_netcdf_names`` delegate to ``NetCDFInterrogator`` and + ``util.get_netcdf_names`` delegate to ``NetCDFInspector`` and removing the duplicate -- remains a TODO. - WRF raw output requires ``MetPreprocessor`` for string-time decoding (``Times`` character array, not a numeric CF time variable) and @@ -484,13 +484,13 @@ Adding a new met forcing field (e.g. precipitation) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1. Add the role and unit to ``GEOCLAW_NETCDF_UNITS`` in ``units.py`` -2. Add detection logic to ``MetInterrogator`` for the new variable +2. Add detection logic to ``MetInspector`` for the new variable 3. ``DescriptorWriter`` requires no change -- new ``&variable_info`` - blocks are written automatically for any role the interrogator + blocks are written automatically for any role the inspector returns 4. Add the Fortran reader logic in the storm NetCDF module to consume the new ``geoclaw_role`` value -5. Add unit tests in ``tests/netcdf/test_met_interrogator.py`` +5. Add unit tests in ``tests/netcdf/test_met_inspector.py`` 6. Add a regression test in the appropriate storm surge example Test structure @@ -500,9 +500,9 @@ Test structure tests/netcdf/ unit tests for netcdf_utils.py conftest.py in-memory NetCDF fixtures - test_base_interrogator.py - test_topo_interrogator.py - test_met_interrogator.py + test_base_inspector.py + test_topo_inspector.py + test_met_inspector.py test_descriptor_writer.py test_cf_normalizer.py @@ -521,19 +521,19 @@ Next steps ---------- - **Consolidate util.get_netcdf_names**: ``Storm.write_data`` now - calls ``MetInterrogator`` and ``DescriptorWriter`` for NetCDF met + calls ``MetInspector`` and ``DescriptorWriter`` for NetCDF met forcing; the remaining step is to make ``util.get_netcdf_names`` - delegate to ``NetCDFInterrogator`` and remove the duplicate discovery + delegate to ``NetCDFInspector`` and remove the duplicate discovery logic in ``storm.py`` - **WRF support**: implement ``MetPreprocessor`` to handle string-time axis and curvilinear grid; the skip-marked test stub in ``test_storm.py`` documents the expected interface - **Precipitation field**: reserved in ``GEOCLAW_NETCDF_UNITS``, add - ``MetInterrogator`` detection and Fortran consumer + ``MetInspector`` detection and Fortran consumer - **Friction field**: same pattern as precipitation - **CLI tool**: expose ``CFNormalizer`` as a command-line utility for users who want to check or repair their NetCDF files before running GeoClaw - **Duplicate code audit**: other places in the GeoClaw codebase that do ad-hoc NetCDF variable discovery should be identified and migrated - to ``NetCDFInterrogator`` + to ``NetCDFInspector`` diff --git a/doc/netcdf_utils_module.rst b/doc/netcdf_utils_module.rst new file mode 100644 index 0000000..a760778 --- /dev/null +++ b/doc/netcdf_utils_module.rst @@ -0,0 +1,39 @@ + +.. _netcdf_utils_module: + +netcdf_utils module for NetCDF input +===================================== + +.. seealso:: + - :ref:`netcdf_input` + - :ref:`topotools_module` + - :ref:`storm_module` + +The ``netcdf_utils`` module provides the Python layer for reading NetCDF +topography and meteorological forcing files in GeoClaw. It handles CF +attribute parsing, coordinate convention detection, unit conversion, and +descriptor writing so that Fortran only needs to open the file and read +pre-validated indices. + +The main classes are: + +- :class:`~clawpack.geoclaw.netcdf_utils.NetCDFInspector` — base class; + coordinate discovery, fill-value resolution, crop-bound validation. +- :class:`~clawpack.geoclaw.netcdf_utils.TopoInspector` — topography + subclass; fill-value checking within the crop region, unit enforcement. +- :class:`~clawpack.geoclaw.netcdf_utils.MetInspector` — meteorological + subclass; wind/pressure variable discovery, CF time decoding, unit + conversion. +- :class:`~clawpack.geoclaw.netcdf_utils.CFNormalizer` — adds or repairs + CF attributes in place without modifying data values. +- :class:`~clawpack.geoclaw.netcdf_utils.DescriptorWriter` — writes the + key=value topo descriptor lines or the ``&file_info`` / ``&variable_info`` + Fortran namelist blocks for met forcing. + +Documentation auto-generated from the module docstrings +-------------------------------------------------------- + +.. automodule:: clawpack.geoclaw.netcdf_utils + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/setrun_geoclaw.rst b/doc/setrun_geoclaw.rst index e67bc2b..0f94d95 100644 --- a/doc/setrun_geoclaw.rst +++ b/doc/setrun_geoclaw.rst @@ -453,7 +453,16 @@ Storm Specification Data .. attribute:: rundata.surge_data.storm_file : string Specifies the path to the storm data. IF `storm_specification_type > 0` then - this should point to a GeoClaw formatted storm file (see :ref:`storm_module` for + this should point to a GeoClaw formatted storm file (see :ref:`storm_module` for details). If `storm_specification < 0` then this should either specify a path to files specifying the storm or a single file depending on the type of input data. + +.. attribute:: rundata.surge_data.storm_time_scale : float + + A multiplicative scale factor applied to the storm's time axis. Values + greater than 1 slow the storm down (it takes longer to traverse the + domain); values less than 1 speed it up. The default is ``1.0`` (no + scaling). This is primarily useful for sensitivity studies and for + creating synthetic storms whose tracks are spatially identical to a + historical event but travel at a different speed. From 4441249b5c69fde2029232e1c65522d665764635 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Sun, 31 May 2026 20:20:33 -0600 Subject: [PATCH 09/16] Add documentation related to topography data This is an update that reflects adding additional parameters to the Topography class in topotools.py. This now includes additional parameters that impact both netCDF and ASCII topography files. --- doc/changes_to_master.rst | 41 ++++++ doc/geoclaw.rst | 1 + doc/setrun_geoclaw.rst | 102 +++++++++++--- doc/topo_order.rst | 25 ++-- doc/topodata_format.rst | 283 ++++++++++++++++++++++++++++++++++++++ doc/topotools.rst | 74 ++++++++++ 6 files changed, 498 insertions(+), 28 deletions(-) create mode 100644 doc/topodata_format.rst diff --git a/doc/changes_to_master.rst b/doc/changes_to_master.rst index 7443453..d174d42 100644 --- a/doc/changes_to_master.rst +++ b/doc/changes_to_master.rst @@ -64,6 +64,47 @@ See `amrclaw diffs Changes to geoclaw ------------------ +- **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. The first file listed in ``topo.data`` is assigned rank 1 + (highest priority) by Fortran; no Fortran-side sorting occurs. + ``rundata.topo_data.override_order = True`` preserves user-specified list + order; when used, the finest (highest-resolution) file should be listed + first. 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:: + + t = Topography() + t.read('old.tt1', topo_type=1) # DeprecationWarning + t.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. See `geoclaw diffs `_ diff --git a/doc/geoclaw.rst b/doc/geoclaw.rst index 4dabffe..a206288 100644 --- a/doc/geoclaw.rst +++ b/doc/geoclaw.rst @@ -49,6 +49,7 @@ More will eventually appear in the :ref:`apps`. topo grid_registration topo_order + topodata_format topotools dtopo kmltools_module diff --git a/doc/setrun_geoclaw.rst b/doc/setrun_geoclaw.rst index 0f94d95..d1a9029 100644 --- a/doc/setrun_geoclaw.rst +++ b/doc/setrun_geoclaw.rst @@ -155,26 +155,88 @@ See :ref:`topo` for more information about specifying topography (and bathymetry) data files in GeoClaw. -.. attribute:: rundata.topo_data.topofiles : list of lists - - *topofiles* should be a list of the form *[file1info, file2info, etc.]* - where each element is itself a list of the form - - [topotype, fname] - - with values - - *topotype* : integer - - 1,2 or 3 depending on the format of the file (see :ref:`topo`). - - *fname* : string - - the name of the topo file. - - **Note:** Starting in v5.8.0 implicitly specifying a flag region for - AMR is no longer supported in the specification of a topo file. - For more about controlling AMR in various regions, see :ref:`flagregions`. +.. attribute:: rundata.topo_data.topofiles : list + + A list of topography file entries. The preferred form uses + :class:`~clawpack.geoclaw.topotools.Topography` objects:: + + from clawpack.geoclaw.topotools import Topography + + t = Topography() + t.path = 'etopo1.tt2' + t.topo_type = 2 + rundata.topo_data.topofiles.append(t) + + .. deprecated:: + The ``[topotype, fname]`` list format is deprecated. It still works + but emits a ``DeprecationWarning``. Use ``Topography`` objects as + shown above. + + **Note:** Starting in v5.8.0 implicitly specifying a flag region for + AMR is no longer supported in the specification of a topo file. + For more about controlling AMR in various regions, see :ref:`flagregions`. + + See :ref:`topo` for information about the topo file formats + (``topo_type`` values 2, 3, 4). + + See :ref:`topo_order` for how priority is determined when topo files + overlap, and :ref:`topodata_format` for the ``topo.data`` file format. + +.. _setrun_topo_preprocessing: + +Topography preprocessing attributes +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Each :class:`~clawpack.geoclaw.topotools.Topography` entry in ``topofiles`` +supports seven preprocessing attributes that are applied automatically when +the file is loaded by :meth:`~clawpack.geoclaw.topotools.Topography.read`. +Set them on the object before appending to ``topofiles``: + ++----------------+----------------+-----------+-----------------------------------------------+ +| Attribute | Type | Default | Description | ++================+================+===========+===============================================+ +| ``crop_extent``| list or None | ``None`` | ``[x1, x2, y1, y2]`` crop region. | ++----------------+----------------+-----------+-----------------------------------------------+ +| ``coarsen`` | int | ``1`` | Stride-subsampling factor (1 = no coarsen). | ++----------------+----------------+-----------+-----------------------------------------------+ +| ``buffer`` | int | ``0`` | Grid-point margin outside crop region. | ++----------------+----------------+-----------+-----------------------------------------------+ +| ``align`` | tuple or None | ``None`` | ``(x, y)`` alignment for coarsened grids. | ++----------------+----------------+-----------+-----------------------------------------------+ +| ``x_shift`` | float | ``0.0`` | Constant added to all x coordinates. | ++----------------+----------------+-----------+-----------------------------------------------+ +| ``z_shift`` | float | ``0.0`` | Constant added to non-missing Z values. | ++----------------+----------------+-----------+-----------------------------------------------+ +| ``negate_z`` | bool | ``False`` | If ``True``, flip the sign of all Z values. | ++----------------+----------------+-----------+-----------------------------------------------+ + +.. note:: + + **Non-obvious behaviors:** + + - ``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. + + - ``buffer`` is in grid points, not degrees. Float values are truncated + via ``int()``, so ``buffer=0.5`` is equivalent to ``buffer=0``. + + - ``negate_z`` is independent of the ``topo_type < 0`` sign convention. + 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. + +Example using several attributes:: + + t = Topography() + t.path = 'gebco_2023.nc' + t.topo_type = 4 + t.crop_extent = [-100., -60., 10., 50.] + t.coarsen = 2 + t.z_shift = 10.0 + rundata.topo_data.topofiles.append(t) .. attribute:: rundata.dtopo_data.dtopofiles : list of lists diff --git a/doc/topo_order.rst b/doc/topo_order.rst index 6fe2020..cbc740b 100644 --- a/doc/topo_order.rst +++ b/doc/topo_order.rst @@ -21,6 +21,11 @@ regular spacing `dx` and `dy`, the resolutions are compared by computing the cell areas `dx*dy`. If two or more topofiles overlap a point in the cell then the finer one is used. +This ordering is computed entirely in Python by +``TopographyData._compute_priority_order`` and written into ``topo.data`` +before the Fortran run. The first file listed in ``topo.data`` is assigned +rank 1 (highest priority) by Fortran. No Fortran-side sorting occurs. + This usually works fine but may not always be what is wanted. In particular there may be two DEMs at essentially the same resolution, but one is a more recent and hence a better dataset to use @@ -31,8 +36,8 @@ resolutions are essentially identical, `dx*dy` maybe be slightly different (perhaps at the roundoff level) and so it could be that the older one is viewed as "finer" and would be used preferentially. Even if `dx*dy` is exactly the same, it is not obvious which one would be used in GeoClaw. -(The one listed later in `rundata.topo_data.topofiles` -would be considered finer in this case, with a relative tolerance of +(The one listed first in `rundata.topo_data.topofiles` +is given higher priority in this case, with a relative tolerance of 0.001 used to judge whether the resolutions are the same). .. warning :: This relative tolerance was introduced in v5.14.0, @@ -48,12 +53,16 @@ If `setrun.py` specifies:: rundata.topo_data.override_order = True -then the preference order is no longer based on the computed -resolution. Instead it is assumed that the topo files in -`rundata.topo_data.topofiles` are ordered from least preferable to most -preferable (so normally you would list the coarsest DEM first, going down to -the finest DEM, in order to reproduce the default behavior). -This allows precise specification of the order of preference. +then the preference order is no longer based on the computed resolution. +Instead the list order in `rundata.topo_data.topofiles` is used as-is: the +first file in the list is given the highest priority (rank 1). Normally you +should list the finest (highest-resolution) DEM first and the coarsest last, +to reproduce the default behavior. This allows precise specification of the +order of preference when automatic sorting by resolution is not appropriate. + +``override_order`` is a Python-only attribute of ``rundata.topo_data``. +It controls how Python orders the list before writing ``topo.data`` but is +not itself written to ``topo.data`` and has no Fortran-side counterpart. Resolution of dtopo files ------------------------- diff --git a/doc/topodata_format.rst b/doc/topodata_format.rst new file mode 100644 index 0000000..08ac60c --- /dev/null +++ b/doc/topodata_format.rst @@ -0,0 +1,283 @@ +.. _topodata_format: + +topo.data File Format +===================== + +.. note:: + ``topo.data`` is generated automatically by ``make data`` or by calling + :meth:`TopographyData.write() `. + It is not intended to be hand-authored. The format is tied to the GeoClaw + version; files generated with an older build may be incompatible with a + newer Fortran binary. + +Overview +-------- + +``topo.data`` is read by the Fortran routine ``read_topo_settings`` in +``topo_module.f90``. It describes every topography file used in a simulation: +how to locate the file, its format type, and preprocessing parameters to apply +at load time. + +The file has two sections: + +1. A 3-line **global header**. +2. One **per-file block** for each topography file. + + +Global Header (3 lines) +----------------------- + +:: + + # replace no_data_value in topofile + # (Type topography specification) + # number of topo files + +``topo_missing`` + Float sentinel (default ``99999.0``). Fortran replaces each file's own + ``no_data_value`` with this value after loading. + +``test_topography`` + Integer. ``0`` = use the file list that follows. Values ``1``–``3`` + select built-in analytic test bathymetries (jump discontinuity, oceanic + shelf); in those cases no per-file blocks are written. + +``ntopofiles`` + Integer count of per-file blocks that follow. + +.. note:: + ``override_order`` is a Python-only attribute of + :class:`~clawpack.geoclaw.data.TopographyData`; it controls how Python + sorts the files before writing but is **not** written to ``topo.data``. + + +Per-File Block (9 lines, all types) +------------------------------------ + +One block follows for each topography file, in priority order (finest +resolution first; see :ref:`priority_convention` below). + +:: + + '' # topo_path + # topo_type + # crop_extent [x1 x2 y1 y2] + # coarsen + # buffer + # align [x y] + # x_shift + # z_shift + T|F # negate_z + +``topo_path`` + Absolute path to the topography file, single-quoted. + +``topo_type`` + Integer format code: ``2`` (one-value-per-line ASCII), ``3`` (row-major + ASCII), ``4`` (NetCDF), ``5`` (GeoTIFF). Negative values (e.g. ``-2``) + apply an additional sign flip to Z on load (Fortran convention). + +``crop_extent`` + Four floats: ``x1 x2 y1 y2`` defining a bounding box. + Sentinel ``"0. 0. 0. 0."`` means no cropping (sentinel is safe because + a valid extent requires ``x1 < x2`` and ``y1 < y2``). + +``coarsen`` + Integer stride factor for subsampling (``1`` = no coarsening). + +``buffer`` + Integer number of grid points to retain outside the crop region on each + side (``0`` = no buffer). + +``align`` + Two floats: ``x_align y_align`` alignment target for subsampled grids. + Sentinel ``"0. 0."`` means no alignment constraint. + +``x_shift`` + Float added to all x coordinates after loading (``0`` = no shift). + +``z_shift`` + Float added to all non-missing elevation values (``0`` = no shift). + +``negate_z`` + ``T`` or ``F``. If ``T``, all elevation values are negated after loading, + independently of the ``topo_type`` sign convention. + +.. _netcdf_descriptor_block: + +NetCDF Descriptor Block (type 4 only) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For ``topo_type = 4`` files, the 9-line block is followed immediately by a +CF descriptor block. The descriptor is a sequence of ``key=value`` lines +terminated by a blank line, parsed by Fortran's ``read_netcdf_descriptor``. +Example:: + + var_name=elevation + lon_name=lon + lat_name=lat + lat_order=S_to_N + dim_order=lat,lon + lon_convention=180 + lon_offset=0.0 + source_units=m + fill_action=abort + +``lon_offset`` + Scalar added by Fortran to file longitudes to convert them to domain + coordinates: ``x_domain = x_file + lon_offset``. For files using the + ``[0, 360]`` convention against a ``[-180, 180]`` domain, this is + ``-360.0``. Use + :meth:`TopoInspector.topo_entries() + ` to compute + the correct value automatically. + + +Annotated Examples +------------------ + +Type-2 ASCII file, no preprocessing +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +:: + + 99999.0 # replace no_data_value in topofile + 0 # (Type topography specification) + 1 # number of topo files + + '/data/etopo1.tt2' # topo_path + 2 # topo_type + 0. 0. 0. 0. # crop_extent [x1 x2 y1 y2] + 1 # coarsen + 0 # buffer + 0. 0. # align [x y] + 0 # x_shift + 0 # z_shift + F # negate_z + +Type-4 NetCDF file, with crop and negate +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +:: + + 99999.0 # replace no_data_value in topofile + 0 # (Type topography specification) + 1 # number of topo files + + '/data/gebco.nc' # topo_path + 4 # topo_type + -100. -60. 10. 50. # crop_extent [x1 x2 y1 y2] + 1 # coarsen + 2 # buffer + 0. 0. # align [x y] + 0 # x_shift + 0 # z_shift + T # negate_z + + var_name=elevation + lon_name=lon + lat_name=lat + lat_order=S_to_N + dim_order=lat,lon + lon_convention=180 + lon_offset=0.0 + source_units=m + fill_action=abort + + +.. _priority_convention: + +Priority Convention +------------------- + +The first file block in ``topo.data`` is assigned **rank 1** by Fortran +(highest priority in overlap resolution). Fortran stores this mapping in +``mtopoorder``: ``mtopoorder(1) = 1`` means file-index 1 → rank 1. + +Python's :meth:`~clawpack.geoclaw.data.TopographyData._compute_priority_order` +sorts files by cell area ascending (finest = smallest ``dx * dy``) before +writing, so the finest file is automatically written first. + +To override this sort and use your own ordering, set +``TopographyData.override_order = True``. When ``True``, you are responsible +for placing the finest file first in ``topofiles``. + +.. warning:: + ``override_order`` is a Python-only attribute. It is **not** written to + ``topo.data`` and has no Fortran-side counterpart. + + +Sentinel Values +--------------- + ++----------------+--------------------+--------------------------------------+ +| Parameter | Sentinel value | Meaning | ++================+====================+======================================+ +| ``crop_extent``| ``0. 0. 0. 0.`` | No cropping; use full file domain. | ++----------------+--------------------+--------------------------------------+ +| ``align`` | ``0. 0.`` | No alignment constraint. | ++----------------+--------------------+--------------------------------------+ +| ``negate_z`` | ``F`` | No sign flip. | ++----------------+--------------------+--------------------------------------+ + +Sentinel values are geometrically invalid as real parameters (a valid extent +requires ``x1 < x2``; a valid alignment is non-zero in typical cases) so +Fortran can unambiguously distinguish "not set" from a genuine value. + + +Deprecated Formats +------------------ + +In older setrun.py files, ``topofiles`` entries were plain Python lists or +dicts. These are still accepted by +:class:`~clawpack.geoclaw.data.TopographyData` but emit a +``DeprecationWarning``. + +**Old list format (deprecated):** + +.. code-block:: python + + rundata.topo_data.topofiles.append([2, '/path/to/topo.tt2']) + # or the older 6-element form (level/time info is ignored): + rundata.topo_data.topofiles.append([2, 1, 5, 0.0, 1e10, '/path/to/topo.tt2']) + +**Old dict format (deprecated):** + +.. code-block:: python + + rundata.topo_data.topofiles.append({ + 'topo_type': 2, + 'topo_path': '/path/to/topo.tt2', + 'extent': [-100., -60., 10., 50.], # note: 'extent' → crop_extent + 'coarsen': 2, + }) + +Note: the dict key ``'extent'`` maps to ``Topography.crop_extent``, **not** +to ``Topography.extent`` (which is a read-only property returning loaded-data +bounds). + +**Recommended form:** + +.. code-block:: python + + from clawpack.geoclaw.topotools import Topography + + t = Topography() + t.path = '/path/to/topo.tt2' + t.topo_type = 2 + t.crop_extent = [-100., -60., 10., 50.] + t.coarsen = 2 + rundata.topo_data.topofiles.append(t) + +**topo_type=1 (deprecated):** + +Type-1 files (``x y z`` ASCII, one point per line) are still readable with a +``DeprecationWarning``, but preprocessing attributes (``crop_extent``, +``coarsen``, etc.) raise ``NotImplementedError``. To convert:: + + t = Topography() + t.read('old_file.tt1', topo_type=1) # DeprecationWarning + t.write('new_file.tt2', topo_type=2) # convert to type 2 + +For genuinely unstructured (scattered) point data, use +``scipy.interpolate``, GMT, or a similar tool to grid the data before use. diff --git a/doc/topotools.rst b/doc/topotools.rst index bf6ffb7..389a31b 100644 --- a/doc/topotools.rst +++ b/doc/topotools.rst @@ -4,6 +4,80 @@ 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. + + +.. 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 and must be gridded externally (e.g. ``scipy.interpolate``, GMT) + before use in GeoClaw. + .. toctree:: :maxdepth: 1 From ff47df630182bda20ac07bf862471e3e18248f69 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Tue, 16 Jun 2026 10:38:55 -0400 Subject: [PATCH 10/16] Add docs for ordering of topo and dtopo --- doc/changes_to_master.rst | 11 ++++++----- doc/topo_order.rst | 12 +++++++----- doc/topodata_format.rst | 17 ++++++++++------- 3 files changed, 23 insertions(+), 17 deletions(-) diff --git a/doc/changes_to_master.rst b/doc/changes_to_master.rst index d174d42..688cd31 100644 --- a/doc/changes_to_master.rst +++ b/doc/changes_to_master.rst @@ -86,11 +86,12 @@ Changes to geoclaw - **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. The first file listed in ``topo.data`` is assigned rank 1 - (highest priority) by Fortran; no Fortran-side sorting occurs. - ``rundata.topo_data.override_order = True`` preserves user-specified list - order; when used, the finest (highest-resolution) file should be listed - first. See :ref:`topo_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) diff --git a/doc/topo_order.rst b/doc/topo_order.rst index cbc740b..311394f 100644 --- a/doc/topo_order.rst +++ b/doc/topo_order.rst @@ -23,8 +23,10 @@ then the finer one is used. This ordering is computed entirely in Python by ``TopographyData._compute_priority_order`` and written into ``topo.data`` -before the Fortran run. The first file listed in ``topo.data`` is assigned -rank 1 (highest priority) by Fortran. No Fortran-side sorting occurs. +before the Fortran run. 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. No +Fortran-side sorting occurs. This usually works fine but may not always be what is wanted. In particular there may be two DEMs at essentially the same resolution, but one is a more @@ -36,7 +38,7 @@ resolutions are essentially identical, `dx*dy` maybe be slightly different (perhaps at the roundoff level) and so it could be that the older one is viewed as "finer" and would be used preferentially. Even if `dx*dy` is exactly the same, it is not obvious which one would be used in GeoClaw. -(The one listed first in `rundata.topo_data.topofiles` +(The one listed last in `rundata.topo_data.topofiles` is given higher priority in this case, with a relative tolerance of 0.001 used to judge whether the resolutions are the same). @@ -55,8 +57,8 @@ If `setrun.py` specifies:: then the preference order is no longer based on the computed resolution. Instead the list order in `rundata.topo_data.topofiles` is used as-is: the -first file in the list is given the highest priority (rank 1). Normally you -should list the finest (highest-resolution) DEM first and the coarsest last, +last file in the list is given the highest priority (rank 1). Normally you +should list the coarsest DEM first and the finest (highest-resolution) last, to reproduce the default behavior. This allows precise specification of the order of preference when automatic sorting by resolution is not appropriate. diff --git a/doc/topodata_format.rst b/doc/topodata_format.rst index 08ac60c..c2d7569 100644 --- a/doc/topodata_format.rst +++ b/doc/topodata_format.rst @@ -54,8 +54,8 @@ Global Header (3 lines) Per-File Block (9 lines, all types) ------------------------------------ -One block follows for each topography file, in priority order (finest -resolution first; see :ref:`priority_convention` below). +One block follows for each topography file, in priority order (coarsest +resolution first, finest last; see :ref:`priority_convention` below). :: @@ -190,17 +190,20 @@ Type-4 NetCDF file, with crop and negate Priority Convention ------------------- -The first file block in ``topo.data`` is assigned **rank 1** by Fortran +The **last** file block in ``topo.data`` is assigned **rank 1** by Fortran (highest priority in overlap resolution). Fortran stores this mapping in -``mtopoorder``: ``mtopoorder(1) = 1`` means file-index 1 → rank 1. +``mtopoorder`` using the reversed assignment ``mtopoorder(i) = mtopofiles+1-i``, +so the last file listed maps to rank 1. Python's :meth:`~clawpack.geoclaw.data.TopographyData._compute_priority_order` -sorts files by cell area ascending (finest = smallest ``dx * dy``) before -writing, so the finest file is automatically written first. +sorts files by cell area descending (coarsest = largest ``dx * dy`` first, +finest last) before writing, so the finest file is automatically written +last. This matches the traditional GeoClaw convention of listing topography +files from coarsest to finest. To override this sort and use your own ordering, set ``TopographyData.override_order = True``. When ``True``, you are responsible -for placing the finest file first in ``topofiles``. +for placing the finest file last in ``topofiles``. .. warning:: ``override_order`` is a Python-only attribute. It is **not** written to From 97d3aeedd312daa21caf4e48678f5534a3f7595e Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Tue, 16 Jun 2026 10:39:35 -0400 Subject: [PATCH 11/16] Add pyclaw.log to .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index dd68408..26f2e21 100644 --- a/.gitignore +++ b/.gitignore @@ -63,3 +63,4 @@ Thumbs.db .uv/ .pixi/ conda-meta/ +doc/pyclaw.log From 6426c15f34acf4742e8e2ff4e0f7036873dd5535 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Tue, 16 Jun 2026 20:37:14 -0400 Subject: [PATCH 12/16] Update with unstructured info --- doc/topodata_format.rst | 7 ++++-- doc/topotools.rst | 49 +++++++++++++++++++++++++++++++++++++++-- 2 files changed, 52 insertions(+), 4 deletions(-) diff --git a/doc/topodata_format.rst b/doc/topodata_format.rst index c2d7569..5ccf3c2 100644 --- a/doc/topodata_format.rst +++ b/doc/topodata_format.rst @@ -282,5 +282,8 @@ Type-1 files (``x y z`` ASCII, one point per line) are still readable with a t.read('old_file.tt1', topo_type=1) # DeprecationWarning t.write('new_file.tt2', topo_type=2) # convert to type 2 -For genuinely unstructured (scattered) point data, use -``scipy.interpolate``, GMT, or a similar tool to grid the data before use. +For genuinely unstructured (scattered) point data, grid it onto a regular, +logically rectangular grid before use -- either with an external tool +(``scipy.interpolate``, GMT) or with +:meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` (see +:ref:`topotools`). diff --git a/doc/topotools.rst b/doc/topotools.rst index 389a31b..57a44b3 100644 --- a/doc/topotools.rst +++ b/doc/topotools.rst @@ -75,8 +75,53 @@ to determine file resolution without loading elevation data. t.write('new.tt2', topo_type=2) # save as type 2 Genuinely unstructured (scattered) point data cannot be converted this - way and must be gridded externally (e.g. ``scipy.interpolate``, GMT) - before use in GeoClaw. + 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. .. toctree:: From e9ac520ceea8e6c60b99dd601d33b89851763ccc Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Fri, 19 Jun 2026 09:58:05 -0400 Subject: [PATCH 13/16] Add docs for datum and formatting guidance --- doc/topodata_format.rst | 56 ++++++++++++++++++++++++++--------------- doc/topotools.rst | 20 +++++++++++++++ 2 files changed, 56 insertions(+), 20 deletions(-) diff --git a/doc/topodata_format.rst b/doc/topodata_format.rst index 5ccf3c2..e849578 100644 --- a/doc/topodata_format.rst +++ b/doc/topodata_format.rst @@ -51,7 +51,7 @@ Global Header (3 lines) sorts the files before writing but is **not** written to ``topo.data``. -Per-File Block (9 lines, all types) +Per-File Block (10 lines, all types) ------------------------------------ One block follows for each topography file, in priority order (coarsest @@ -66,6 +66,7 @@ resolution first, finest last; see :ref:`priority_convention` below). # buffer # align [x y] # x_shift + # y_shift # z_shift T|F # negate_z @@ -95,6 +96,11 @@ resolution first, finest last; see :ref:`priority_convention` below). ``x_shift`` Float added to all x coordinates after loading (``0`` = no shift). + A registration offset: ``x_domain = x_file + x_shift``. + +``y_shift`` + Float added to all y coordinates after loading (``0`` = no shift), + the y-direction counterpart of ``x_shift``. ``z_shift`` Float added to all non-missing elevation values (``0`` = no shift). @@ -108,26 +114,36 @@ resolution first, finest last; see :ref:`priority_convention` below). NetCDF Descriptor Block (type 4 only) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -For ``topo_type = 4`` files, the 9-line block is followed immediately by a +For ``topo_type = 4`` files, the 10-line block is followed immediately by a CF descriptor block. The descriptor is a sequence of ``key=value`` lines terminated by a blank line, parsed by Fortran's ``read_netcdf_descriptor``. +The coordinate keys are projection-agnostic (``x_name``/``y_name``, not +``lon``/``lat``); only the ``lon_wrap_offset`` longitude wrap is geographic. Example:: var_name=elevation - lon_name=lon - lat_name=lat - lat_order=S_to_N - dim_order=lat,lon - lon_convention=180 - lon_offset=0.0 - source_units=m + x_name=lon + y_name=lat + y_increasing=True + dim_order=y,x + lon_wrap_offset=0.0 fill_action=abort -``lon_offset`` +``x_name`` / ``y_name`` + Names of the x and y coordinate variables in the file (auto-detected by + the inspector via CF ``axis``/``standard_name``/common names). + +``y_increasing`` + ``True`` if the y coordinate increases with array index, else ``False``. + Informational: the y-axis direction is re-detected from the coordinates + at read time. + +``lon_wrap_offset`` Scalar added by Fortran to file longitudes to convert them to domain - coordinates: ``x_domain = x_file + lon_offset``. For files using the - ``[0, 360]`` convention against a ``[-180, 180]`` domain, this is - ``-360.0``. Use + coordinates: ``x_domain = x_file + lon_wrap_offset``. For a geographic + file using the ``[0, 360]`` convention against a ``[-180, 180]`` domain, + this is ``-360.0``; for a non-geographic (projected) x axis it is ``0.0`` + (the wrap is skipped). Use :meth:`TopoInspector.topo_entries() ` to compute the correct value automatically. @@ -152,6 +168,7 @@ Type-2 ASCII file, no preprocessing 0 # buffer 0. 0. # align [x y] 0 # x_shift + 0 # y_shift 0 # z_shift F # negate_z @@ -171,17 +188,16 @@ Type-4 NetCDF file, with crop and negate 2 # buffer 0. 0. # align [x y] 0 # x_shift + 0 # y_shift 0 # z_shift T # negate_z var_name=elevation - lon_name=lon - lat_name=lat - lat_order=S_to_N - dim_order=lat,lon - lon_convention=180 - lon_offset=0.0 - source_units=m + x_name=lon + y_name=lat + y_increasing=True + dim_order=y,x + lon_wrap_offset=0.0 fill_action=abort diff --git a/doc/topotools.rst b/doc/topotools.rst index 57a44b3..1bff5a1 100644 --- a/doc/topotools.rst +++ b/doc/topotools.rst @@ -62,6 +62,26 @@ 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 From 8df638cfd7a165d3da4ff9a7f660d89f30ee5cf0 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Sat, 11 Jul 2026 15:57:43 -0400 Subject: [PATCH 14/16] Update docs for netCDF dtopo support Also updates other topo, dtopo, and met netCDF input file support. --- doc/changes_to_master.rst | 34 +++++++ doc/dtopo.rst | 14 ++- doc/netcdf.rst | 186 +++++++++++++++++++++++++++--------- doc/netcdf_utils_module.rst | 10 +- 4 files changed, 192 insertions(+), 52 deletions(-) diff --git a/doc/changes_to_master.rst b/doc/changes_to_master.rst index 688cd31..256d488 100644 --- a/doc/changes_to_master.rst +++ b/doc/changes_to_master.rst @@ -21,6 +21,19 @@ the menu on the left hand side of this page. Changes that are not backward compatible ---------------------------------------- +- **NetCDF input units are now required and enforced (geoclaw).** NetCDF + topography and dtopo elevation must declare CF ``units`` of meters, and + gridded meteorological forcing must be wind ``m/s`` / pressure ``Pa``. A + missing ``units`` attribute, or a non-contract unit (e.g. ``km``, ``hPa``, + ``mbar``, ``knots``), now raises a ``ValueError`` instead of being assumed + or silently converted -- GeoClaw does not convert units on the read path. + Pre-convert to the contract units, or pass ``assume_units`` + (``TopoInspector`` / ``Topography.read`` via ``nc_params``, or + ``MetInspector(assume_units=True)``) for a file known to be in contract + units but lacking the attribute. Met NetCDF forcing additionally requires + a CF *absolute* datetime time axis; a bare numeric time axis is rejected. + See :ref:`netcdf_input`. + General changes --------------- @@ -107,6 +120,27 @@ Changes to geoclaw 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. + See `geoclaw diffs `_ diff --git a/doc/dtopo.rst b/doc/dtopo.rst index 347ada0..7ad42d3 100644 --- a/doc/dtopo.rst +++ b/doc/dtopo.rst @@ -115,7 +115,8 @@ on a uniform grid at a sequence of one or more times. dtopo file formats ------------------ -Currently two formats are supported for this file: +Three formats are supported for this file: the two ASCII types below and the +CF-compliant NetCDF ``dtopotype=4`` (see :ref:`netcdf_input`): **dtopotype=1:** @@ -135,6 +136,17 @@ Currently two formats are supported for this file: dx, dy*, and *dt*. These are followed by *mt* sets of *my* lines, each line containing *mx* values of *dz*. + **dtopotype=4:** + + A CF-compliant NetCDF format. ``dtopotools.DTopography`` reads and writes + this type (``dtopo_type=4``) and GeoClaw reads it directly at runtime. + The deformation is in meters; the time axis may be a bare CF duration + (``units = "seconds"``, simulation-relative -- the default on write) or a + CF datetime axis (``units = "seconds since "``, written when + ``DTopography.time_reference`` is set to an epoch such as the earthquake + origin time), and is interpreted on read by its CF ``units``. See + :ref:`netcdf_input` for details, including the ``dz_dtype`` write option. + Here is an example header for a dtopofile:: 721 mx diff --git a/doc/netcdf.rst b/doc/netcdf.rst index 173fd1c..8949a94 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -1,23 +1,10 @@ -.. toctree:: - :maxdepth: 2 - - user_guide - developer_reference - descriptor_format - unit_contract - coordinate_normalization - technical_debt - adding_met_field - test_structure - next_steps - .. _netcdf_input: GeoClaw NetCDF Input System =========================== This document covers the NetCDF input pipeline introduced in the -``refactor-netcdf-support`` PR (VERSION ?5.15?). It has two sections: a user +``refactor-netcdf-support`` PR (first released after v5.14.0). 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. @@ -46,13 +33,24 @@ What GeoClaw handles automatically resolved automatically. If a fill value is found within your simulation domain, GeoClaw will abort with an error -- this is intentional, as missing bathymetry is a silent correctness hazard. +- File extension: the NetCDF backend engine is selected explicitly, so a + valid NetCDF file opens even when its name uses a non-standard + extension (e.g. ``.dtt3``) that xarray's extension-based engine + guessing would otherwise fail to recognize. What your file must provide ^^^^^^^^^^^^^^^^^^^^^^^^^^^ - A 2D variable containing elevation in **meters** (positive up, - negative for ocean). If your file uses feet or another unit, use - ``CFNormalizer`` to convert before running (see below). + negative for ocean), carrying a CF ``units`` attribute equal to meters + (e.g. ``units = "m"``). Units are **required** and never assumed: a + variable with no ``units`` attribute, or with a non-meter unit (even a + recognized one such as ``km``), is **rejected** with an error rather + than being silently misread -- GeoClaw does not convert units on the + read path. Pre-convert non-meter data to meters first. For a file + that is genuinely in meters but merely omits the attribute, you can opt + in explicitly by passing ``assume_units='m'`` to ``TopoInspector`` (or + ``nc_params={'assume_units': 'm'}`` to ``Topography.read``). - 1D coordinate variables for latitude and longitude with recognizable names (``lat``/``latitude``/``y`` and ``lon``/``longitude``/``x`` are all detected). Curvilinear (2D coordinate) grids are not currently @@ -142,6 +140,56 @@ If you are unsure whether your file will be read correctly, the attributes, and resolves ``_FillValue``/``missing_value`` conflicts. It does not resample or reproject. +Writing NetCDF topo files +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +A ``topotools.Topography`` object can be written as a CF-compliant NetCDF +file (``topo_type=4``) that GeoClaw reads back through the same descriptor +mechanism:: + + topo.write('bathy.nc', topo_type=4) # elevation float32 + topo.write('bathy.nc', topo_type=4, z_dtype='float64') # full precision + +The elevation variable is written with ``units = "m"``. It is stored on +disk as ``float32`` by default -- sub-millimeter precision for Earth +topography (``abs(Z) < 10000`` m) at half the file size -- and +``z_dtype='float64'`` selects full double precision. + +-------------- + +Seafloor deformation (dtopo) from NetCDF +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Time-dependent seafloor deformation (dtopo) can also be read from and +written to CF-compliant NetCDF using ``dtopo_type=4`` (see :ref:`dtopo`): + +.. code:: python + + from clawpack.geoclaw import dtopotools + dtopo = dtopotools.DTopography('deformation.nc', dtopo_type=4) # read + dtopo.write('out.nc', dtopo_type=4) # write + dtopo.write('out.nc', dtopo_type=4, dz_dtype='float64') # full precision + +As with topography, the deformation is in meters and stored ``float32`` by +default (pass ``dz_dtype='float64'`` for full double precision). + +**Time axis.** By default the time coordinate is written as a bare CF +duration, ``units = "seconds"``, holding simulation-relative times. If you +set the optional ``time_reference`` attribute to a real-world epoch (for +example the earthquake origin time), the file is instead written with a CF +datetime axis, ``units = "seconds since "``:: + + dtopo.time_reference = '2011-03-11T05:46:00' + dtopo.write('out.nc', dtopo_type=4) + +A datetime axis is more interoperable (a plain ``xarray.open_dataset`` or a +GIS tool decodes it to real timestamps) and round-trips back to the same +simulation-relative times. On read, the time axis is interpreted using its +CF ``units``, so an axis in ``"minutes"``, ``"hours"``, or a +``" since "`` datetime is scaled to seconds correctly rather +than assumed to already be in seconds. Unlike met forcing, a bare +``"seconds"`` (relative) dtopo axis is allowed. + -------------- Storm surge met forcing from NetCDF @@ -164,21 +212,31 @@ curvilinear grid and string-encoded time axis). Required variables and units ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -GeoClaw expects the following in the NetCDF file **after any unit -conversion** (conversion happens automatically during Python -preprocessing): +GeoClaw expects the following variables in the NetCDF file, already in +the contract units shown (units are verified, not converted): -============================= ================================ +============================= ================================== Variable Unit -============================= ================================ +============================= ================================== Wind (u-component, eastward) m/s Wind (v-component, northward) m/s Surface pressure Pa -Time seconds from user-defined offset -============================= ================================ - -If your file uses hPa/mbar for pressure or knots for wind, -``MetInspector`` will convert automatically. +Time seconds since a reference date +============================= ================================== + +Units are **required** and never assumed. A variable with no ``units`` +attribute, or in a non-contract unit (e.g. ``hPa``/``mbar`` for pressure +or ``knots`` for wind), is **rejected** with an error rather than being +converted -- GeoClaw does not convert units on the read path. Pre-convert +such a file to ``m/s`` / ``Pa`` first (or, for a file known to be in +contract units but lacking the attributes, pass ``assume_units=True`` to +``MetInspector``). The time coordinate must be an **absolute time in +seconds** -- CF ``units`` of the form ``"seconds since "`` (e.g. +``"seconds since 2020-01-01"``). Fortran reads the raw time values as +integer seconds and does **not** scale by the CF unit, so a file whose time +is in ``hours``/``days`` since a reference (such as a raw ERA5 file) must be +converted to ``"seconds since ..."`` first; a bare numeric/duration time +axis with no reference date is rejected outright (see `Time handling`_ below). Registering a NetCDF storm file ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -311,6 +369,17 @@ days→seconds) is performed in Fortran: the descriptor stores ``write_data`` passes ``time_reference=storm.time_offset`` to ``MetInspector`` automatically; no additional user action is needed. +The met time axis must be an absolute time in **seconds** -- CF ``units`` of +the form ``"seconds since "`` (e.g. ``"seconds since 2020-01-01"``), +which xarray decodes to datetimes for the offset computation above. Because +Fortran reads the raw time values as integer seconds and does **not** scale +by the CF unit, a file storing time in ``hours``/``days`` since a reference +(such as a raw ERA5 file) must be converted to ``"seconds since ..."`` first. +A bare numeric/duration axis with no reference date is **rejected** outright. +This differs from dtopo NetCDF, where Python collapses the time axis to +``(t0, dt)`` in seconds and therefore any CF duration or datetime axis is +accepted (see `Seafloor deformation (dtopo) from NetCDF`_). + Longitude convention ^^^^^^^^^^^^^^^^^^^^^^ @@ -331,7 +400,8 @@ Architecture overview The system has a strict Python/Fortran split: **Python** handles: file inspection, CF attribute parsing, coordinate -convention detection, fill value resolution, unit conversion, +convention detection, fill value resolution, unit verification +(non-contract units are rejected, not converted), ``nc_time_offset`` computation (elapsed seconds from ``storm.time_offset`` to the first record), crop bound validation, and descriptor writing. @@ -360,13 +430,13 @@ Class hierarchy (``netcdf_utils.py``) TopoInspector(NetCDFInspector) - detect fill values within crop region (hard error) - - verify and convert units to contract + - verify units against the contract (reject non-contract; no conversion) - no multi-file coverage logic (Fortran handles compositing) MetInspector(NetCDFInspector) - check wind_u, wind_v, pressure present and on same grid/time axis - - convert units to contract - - decode CF time to seconds from offset + - verify units against the contract (reject non-contract; no conversion) + - decode CF datetime axis to seconds from offset (reject bare-numeric) - detect ensemble/member dimensions (hard error if non-singleton) - partial domain coverage is allowed (Fortran fills edges) @@ -389,10 +459,10 @@ Topo (lines in topo.data after topo_type) :: var_name = z - lon_name = longitude - lat_name = latitude - lon_offset = 0.0 - lat_order = S_to_N + x_name = longitude + y_name = latitude + lon_wrap_offset = 0.0 + y_increasing = True dim_order = lat,lon fill_value = -9999.0 fill_action = abort @@ -401,29 +471,50 @@ Topo (lines in topo.data after topo_type) ``fill_action = abort`` is the only supported value for topography. ``crop_bounds`` is omitted if no crop is specified. +dtopo (lines in dtopo.data after the per-file block) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +:: + + var_name = dz + x_name = longitude + y_name = latitude + time_name = time + lon_wrap_offset = 0.0 + y_increasing = True + dim_order = time,lat,lon + t0 = 0.0 + dt = 10.0 + +``t0`` and ``dt`` give the first time and the uniform time step in simulation +**seconds**; Python collapses the CF time axis to these (scaling by the file's +CF ``units``), so Fortran never parses CF time for dtopo. ``time_name`` is the +name of the time coordinate in the file. + Met forcing (body of \*.storm after format header) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ :: - # format: netcdf &file_info - source_file = /path/to/forcing.nc - lon_name = longitude - lat_name = latitude + x_name = longitude + y_name = latitude time_name = valid_time dim_order = time,lat,lon - lon_convention = 360 - lat_order = S_to_N + lon_wrap = 360 + y_increasing = True fill_value = -9999.0 fill_action = warn time_offset = 0.0 - crop_bounds = -100.0 -80.0 20.0 35.0 / &variable_info var_name=u10 geoclaw_role=wind_u / &variable_info var_name=v10 geoclaw_role=wind_v / &variable_info var_name=msl geoclaw_role=pressure / +(``lon_wrap`` is written only for a geographic longitude axis; ``fill_value`` +only when the file declares one; ``crop_bounds`` is carried by the top-level +``crop_extent`` line rather than the descriptor.) + The ``&variable_info`` blocks use a manually parsed key=value format (not true repeated Fortran namelists) for compiler portability. Fortran reads them in a loop until EOF. Adding new roles (precipitation, @@ -444,16 +535,19 @@ Defined in ``units.py``: "time": "s", } -All conversion happens in Python before the descriptor is written. -Fortran trusts the descriptor and never checks units. If you add a new -``geoclaw_role``, add its contract unit here first. +Python verifies units against this contract and **rejects** files whose +variables are in non-contract units or lack a ``units`` attribute; no unit +conversion is performed on the read path. Fortran trusts the descriptor +and never checks units. If you add a new ``geoclaw_role``, add its contract +unit here first. Fortran coordinate normalization ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Fortran applies coordinate normalization via index arithmetic only -- no -arrays are duplicated. The descriptor provides ``lon_convention`` and -``lat_order``; the reader uses these to compute correct indices when +arrays are duplicated. The descriptor provides the longitude-wrap field +(``lon_wrap`` for met, ``lon_wrap_offset`` for topo) and the latitude order +(``y_increasing``); the reader uses these to compute correct indices when calling ``nf90_get_var``. For crop bounds, ``start`` and ``count`` are computed from the coordinate arrays at runtime. diff --git a/doc/netcdf_utils_module.rst b/doc/netcdf_utils_module.rst index a760778..359d036 100644 --- a/doc/netcdf_utils_module.rst +++ b/doc/netcdf_utils_module.rst @@ -11,9 +11,9 @@ netcdf_utils module for NetCDF input The ``netcdf_utils`` module provides the Python layer for reading NetCDF topography and meteorological forcing files in GeoClaw. It handles CF -attribute parsing, coordinate convention detection, unit conversion, and -descriptor writing so that Fortran only needs to open the file and read -pre-validated indices. +attribute parsing, coordinate convention detection, unit verification +(non-contract units are rejected, not converted), and descriptor writing so +that Fortran only needs to open the file and read pre-validated indices. The main classes are: @@ -22,8 +22,8 @@ The main classes are: - :class:`~clawpack.geoclaw.netcdf_utils.TopoInspector` — topography subclass; fill-value checking within the crop region, unit enforcement. - :class:`~clawpack.geoclaw.netcdf_utils.MetInspector` — meteorological - subclass; wind/pressure variable discovery, CF time decoding, unit - conversion. + subclass; wind/pressure variable discovery, CF datetime decoding, unit + verification (non-contract or missing units are rejected). - :class:`~clawpack.geoclaw.netcdf_utils.CFNormalizer` — adds or repairs CF attributes in place without modifying data values. - :class:`~clawpack.geoclaw.netcdf_utils.DescriptorWriter` — writes the From b0ead5c9992277f42484a71d9c0f59226d485c3c Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Sun, 12 Jul 2026 13:39:51 -0400 Subject: [PATCH 15/16] Update docs for non-contract unit handling Addresses a slight regression in capability related to units that are incorrect or missing when it goes to Fortran. --- doc/changes_to_master.rst | 31 +++--- doc/netcdf.rst | 181 +++++++++++++++++++++++------------- doc/netcdf_utils_module.rst | 11 ++- 3 files changed, 140 insertions(+), 83 deletions(-) diff --git a/doc/changes_to_master.rst b/doc/changes_to_master.rst index 256d488..a63f05e 100644 --- a/doc/changes_to_master.rst +++ b/doc/changes_to_master.rst @@ -21,18 +21,25 @@ the menu on the left hand side of this page. Changes that are not backward compatible ---------------------------------------- -- **NetCDF input units are now required and enforced (geoclaw).** NetCDF - topography and dtopo elevation must declare CF ``units`` of meters, and - gridded meteorological forcing must be wind ``m/s`` / pressure ``Pa``. A - missing ``units`` attribute, or a non-contract unit (e.g. ``km``, ``hPa``, - ``mbar``, ``knots``), now raises a ``ValueError`` instead of being assumed - or silently converted -- GeoClaw does not convert units on the read path. - Pre-convert to the contract units, or pass ``assume_units`` - (``TopoInspector`` / ``Topography.read`` via ``nc_params``, or - ``MetInspector(assume_units=True)``) for a file known to be in contract - units but lacking the attribute. Met NetCDF forcing additionally requires - a CF *absolute* datetime time axis; a bare numeric time axis is rejected. - See :ref:`netcdf_input`. +- **NetCDF input units are now required, and recognized units are converted + automatically (geoclaw).** NetCDF input variables must declare a CF + ``units`` attribute. A variable already in the contract unit (topo/dtopo + ``m``, wind ``m/s``, pressure ``Pa``) passes through unchanged; a + *recognized* non-contract unit (e.g. ``km``/``cm``, ``hPa``/``mbar``, + ``knots``, or a ``"hours since ..."`` time axis) is now **converted + automatically** -- Python resolves it to a multiplicative scale factor that + Fortran applies on read. A **missing** ``units`` attribute or an + **unrecognized** unit still raises a ``ValueError`` rather than being + silently misread (met forcing falls back to the storm format's documented + unit when the attribute is absent, e.g. NWS13/OWI pressure ``mbar``). Pass + ``assume_units`` (``TopoInspector`` / ``Topography.read`` via ``nc_params``, + or ``MetInspector(assume_units=True)``) for a file that omits the attribute. + A resolved field is also magnitude-checked: a pressure field mislabeled + ``Pa`` but really ``hPa``/``mbar`` is auto-corrected with a warning, and a + physically implausible field raises (bypass with ``skip_sanity_check``). + Met NetCDF forcing requires a CF *absolute* datetime time axis; a bare + numeric time axis with no reference date is rejected. See + :ref:`netcdf_input`. General changes diff --git a/doc/netcdf.rst b/doc/netcdf.rst index 8949a94..dcc4947 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -24,6 +24,12 @@ NOAA coastal DEMs, and any file produced by standard GIS tools. What GeoClaw handles automatically ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +- Units: a recognized non-contract unit (e.g. ``km``/``cm`` for elevation, + ``hPa``/``mbar`` for pressure, ``knots`` for wind) is converted + automatically. Python resolves the unit to a single multiplicative scale + factor; Fortran applies it on read (in-memory Python reads convert + directly). A missing or unrecognized unit is still an error -- see the + per-variable notes below. - Longitude conventions: both [-180, 180] and [0, 360] are detected and normalized at runtime. You do not need to preprocess the file. - Latitude ordering: both S-to-N and N-to-S are handled correctly. @@ -41,16 +47,18 @@ What GeoClaw handles automatically What your file must provide ^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- A 2D variable containing elevation in **meters** (positive up, - negative for ocean), carrying a CF ``units`` attribute equal to meters - (e.g. ``units = "m"``). Units are **required** and never assumed: a - variable with no ``units`` attribute, or with a non-meter unit (even a - recognized one such as ``km``), is **rejected** with an error rather - than being silently misread -- GeoClaw does not convert units on the - read path. Pre-convert non-meter data to meters first. For a file - that is genuinely in meters but merely omits the attribute, you can opt - in explicitly by passing ``assume_units='m'`` to ``TopoInspector`` (or - ``nc_params={'assume_units': 'm'}`` to ``Topography.read``). +- A 2D variable containing elevation (positive up, negative for ocean), + carrying a CF ``units`` attribute. Meters (e.g. ``units = "m"``) is the + contract unit; a variable in a **recognized** non-meter unit such as + ``km`` or ``cm`` is **converted automatically** on read (a warning is + emitted noting the conversion). Units are still **required** and never + assumed: a variable with **no** ``units`` attribute, or one whose unit is + unrecognized, is rejected with an error rather than being silently + misread. For a file that is genuinely in meters but merely omits the + attribute, opt in explicitly by passing ``assume_units='m'`` to + ``TopoInspector`` (or ``nc_params={'assume_units': 'm'}`` to + ``Topography.read``); ``assume_units`` also accepts a non-meter unit + (e.g. ``'km'``), which is then converted. - 1D coordinate variables for latitude and longitude with recognizable names (``lat``/``latitude``/``y`` and ``lon``/``longitude``/``x`` are all detected). Curvilinear (2D coordinate) grids are not currently @@ -170,8 +178,11 @@ written to CF-compliant NetCDF using ``dtopo_type=4`` (see :ref:`dtopo`): dtopo.write('out.nc', dtopo_type=4) # write dtopo.write('out.nc', dtopo_type=4, dz_dtype='float64') # full precision -As with topography, the deformation is in meters and stored ``float32`` by -default (pass ``dz_dtype='float64'`` for full double precision). +As with topography, the deformation contract unit is meters and it is stored +``float32`` by default (pass ``dz_dtype='float64'`` for full double +precision). A recognized non-meter unit (e.g. ``km``) is converted +automatically on read, exactly as for topography; a missing or unrecognized +unit raises. (ASCII dtopo types remain meters-implied, unchanged.) **Time axis.** By default the time coordinate is written as a bare CF duration, ``units = "seconds"``, holding simulation-relative times. If you @@ -212,31 +223,45 @@ curvilinear grid and string-encoded time axis). Required variables and units ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -GeoClaw expects the following variables in the NetCDF file, already in -the contract units shown (units are verified, not converted): +GeoClaw uses the following variables, whose **contract units** are shown. +A variable already in the contract unit passes through unchanged; a variable +in a recognized non-contract unit is converted automatically (see below): ============================= ================================== -Variable Unit +Variable Contract unit ============================= ================================== Wind (u-component, eastward) m/s Wind (v-component, northward) m/s Surface pressure Pa -Time seconds since a reference date +Time CF datetime (`` since ``) ============================= ================================== -Units are **required** and never assumed. A variable with no ``units`` -attribute, or in a non-contract unit (e.g. ``hPa``/``mbar`` for pressure -or ``knots`` for wind), is **rejected** with an error rather than being -converted -- GeoClaw does not convert units on the read path. Pre-convert -such a file to ``m/s`` / ``Pa`` first (or, for a file known to be in -contract units but lacking the attributes, pass ``assume_units=True`` to -``MetInspector``). The time coordinate must be an **absolute time in -seconds** -- CF ``units`` of the form ``"seconds since "`` (e.g. -``"seconds since 2020-01-01"``). Fortran reads the raw time values as -integer seconds and does **not** scale by the CF unit, so a file whose time -is in ``hours``/``days`` since a reference (such as a raw ERA5 file) must be -converted to ``"seconds since ..."`` first; a bare numeric/duration time -axis with no reference date is rejected outright (see `Time handling`_ below). +**Unit conversion.** A recognized non-contract unit -- ``hPa``/``mbar`` for +pressure, ``knots`` for wind -- is converted automatically: Python computes a +multiplicative scale factor and Fortran applies it on read (a warning notes +the conversion). An unrecognized or dimensionally incompatible unit raises. + +**Missing units.** When a variable has **no** ``units`` attribute, GeoClaw +falls back to the unit documented by the storm *format* if it is known: for +NWS13/OWI files, pressure is taken as ``mbar`` (and converted to ``Pa``) and +wind as ``m/s``. Otherwise a missing unit raises, unless you pass +``assume_units=True`` to ``MetInspector`` to declare that the variables are +already in contract units. + +**Magnitude sanity check.** After units resolve, GeoClaw runs a bounded +min/max check on wind and pressure. A pressure field that is ~1000x too +small (unmistakably ``hPa``/``mbar`` that was mislabeled ``Pa``) is +auto-corrected with a warning; a pressure, wind, or elevation field that is +otherwise physically implausible raises. Pass ``skip_sanity_check=True`` to +the inspector to bypass this for an exotic-but-valid file. + +**Time axis.** The time coordinate must be an absolute CF datetime axis -- +``units`` of the form ``" since "`` (e.g. ``"seconds since +2020-01-01"`` or ``"hours since 2020-01-01"``). A non-second axis such as a +raw ERA5 ``"hours since ..."`` is scaled to seconds automatically (Python +records a ``time_scale`` that Fortran applies), so it no longer needs +pre-conversion. A bare numeric/duration axis with no reference date is +rejected outright (see `Time handling`_ below). Registering a NetCDF storm file ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -359,26 +384,26 @@ Fortran converts raw time values using: .. code:: - storm_time[i] = raw[i] − raw[0] + nc_time_offset + storm_time[i] = nint((raw[i] − raw[0]) * time_scale) + nint(nc_time_offset) Subtracting ``raw[0]`` converts any absolute encoding to -elapsed-since-first-record; adding ``nc_time_offset`` then anchors that -to ``storm.time_offset``. No unit conversion (hours→seconds, -days→seconds) is performed in Fortran: the descriptor stores -``nc_time_offset`` in seconds and the raw values are read as stored. -``write_data`` passes ``time_reference=storm.time_offset`` to -``MetInspector`` automatically; no additional user action is needed. - -The met time axis must be an absolute time in **seconds** -- CF ``units`` of -the form ``"seconds since "`` (e.g. ``"seconds since 2020-01-01"``), -which xarray decodes to datetimes for the offset computation above. Because -Fortran reads the raw time values as integer seconds and does **not** scale -by the CF unit, a file storing time in ``hours``/``days`` since a reference -(such as a raw ERA5 file) must be converted to ``"seconds since ..."`` first. -A bare numeric/duration axis with no reference date is **rejected** outright. -This differs from dtopo NetCDF, where Python collapses the time axis to -``(t0, dt)`` in seconds and therefore any CF duration or datetime axis is -accepted (see `Seafloor deformation (dtopo) from NetCDF`_). +elapsed-since-first-record; multiplying by ``time_scale`` converts the +file's CF time unit to seconds (``1.0`` for ``"seconds since ..."``, +``3600.0`` for ``"hours since ..."``, etc.); adding ``nc_time_offset`` then +anchors that to ``storm.time_offset``. Both factors are resolved in Python +and written to the descriptor, so Fortran never parses CF time units -- +it just multiplies. ``write_data`` passes +``time_reference=storm.time_offset`` to ``MetInspector`` automatically; no +additional user action is needed. + +The met time axis must be an absolute CF datetime axis -- ``units`` of the +form ``" since "`` -- which xarray decodes to datetimes for the +offset computation above. A non-second axis (e.g. a raw ERA5 ``"hours +since ..."``) is scaled to seconds via ``time_scale`` and no longer needs +pre-conversion. A bare numeric/duration axis with no reference date is +**rejected** outright. This parallels dtopo NetCDF, where Python collapses +the time axis to ``(t0, dt)`` in seconds (see +`Seafloor deformation (dtopo) from NetCDF`_). Longitude convention ^^^^^^^^^^^^^^^^^^^^^^ @@ -400,20 +425,26 @@ Architecture overview The system has a strict Python/Fortran split: **Python** handles: file inspection, CF attribute parsing, coordinate -convention detection, fill value resolution, unit verification -(non-contract units are rejected, not converted), -``nc_time_offset`` computation (elapsed seconds from ``storm.time_offset`` -to the first record), crop bound validation, and descriptor writing. +convention detection, fill value resolution, unit resolution (a recognized +non-contract unit becomes a multiplicative ``scale_factor``, and a +non-second time axis a ``time_scale``; missing/unrecognized units raise), +the magnitude sanity check, ``nc_time_offset`` computation (elapsed seconds +from ``storm.time_offset`` to the first record), crop bound validation, and +descriptor writing. **Fortran** handles: opening the NetCDF file at runtime using information from the descriptor, index arithmetic for coordinate normalization (no data copies), domain subsetting via ``start``/``count`` arguments to ``nf90_get_var``, time-slice reads -for met forcing, and converting raw time values to seconds from -``storm.time_offset`` via ``raw[i] − raw[0] + nc_time_offset``. +for met forcing, multiplying each variable by its ``scale_factor`` on read, +and converting raw time values to seconds from ``storm.time_offset`` via +``nint((raw[i] − raw[0]) * time_scale) + nint(nc_time_offset)``. -Fortran assumes the unit contract from ``GEOCLAW_NETCDF_UNITS`` in -``units.py`` without checking. Python enforces it. +All unit logic lives in Python; Fortran is a dumb multiplier. Python +resolves every unit against the contract in ``GEOCLAW_NETCDF_UNITS`` +(``units.py``) to a single ``scale_factor`` (defaulting to ``1.0``, so +contract-unit files are bit-identical) and Fortran applies it without +checking. Class hierarchy (``netcdf_utils.py``) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -430,12 +461,15 @@ Class hierarchy (``netcdf_utils.py``) TopoInspector(NetCDFInspector) - detect fill values within crop region (hard error) - - verify units against the contract (reject non-contract; no conversion) + - resolve units -> scale_factor (convert recognized; reject missing/unknown) + - magnitude sanity check on resolved elevation - no multi-file coverage logic (Fortran handles compositing) MetInspector(NetCDFInspector) - check wind_u, wind_v, pressure present and on same grid/time axis - - verify units against the contract (reject non-contract; no conversion) + - resolve units -> per-variable scale_factor (+ format-unit fallback) + - resolve non-second time axis -> time_scale + - magnitude sanity check on resolved wind/pressure - decode CF datetime axis to seconds from offset (reject bare-numeric) - detect ensemble/member dimensions (hard error if non-singleton) - partial domain coverage is allowed (Fortran fills edges) @@ -464,11 +498,14 @@ Topo (lines in topo.data after topo_type) lon_wrap_offset = 0.0 y_increasing = True dim_order = lat,lon + scale_factor = 1.0 fill_value = -9999.0 fill_action = abort crop_bounds = -100.0 -80.0 20.0 35.0 ``fill_action = abort`` is the only supported value for topography. +``scale_factor`` is the multiplier Fortran applies to the elevation on read +(``1.0`` for a meters file; e.g. ``1000.0`` for a ``km`` file). ``crop_bounds`` is omitted if no crop is specified. dtopo (lines in dtopo.data after the per-file block) @@ -481,6 +518,7 @@ dtopo (lines in dtopo.data after the per-file block) y_name = latitude time_name = time lon_wrap_offset = 0.0 + scale_factor = 1.0 y_increasing = True dim_order = time,lat,lon t0 = 0.0 @@ -489,7 +527,8 @@ dtopo (lines in dtopo.data after the per-file block) ``t0`` and ``dt`` give the first time and the uniform time step in simulation **seconds**; Python collapses the CF time axis to these (scaling by the file's CF ``units``), so Fortran never parses CF time for dtopo. ``time_name`` is the -name of the time coordinate in the file. +name of the time coordinate in the file. ``scale_factor`` is the multiplier +Fortran applies to the deformation on read (``1.0`` for a meters file). Met forcing (body of \*.storm after format header) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -506,14 +545,19 @@ Met forcing (body of \*.storm after format header) fill_value = -9999.0 fill_action = warn time_offset = 0.0 + time_scale = 1.0 / - &variable_info var_name=u10 geoclaw_role=wind_u / - &variable_info var_name=v10 geoclaw_role=wind_v / - &variable_info var_name=msl geoclaw_role=pressure / + &variable_info var_name=u10 geoclaw_role=wind_u scale_factor=1.0 / + &variable_info var_name=v10 geoclaw_role=wind_v scale_factor=1.0 / + &variable_info var_name=msl geoclaw_role=pressure scale_factor=1.0 / (``lon_wrap`` is written only for a geographic longitude axis; ``fill_value`` only when the file declares one; ``crop_bounds`` is carried by the top-level -``crop_extent`` line rather than the descriptor.) +``crop_extent`` line rather than the descriptor.) ``time_scale`` is the +seconds-per-file-time-unit that Fortran multiplies the elapsed time by +(``1.0`` for ``"seconds since ..."``, ``3600.0`` for ``"hours since ..."``); +each ``&variable_info`` carries the per-variable ``scale_factor`` Fortran +applies to that field (e.g. ``100.0`` for an ``hPa``/``mbar`` pressure). The ``&variable_info`` blocks use a manually parsed key=value format (not true repeated Fortran namelists) for compiler portability. Fortran @@ -535,11 +579,14 @@ Defined in ``units.py``: "time": "s", } -Python verifies units against this contract and **rejects** files whose -variables are in non-contract units or lack a ``units`` attribute; no unit -conversion is performed on the read path. Fortran trusts the descriptor -and never checks units. If you add a new ``geoclaw_role``, add its contract -unit here first. +Python resolves each variable's ``units`` attribute against this contract: +a matching unit passes through (``scale_factor = 1.0``), a recognized +non-contract unit yields the multiplicative ``scale_factor`` written to the +descriptor, and a missing or unrecognized unit raises (met additionally +falls back to the storm format's documented unit when the attribute is +absent). Fortran trusts the descriptor, multiplies by ``scale_factor``, and +never checks units. If you add a new ``geoclaw_role``, add its contract unit +here first. Fortran coordinate normalization ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/doc/netcdf_utils_module.rst b/doc/netcdf_utils_module.rst index 359d036..453720e 100644 --- a/doc/netcdf_utils_module.rst +++ b/doc/netcdf_utils_module.rst @@ -11,9 +11,10 @@ netcdf_utils module for NetCDF input The ``netcdf_utils`` module provides the Python layer for reading NetCDF topography and meteorological forcing files in GeoClaw. It handles CF -attribute parsing, coordinate convention detection, unit verification -(non-contract units are rejected, not converted), and descriptor writing so -that Fortran only needs to open the file and read pre-validated indices. +attribute parsing, coordinate convention detection, unit resolution +(recognized non-contract units are converted via a scale factor; missing or +unrecognized units raise), and descriptor writing so that Fortran only needs +to open the file and read pre-validated indices. The main classes are: @@ -23,7 +24,9 @@ The main classes are: subclass; fill-value checking within the crop region, unit enforcement. - :class:`~clawpack.geoclaw.netcdf_utils.MetInspector` — meteorological subclass; wind/pressure variable discovery, CF datetime decoding, unit - verification (non-contract or missing units are rejected). + resolution (recognized non-contract units converted, with a storm-format + fallback for missing units; unrecognized units raise), and a magnitude + sanity check. - :class:`~clawpack.geoclaw.netcdf_utils.CFNormalizer` — adds or repairs CF attributes in place without modifying data values. - :class:`~clawpack.geoclaw.netcdf_utils.DescriptorWriter` — writes the From fb8dec6869cfc2c603a048667b4bc44b4fff8402 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Sun, 12 Jul 2026 18:17:55 -0400 Subject: [PATCH 16/16] Add docs on netcdf compression support --- doc/netcdf.rst | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/doc/netcdf.rst b/doc/netcdf.rst index dcc4947..bd517f2 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -157,12 +157,22 @@ mechanism:: topo.write('bathy.nc', topo_type=4) # elevation float32 topo.write('bathy.nc', topo_type=4, z_dtype='float64') # full precision + topo.write('bathy.nc', topo_type=4, compression=True) # zlib-compressed The elevation variable is written with ``units = "m"``. It is stored on disk as ``float32`` by default -- sub-millimeter precision for Earth topography (``abs(Z) < 10000`` m) at half the file size -- and ``z_dtype='float64'`` selects full double precision. +Pass ``compression=True`` to zlib-compress the elevation variable (zlib level +1 with the byte ``shuffle`` filter). The compressed file is a normal, +randomly-readable NetCDF -- the netCDF library decompresses on read, so no +reader or Fortran change is needed -- and is typically much smaller, +especially for sparse or smooth fields. ``compression`` also accepts an +integer zlib level (``1``--``9``) or a dict of encoding options for full +control; the default ``None`` writes an uncompressed file, bit-identical to +before. + -------------- Seafloor deformation (dtopo) from NetCDF @@ -177,10 +187,14 @@ written to CF-compliant NetCDF using ``dtopo_type=4`` (see :ref:`dtopo`): dtopo = dtopotools.DTopography('deformation.nc', dtopo_type=4) # read dtopo.write('out.nc', dtopo_type=4) # write dtopo.write('out.nc', dtopo_type=4, dz_dtype='float64') # full precision + dtopo.write('out.nc', dtopo_type=4, compression=True) # zlib-compressed As with topography, the deformation contract unit is meters and it is stored ``float32`` by default (pass ``dz_dtype='float64'`` for full double -precision). A recognized non-meter unit (e.g. ``km``) is converted +precision). ``compression=True`` zlib-compresses the deformation variable +(and chunks it one time slice at a time, matching how Fortran reads it), +which shrinks a typical sparse dtopo file several-fold while remaining +directly readable; an integer level or a dict is also accepted. A recognized non-meter unit (e.g. ``km``) is converted automatically on read, exactly as for topography; a missing or unrecognized unit raises. (ASCII dtopo types remain meters-implied, unchanged.)