Warning

Many changes are being implemented in the way topo and dtopo files are handled, in both the Python tools and the Fortran code. See Changes to topo and dtopo handling (planned for v5.15.0) for a summary.

GeoClaw NetCDF Input System

This document covers the NetCDF input pipeline introduced in the refactor-netcdf-support [PR #701](https://github.com/clawpack/geoclaw/pull/701) (merged but not yet released). It has two sections: a user guide for scientists who want to use NetCDF files as input, and a developer reference for those working on the implementation.

---

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:

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:

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:

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:

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 dtopo):

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 <time_reference>":

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 "<unit> since <date>" 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 (<unit> since <date>)

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 "<unit> since <date>" (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:

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:

import numpy as np
from clawpack.geoclaw.met.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:

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:

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 MetInspector explicitly and pass it via met_inspector. When a pre-built inspector is supplied, auto-discovery and var_mapping validation are bypassed entirely:

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:

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:

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 "<unit> since <date>" – 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:

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