diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..26f2e21 --- /dev/null +++ b/.gitignore @@ -0,0 +1,66 @@ +# ----------------------------------------------------------------------------- +# 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/ +doc/pyclaw.log diff --git a/doc/changes_to_master.rst b/doc/changes_to_master.rst index 7443453..a63f05e 100644 --- a/doc/changes_to_master.rst +++ b/doc/changes_to_master.rst @@ -21,6 +21,26 @@ the menu on the left hand side of this page. Changes that are not backward compatible ---------------------------------------- +- **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 --------------- @@ -64,6 +84,69 @@ 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. Files are written coarsest-first (finest last), matching + the traditional GeoClaw listing order; the last file listed in ``topo.data`` + is assigned rank 1 (highest priority) by Fortran, with no Fortran-side + sorting. ``rundata.topo_data.override_order = True`` preserves + user-specified list order; when used, the finest (highest-resolution) file + should be listed last. See :ref:`topo_order`. + +- **topo_type=1 deprecated.** + Reading and writing ``topo_type=1`` (``x y z`` one-point-per-line ASCII) + now emit a ``DeprecationWarning``. Setting any preprocessing attribute + before reading a type-1 file raises ``NotImplementedError``. To convert:: + + 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. + +- **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/geoclaw.rst b/doc/geoclaw.rst index c3ee6d0..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 @@ -69,6 +70,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 fdcb9af..bd517f2 100644 --- a/doc/netcdf.rst +++ b/doc/netcdf.rst @@ -1,12 +1,694 @@ -:orphan: +.. _netcdf_input: -.. _netcdf: +GeoClaw NetCDF Input System +=========================== -.. seealso:: :ref:`topo_netcdf`. +This document covers the NetCDF input pipeline introduced in the +``refactor-netcdf-support`` PR (first released after v5.14.0). It has two sections: a user +guide for scientists who want to use NetCDF files as input, and a developer +reference for those working on the implementation. -========================== -Using NetCDF output -========================== +-------------- -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 +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- 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. +- 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. +- 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 (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 + supported for topography. + +Registering a NetCDF topo file in setrun.py +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +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 you want to specify crop bounds, or if GeoClaw cannot find the elevation +variable automatically, use ``topo_entries()`` to inspect the file: + +.. code:: python + + 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, +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: + +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 = 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. + +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 +``TopoInspector`` and use ``topo_entries()`` to register the file: + +.. code:: python + + 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. + +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. + +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 + 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 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +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 + 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). ``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.) + +**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 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +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 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 Contract unit +============================= ================================== +Wind (u-component, eastward) m/s +Wind (v-component, northward) m/s +Surface pressure Pa +Time CF datetime (`` since ``) +============================= ================================== + +**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 +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +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``. 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 + + 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/era5_forcing.nc'] + storm.write('isaac.storm', file_format='data') + +**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', + var_mapping={'pressure': 'prmsl'}) + +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'}) + +**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.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 MetInspector + + 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:: + + ``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 + ``MetInspector``'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] = nint((raw[i] − raw[0]) * time_scale) + nint(nc_time_offset) + +Subtracting ``raw[0]`` converts any absolute encoding to +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 +^^^^^^^^^^^^^^^^^^^^^^ + +Longitude convention ([0, 360] vs [-180, 180]) is detected and +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. + +-------------- + +Developer Reference +------------------- + +Architecture overview +~~~~~~~~~~~~~~~~~~~~~ + +The system has a strict Python/Fortran split: + +**Python** handles: file inspection, CF attribute parsing, coordinate +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, 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)``. + +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``) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +:: + + NetCDFInspector + - 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 + + TopoInspector(NetCDFInspector) + - detect fill values within crop region (hard error) + - 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 + - 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) + + 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 + x_name = longitude + y_name = latitude + 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) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +:: + + var_name = dz + x_name = longitude + 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 + 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. ``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) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +:: + + &file_info + x_name = longitude + y_name = latitude + time_name = valid_time + dim_order = time,lat,lon + lon_wrap = 360 + y_increasing = True + fill_value = -9999.0 + fill_action = warn + time_offset = 0.0 + time_scale = 1.0 + / + &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.) ``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 +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", + } + +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 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Fortran applies coordinate normalization via index arithmetic only -- no +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. + +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`` (in ``util.py``) remains a parallel + implementation of variable name discovery alongside + ``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 ``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 + 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 ``MetInspector`` for the new variable +3. ``DescriptorWriter`` requires no change -- new ``&variable_info`` + 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_inspector.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_inspector.py + test_topo_inspector.py + test_met_inspector.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 +---------- + +- **Consolidate util.get_netcdf_names**: ``Storm.write_data`` now + calls ``MetInspector`` and ``DescriptorWriter`` for NetCDF met + forcing; the remaining step is to make ``util.get_netcdf_names`` + 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 + ``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 ``NetCDFInspector`` diff --git a/doc/netcdf_utils_module.rst b/doc/netcdf_utils_module.rst new file mode 100644 index 0000000..453720e --- /dev/null +++ b/doc/netcdf_utils_module.rst @@ -0,0 +1,42 @@ + +.. _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 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: + +- :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 datetime decoding, unit + 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 + 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/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. - - diff --git a/doc/setrun_geoclaw.rst b/doc/setrun_geoclaw.rst index e67bc2b..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 +.. attribute:: rundata.topo_data.topofiles : list - *topofiles* should be a list of the form *[file1info, file2info, etc.]* - where each element is itself a list of the form + A list of topography file entries. The preferred form uses + :class:`~clawpack.geoclaw.topotools.Topography` objects:: - [topotype, fname] + from clawpack.geoclaw.topotools import Topography - with values + t = Topography() + t.path = 'etopo1.tt2' + t.topo_type = 2 + rundata.topo_data.topofiles.append(t) - *topotype* : integer + .. deprecated:: + The ``[topotype, fname]`` list format is deprecated. It still works + but emits a ``DeprecationWarning``. Use ``Topography`` objects as + shown above. - 1,2 or 3 depending on the format of the file (see :ref:`topo`). + **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`. - *fname* : string + See :ref:`topo` for information about the topo file formats + (``topo_type`` values 2, 3, 4). - the name of the topo file. + See :ref:`topo_order` for how priority is determined when topo files + overlap, and :ref:`topodata_format` for the ``topo.data`` file format. - **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`. +.. _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 @@ -453,7 +515,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. 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 diff --git a/doc/topo_order.rst b/doc/topo_order.rst index 6fe2020..311394f 100644 --- a/doc/topo_order.rst +++ b/doc/topo_order.rst @@ -21,6 +21,13 @@ 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. 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 recent and hence a better dataset to use @@ -31,8 +38,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 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). .. warning :: This relative tolerance was introduced in v5.14.0, @@ -48,12 +55,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 +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. + +``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..e849578 --- /dev/null +++ b/doc/topodata_format.rst @@ -0,0 +1,305 @@ +.. _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 (10 lines, all types) +------------------------------------ + +One block follows for each topography file, in priority order (coarsest +resolution first, finest last; see :ref:`priority_convention` below). + +:: + + '' # topo_path + # topo_type + # crop_extent [x1 x2 y1 y2] + # coarsen + # buffer + # align [x y] + # x_shift + # y_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). + 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). + +``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 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 + x_name=lon + y_name=lat + y_increasing=True + dim_order=y,x + lon_wrap_offset=0.0 + fill_action=abort + +``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_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. + + +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 # y_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 # y_shift + 0 # z_shift + T # negate_z + + var_name=elevation + x_name=lon + y_name=lat + y_increasing=True + dim_order=y,x + lon_wrap_offset=0.0 + fill_action=abort + + +.. _priority_convention: + +Priority Convention +------------------- + +The **last** file block in ``topo.data`` is assigned **rank 1** by Fortran +(highest priority in overlap resolution). Fortran stores this mapping in +``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 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 last 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, 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 bf6ffb7..1bff5a1 100644 --- a/doc/topotools.rst +++ b/doc/topotools.rst @@ -4,6 +4,145 @@ Python tools for working with topo and dtopo -------------------------------------------- +.. seealso:: + - :ref:`topo` + - :ref:`topo_order` + - :ref:`topodata_format` + + +Preprocessing attributes +~~~~~~~~~~~~~~~~~~~~~~~~~ + +:class:`~clawpack.geoclaw.topotools.Topography` objects support seven +preprocessing attributes that are applied automatically by +:meth:`~clawpack.geoclaw.topotools.Topography.read` in this order: + +1. ``negate_z`` — flip sign of Z (independent of ``topo_type < 0``). +2. ``z_shift`` — add a constant to all non-missing Z values. +3. ``x_shift`` — add a constant to all x coordinates. +4. ``crop_extent``, ``buffer``, ``align``, ``coarsen`` — crop and subsample + via :meth:`~clawpack.geoclaw.topotools.Topography.crop`. + +Set them before calling ``read()``:: + + from clawpack.geoclaw.topotools import Topography + + t = Topography() + t.crop_extent = [-100., -60., 10., 50.] + t.coarsen = 2 + t.z_shift = 10.0 + t.read('bathymetry.nc', topo_type=4) + +See :ref:`setrun_topo_preprocessing` for a full attribute table with types +and defaults, and non-obvious behavior notes. + + +Lazy-load pattern for NetCDF (read_header) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For NetCDF files (``topo_type=4``), +:meth:`~clawpack.geoclaw.topotools.Topography.read_header` uses +:class:`~clawpack.geoclaw.netcdf_utils.TopoInspector` to detect coordinate +variable names via CF conventions (``standard_name``, ``axis``, or common +fallback names such as ``lon``/``lat``, ``x``/``y``). After calling +``read_header()``, ``extent`` and ``delta`` are populated without loading the +elevation array. Accessing ``Z`` triggers a deferred ``read()``: + +.. code-block:: python + + t = Topography() + t.path = 'large_dem.nc' + t.topo_type = 4 + t.read_header() # fast: reads coordinate arrays only + print(t.extent) # available immediately + print(t.delta) # available immediately + z = t.Z # triggers full read on first access + +This pattern is also used internally by ``TopographyData._compute_priority_order`` +to determine file resolution without loading elevation data. + + +Vertical datum metadata +~~~~~~~~~~~~~~~~~~~~~~~~~ + +:class:`~clawpack.geoclaw.topotools.Topography` and +:class:`~clawpack.geoclaw.dtopotools.DTopography` carry an optional ``datum`` +attribute recording the vertical reference of the elevation/deformation data +(e.g. ``'MSL'`` or ``'NAVD88'``). It defaults to ``None``, is populated from a +NetCDF ``vertical_datum`` (or ``datum``) attribute on read, and is written +back out on NetCDF write (``topo_type=4`` / ``dtopo_type=4``). ASCII formats +have no datum field, so the value is not persisted for types 1/2/3. + +The datum is **informational only** -- GeoClaw performs no vertical-datum +transformation and always uses the Z values as given. As a guard, both +:meth:`~clawpack.geoclaw.data.TopographyData.write` (when producing +``topo.data``) and :meth:`~clawpack.geoclaw.data.DTopoData.write` (``dtopo.data``) +issue a warning if the files they are given carry more than one distinct +datum, since mixing vertical references without converting between them is a +likely error. + + +.. deprecated:: + ``topo_type=1`` (three-column ``x y z`` ASCII, one point per line) is + deprecated. Reading emits a ``DeprecationWarning``; writing also emits a + ``DeprecationWarning``; setting any preprocessing attribute before reading + raises ``NotImplementedError``. + + To convert a type-1 file:: + + t = Topography() + t.read('old.tt1', topo_type=1) # DeprecationWarning + t.write('new.tt2', topo_type=2) # save as type 2 + + Genuinely unstructured (scattered) point data cannot be converted this + way. Either grid it externally (e.g. ``scipy.interpolate``, GMT) or use + :meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` (see + `Gridding unstructured (scattered) data`_ below). + + +Gridding unstructured (scattered) data +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +GeoClaw's solvers require topography on a regular, logically rectangular grid +(Cartesian x-y in meters, or lon-lat -- see ``coordinate_system`` in +:ref:`setrun_geoclaw`), but survey or sounding data often comes as scattered +``(x, y, z)`` points. A +:class:`~clawpack.geoclaw.topotools.Topography` constructed with +``unstructured=True`` holds such points in its ``x``, ``y``, and ``z`` +arrays, and +:meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` +interpolates them onto a regular grid, optionally filling gaps from one or +more coarser (structured or unstructured) "fill" topographies: + +.. code-block:: python + + from clawpack.geoclaw.topotools import Topography + + # Scattered survey points (x, y, z) + survey = Topography(unstructured=True) + survey.x = x_points + survey.y = y_points + survey.z = z_values + + # A coarser regional DEM used to fill gaps between survey points + regional = Topography(path='regional.tt3', topo_type=3) + regional.read() + + survey.interp_unstructured(regional, extent=[x1, x2, y1, y2], + proximity_radius=100.0) + survey.unstructured # now False -- it holds a regular grid + survey.write('combined.tt3', topo_type=3) + +The grid spacing is taken from the minimum spacing of the scattered points +(bounded below by ``delta_limit`` meters) unless set explicitly with +``delta``. A fill point is used only where it lies inside ``extent``, carries +valid data, and is more than ``proximity_radius`` meters from every scattered +point, so the higher-resolution survey data is preferred wherever it exists. +Missing fill values (``NaN`` in memory, or a numeric ``no_data_value``) are +dropped. See +:meth:`~clawpack.geoclaw.topotools.Topography.interp_unstructured` for the +full set of parameters. + .. toctree:: :maxdepth: 1