From c893f238b38b32771e49507824df3c1e76d9a592 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Fri, 4 Sep 2026 18:13:02 -0400 Subject: [PATCH 1/7] Make crop_extent and antimeridian failures loud instead of silent Eleven silent failures in crop_extent handling produced wrong answers with no message: wrapped spellings gave empty grids, sub-cell and non-overlapping crops gave the full file, a cropped type-4 read reported the full-file extent, and buffer was dropped for every descriptor-cropped NetCDF file. Adds the first Fortran regression coverage for the descriptor-crop path and documents what a Topography represents across the antimeridian. Signed-off-by: Kyle Mandli Assisted-by: claude claude-opus-5[1m] --- src/2d/shallow/topo_module.f90 | 5 + src/python/geoclaw/data.py | 45 ++- src/python/geoclaw/topotools.py | 170 ++++++++++- tests/regression/topo_crop/Makefile | 51 ++++ tests/regression/topo_crop/setrun.py | 74 +++++ .../topo_crop/test_topo_descriptor_crop.py | 282 ++++++++++++++++++ tests/test_topotools_preprocessing.py | 229 +++++++++++++- 7 files changed, 846 insertions(+), 10 deletions(-) create mode 100644 tests/regression/topo_crop/Makefile create mode 100644 tests/regression/topo_crop/setrun.py create mode 100644 tests/regression/topo_crop/test_topo_descriptor_crop.py diff --git a/src/2d/shallow/topo_module.f90 b/src/2d/shallow/topo_module.f90 index 7a85b0be6..63e2f30e4 100644 --- a/src/2d/shallow/topo_module.f90 +++ b/src/2d/shallow/topo_module.f90 @@ -1481,6 +1481,11 @@ subroutine read_topo_header(fname,topo_type,mx,my,xll,yll,xhi,yhi,dx,dy,topo_idx (xlocs <= nc_crop_bounds(2, topo_idx)) y_in_dom = (ylocs >= nc_crop_bounds(3, topo_idx)) .and. & (ylocs <= nc_crop_bounds(4, topo_idx)) + ! buffer applies to a descriptor crop exactly as it does to + ! topo_crop_extent below; leaving nbuf4 = 0 here silently + ! dropped topo_buffer for every file written by + ! TopoInspector.topo_entries(), which always sets crop_bounds. + nbuf4 = topo_buffer(topo_idx) else if (present(topo_idx) .and. & any(topo_crop_extent(:,topo_idx) /= 0.0d0)) then ! topo_crop_extent is in domain coordinates (after nc_lon_wrap_offset diff --git a/src/python/geoclaw/data.py b/src/python/geoclaw/data.py index 47300c3ed..488b31e21 100755 --- a/src/python/geoclaw/data.py +++ b/src/python/geoclaw/data.py @@ -378,8 +378,51 @@ def write(self, data_source='setrun.py', out_file='topo.data'): os.path.join(os.path.dirname(out_file), topo.path)) f = self._out_file + + # topo_type may still be None when a Topography was built by + # hand and never read; the :3d format below would raise a + # TypeError naming neither the file nor the attribute. + _topo_type = topo.topo_type + if _topo_type is None: + from clawpack.geoclaw import topotools + _topo_type = topotools.determine_topo_type( + topo.path, default=None) + if _topo_type is None: + raise ValueError( + f"topo_type is not set for {topo.path} and cannot " + f"be inferred from its extension. Set topo.topo_type " + f"explicitly, or pass topo_type= to Topography().") + topo.topo_type = _topo_type + + # Type 5 (GeoTIFF) is readable in Python but has no case(5) in + # the Fortran reader, which aborts with "Unrecognized topo_type". + if abs(_topo_type) == 5: + warnings.warn( + f"{topo.path} is written to {out_file} as topo_type=5 " + f"(GeoTIFF). GeoClaw's Fortran reader does not support " + f"type 5 and will abort; convert to topo_type 3 or 4 " + f"(Topography.write) before running.", + UserWarning, stacklevel=2) + + # A crop that crosses the antimeridian cannot be expressed as a + # single descriptor: crop_bounds would read "170.0 -170.0", + # which Fortran resolves to mx=0, my=0 -- an empty topo file, + # silently. TopoInspector.topo_entries() splits such a crop + # into two entries with the appropriate lon_wrap_offset. + if (topo.crop_extent is not None + and topo.crop_extent[0] >= topo.crop_extent[1] + and getattr(topo, '_netcdf_meta', None) is None): + raise ValueError( + f"crop_extent {list(topo.crop_extent)} for {topo.path} " + f"is descending in longitude, i.e. it crosses the " + f"antimeridian. Written directly it would produce an " + f"empty grid in Fortran with no error. Build the entries " + f"with TopoInspector.topo_entries(), which emits one " + f"entry per side of the seam with the correct " + f"lon_wrap_offset, and append those to topofiles.") + f.write(f"\n'{fname}' # topo_path\n") - f.write(f"{topo.topo_type:3d} # topo_type\n") + f.write(f"{_topo_type:3d} # topo_type\n") _write_preprocessing_block(f, topo) # For NetCDF (type 4): write the CF descriptor block that diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py index a6343f4ba..eabfe8808 100644 --- a/src/python/geoclaw/topotools.py +++ b/src/python/geoclaw/topotools.py @@ -29,6 +29,7 @@ """ import os +import warnings import numpy @@ -116,7 +117,34 @@ def _crop_indices(x, y, crop_extent, coarsen, buffer, align): ``None`` means no phase snap (start at the crop window). ``coarsen`` and ``buffer`` are assumed already ``int()``-coerced. + + :Raises: + - *ValueError* if *crop_extent* is not increasing in both coordinates. A + descending longitude pair is the natural way to *spell* a crop across the + antimeridian, and treating it as an ordinary window silently produced an + empty grid, so it is rejected explicitly. + - *ValueError* if *crop_extent* overlaps the data but contains no grid + point (a window narrower than one cell, falling between two points). + + Warns when *crop_extent* extends beyond the data and is therefore clipped: + a crop is never wrapped, only reduced. """ + # A descending pair is not an empty window; it is almost always an attempt + # to cross the antimeridian. Both index lookups below succeed when + # crop_extent[0] > crop_extent[1], giving iupper < ilower and so a zero-size + # slice -- an empty Topography whose `.extent` then raises something opaque + # far from the cause. Fail here instead. + if crop_extent[0] >= crop_extent[1] or crop_extent[2] >= crop_extent[3]: + raise ValueError( + f"crop_extent must increase in both coordinates, got " + f"{list(crop_extent)}. To cross the antimeridian use the " + f"continuous spelling (e.g. [-211, -99]) rather than the wrapped " + f"one ([170, -170]) -- and note that a Topography does not wrap on " + f"its own: longitude wrapping is applied by the Fortran reader " + f"from a descriptor written by TopoInspector.topo_entries(), which " + f"emits two entries for a cross-seam crop. See the 'Region " + f"terminology' note in the Topography docstring.") + # dx/dy computed the same way as the `delta` property (round to 15 places), # so the align fractional-offset search matches crop() bit-for-bit. dx = numpy.round(abs(x[1] - x[0]), 15) @@ -133,6 +161,38 @@ def _crop_indices(x, y, crop_extent, coarsen, buffer, align): except IndexError: # crop_extent does not overlap the data return None + if iupper < ilower or jupper < jlower: + # The window overlaps the extent but falls strictly between two grid + # points, so it contains no data. Unlike a genuine non-overlap (which + # has a documented full-file fallback) nothing pins this case, and + # returning the full file for a crop the user narrowed *too far* is + # never what was meant. + raise ValueError( + f"crop_extent {list(crop_extent)} lies between grid points and " + f"contains no data: the grid spacing is dx={dx}, dy={dy}. Widen " + f"the crop to at least one cell, or use buffer= to include the " + f"surrounding points.") + + # Silent clipping is the other half of the antimeridian confusion: a crop + # written in continuous coordinates ([-211, -99] for a file on [-180, 180]) + # is not wrapped, it is quietly reduced to the part that exists. Say so + # when more than a cell is dropped; the result is still returned, because + # over-wide crops are a legitimate and common way to say "all of this". + _clipped = [] + if crop_extent[0] < x[0] - dx: + _clipped.append(f"x1 {crop_extent[0]} -> {x[0]}") + if crop_extent[1] > x[-1] + dx: + _clipped.append(f"x2 {crop_extent[1]} -> {x[-1]}") + if crop_extent[2] < y[0] - dy: + _clipped.append(f"y1 {crop_extent[2]} -> {y[0]}") + if crop_extent[3] > y[-1] + dy: + _clipped.append(f"y2 {crop_extent[3]} -> {y[-1]}") + if _clipped: + warnings.warn( + f"crop_extent {list(crop_extent)} extends past the data, which " + f"covers x=[{x[0]}, {x[-1]}], y=[{y[0]}, {y[-1]}]; it was clipped " + f"({', '.join(_clipped)}). A Topography does not wrap: to cross " + f"the antimeridian use TopoInspector.topo_entries().") # Shift indices if needed for alignment (matches crop() lines historically # at 2085-2099: pick the low index whose coord best lands on `align`). @@ -396,7 +456,6 @@ def _resolve_crop_extent(crop_extent, deprecated): value as ``crop_extent``; raise ``TypeError`` if ``crop_extent`` is also supplied (ambiguous). """ - import warnings for name, value in deprecated.items(): if value is _CROP_EXTENT_UNSET: continue @@ -466,6 +525,52 @@ class Topography(object): Convention: the ``_extent`` suffix denotes domain coordinates; ``_bounds`` denotes file coordinates. + :The antimeridian, and what a cropped Topography represents: + + **A Topography never wraps.** It holds one ascending ``x`` array and one + ascending ``y`` array, so it can represent a rectangle in a continuous + coordinate frame and nothing else. Wrapping is not a property of the + object; it is applied by the *Fortran* reader, from a ``lon_wrap_offset`` + that only exists in a NetCDF descriptor. This is the distinction behind + most antimeridian confusion, so it is worth being concrete about the three + cases: + + 1. **Wrapped spelling, e.g. ``crop_extent=[170, -170, ...]``.** Rejected + with a ``ValueError``. There is no ascending window it could mean, and + silently producing an empty grid (which is what a descending pair used + to do) is worse than refusing. + + 2. **Continuous spelling, e.g. ``crop_extent=[-211, -99, ...]`` for a file + on ``[-180, 180]``.** Accepted, but the part that lies off the file is + *clipped, not wrapped* -- you get ``[-180, -99]`` and a + ``UserWarning`` saying so. It is a legitimate way to say "as far west + as this file goes"; it is not a way to cross the seam. + + 3. **Genuinely crossing the seam.** Not expressible on a single + ``Topography``, and not expressible in a single ``topo.data`` entry + either. Use :meth:`netcdf_utils.TopoInspector.topo_entries`, which + returns one entry per side of the cut with the appropriate + ``lon_wrap_offset`` and file-coordinate ``crop_bounds``; append those to + ``rundata.topo_data.topofiles``. Writing a descending ``crop_extent`` + to ``topo.data`` directly raises, because Fortran would resolve + ``crop_bounds = 170.0 -170.0`` to ``mx=0, my=0``: an empty topography, + with no error. + + The general shape of it: **Python is the single-rectangle case; the wrap + lives in the Fortran interface.** The same split explains why + ``coordinate_system`` gates wrapping (a projected x axis in meters has no + seam to cross) and why ``crop_bounds`` is in file coordinates while + ``crop_extent`` is in domain coordinates. + + :Order of preprocessing operations: + + ``crop_extent`` -> ``align`` -> ``buffer`` -> ``coarsen``, applied in that + order by both :meth:`crop` and the Fortran reader + (``apply_align_buffer_coarsen``). Because the strided subsample is last, + **``buffer`` counts coarsened output points, not native file points**: the + index window is widened by ``buffer * coarsen`` native points on each side, + so ``buffer=2, coarsen=4`` adds 2 points to each edge of the result, not 8. + """ @property @@ -853,7 +958,10 @@ class docstring). Passing it here is equivalent to setting the ``align=[integer_lon, integer_lat]`` to lock the coarsened grid to a fixed lattice regardless of the requested ``crop_extent``. - *buffer* (int) - grid points to keep outside ``crop_extent`` on each - side; see :meth:`crop`. + side. These are *output* points: with ``coarsen > 1`` the window + grows by ``buffer * coarsen`` native points, so the result gains + ``buffer`` points per edge, not ``buffer * coarsen``. See + :meth:`crop`. - *stride* (list or int) - **Deprecated**: use ``coarsen`` instead. A NetCDF-only knob that silently did nothing for ASCII reads and used a different alignment convention. A scalar (or equal-valued list) is @@ -897,7 +1005,6 @@ class docstring). Passing it here is equivalent to setting the # for ASCII reads and used a different alignment convention than # crop()/coarsen. Map it onto the unified scalar `coarsen`. if stride is not None: - import warnings warnings.warn( "The 'stride' argument to Topography.read() is deprecated; use " "'coarsen' (a scalar subsampling factor) instead. 'coarsen' is " @@ -951,6 +1058,24 @@ class docstring). Passing it here is equivalent to setting the raise ValueError("topo_type must be specified") if self.unstructured: + # crop() already refuses unstructured data (NotImplementedError); + # the read path used to attempt its own filter and then index a + # Python list as if it were an array, dying with an unrelated + # TypeError far from the cause. Refuse consistently and early. + _unstructured_preprocessing = ( + self.crop_extent is not None + or self.coarsen != 1 + or self.buffer != 0 + or self.align is not None + ) + if _unstructured_preprocessing: + raise NotImplementedError( + "Preprocessing attributes (crop_extent, coarsen, buffer, " + "align) are not supported for unstructured data; " + "Topography.crop() refuses them for the same reason. Grid " + "the data first (e.g. interp_unstructured), then crop the " + "result.") + # Read in the data as series of tuples data = numpy.loadtxt(self.path) points = [] @@ -981,7 +1106,6 @@ class docstring). Passing it here is equivalent to setting the else: # Data is in one of the GeoClaw supported formats if abs(self.topo_type) == 1: - import warnings warnings.warn( "topo_type=1 is deprecated. Convert to topo_type=2, 3, or 4:\n" " topo.read() # load the type-1 file\n" @@ -1144,6 +1268,11 @@ class docstring). Passing it here is equivalent to setting the # crop_extent misses the file: fall back to the full grid # at native resolution (mirrors crop() returning None -> # no-op), rather than coarsening the whole file. + warnings.warn( + f"crop_extent {list(_ce)} did not overlap " + f"{self.path}; the full file was read uncropped at " + f"native resolution. The Fortran reader treats a " + f"non-overlapping crop as fatal.") _il, _iu, _jl, _ju = 0, _nx, 0, _ny _step = 1 else: @@ -1231,6 +1360,13 @@ class docstring). Passing it here is equivalent to setting the # Make sure these are set to None to force re-generating: self._X = None self._Y = None + # _extent and _delta are derived from _x/_y, which were just + # replaced. read_header() populates them from the file header, so + # without this a cropped topo_type=4 read (whose crop is applied + # while reading the hyperslab, bypassing the property setters) + # reported the *full file* extent alongside cropped data. + self._extent = None + self._delta = None # Normalize missing data to NaN in memory. The numeric # self.no_data_value is only the on-file/Fortran sentinel (written @@ -1278,6 +1414,14 @@ class docstring). Passing it here is equivalent to setting the buffer=int(self.buffer), align=self.align, ) + if _cropped is None: + # crop() already warned about the non-overlap; say what the + # consequence is here, because keeping the *full* file is + # surprising and the Fortran reader would instead abort. + warnings.warn( + f"crop_extent {list(self.crop_extent)} did not overlap " + f"{self.path}; the full file was read uncropped. The " + f"Fortran reader treats this as fatal.") if _cropped is not None: self._x = _cropped._x self._y = _cropped._y @@ -1541,7 +1685,6 @@ def write(self, path, topo_type=None, no_data_value=None, fill_value=None, outfile.write("%s %s %s\n" % (self.x[i], self.y[i], topo)) elif topo_type == 1: - import warnings warnings.warn( "Writing topo_type=1 is deprecated. Prefer topo_type=2 or 3 for ASCII " "output, or topo_type=4 for NetCDF. Type-1 output will be removed in " @@ -2149,7 +2292,13 @@ def crop(self, crop_extent=None, coarsen=1, buffer=0, align=None, - *buffer* (int): integer number of grid points to keep on each side of *crop_extent* (when possible) -- NOT a coordinate distance (cf. ``interp_unstructured``'s ``buffer_length``, which is in meters). - Truncated to an integer via ``int()``. + Truncated to an integer via ``int()``. Counted in *coarsened + output* points: the operations apply in the order crop -> align -> + buffer -> coarsen, so the native index window is widened by + ``buffer * coarsen`` and the strided subsample then keeps + ``buffer`` of those per edge. ``buffer=2, coarsen=4`` therefore + adds 2 points to each edge of the result, not 8. Fortran + ``apply_align_buffer_coarsen`` does the same arithmetic. - *align* (tuple): (xalign,yalign) = desired alignment if coarsening Setting *buffer > 0* may be useful to insure that the @@ -2218,7 +2367,13 @@ def crop(self, crop_extent=None, coarsen=1, buffer=0, align=None, # read path so ASCII and NetCDF reads of the same data match exactly). idx = _crop_indices(self.x, self.y, crop_extent, coarsen, buffer, align) if idx is None: - print('*** crop_extent does not overlap topo') + # Warned rather than printed so a caller can catch, filter or + # escalate it; the Fortran reader treats the same condition as + # fatal (topo_module.f90: "does not overlap topo file", stop 1), + # so a run that ignores this here will fail there. + warnings.warn( + f"crop_extent {list(crop_extent)} does not overlap this " + f"topography (extent {list(self.extent)}); no crop applied.") return None ilower, iupper, jlower, jupper = idx @@ -2486,7 +2641,6 @@ def read_netcdf(path, zvar=None, extent='all', coarsen=1, return_topo=True, depending on ``return_topo`` / ``return_xarray`` (unchanged contract). """ - import warnings warnings.warn( "topotools.read_netcdf is deprecated; use " "topotools.fetch_remote_topo (or Topography.read(topo_type=4)) instead.", diff --git a/tests/regression/topo_crop/Makefile b/tests/regression/topo_crop/Makefile new file mode 100644 index 000000000..0c55af778 --- /dev/null +++ b/tests/regression/topo_crop/Makefile @@ -0,0 +1,51 @@ + +# Makefile for the topo descriptor-crop regression case. +# Builds a standard single-layer GeoClaw executable; identical to the +# met_forcing Makefile, since what is under test is the topo reader, not the +# forcing. NetCDF is required (the descriptor path is topo_type=4 only). +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +CLAW_PKG = geoclaw # Clawpack package to use +EXE = xgeoclaw # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +FFLAGS ?= + +# NetCDF support is opt-in so the parametric (Holland) build stays dependency +# free. The gridded netCDF aux case builds with: +# make new USE_NETCDF=1 NETCDF_FFLAGS="$(nf-config --fflags)" \ +# NETCDF_LFLAGS="$(nf-config --flibs)" +# (the test injects these automatically via build_executable make_vars). +ifeq ($(USE_NETCDF),1) +FFLAGS += -DNETCDF $(NETCDF_FFLAGS) +LFLAGS += $(FFLAGS) $(NETCDF_LFLAGS) +endif + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow +include $(GEOLIB)/Makefile.geoclaw + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + +RIEMANN = $(CLAW)/riemann/src + +MODULES = \ + +SOURCES = \ + $(RIEMANN)/rpn2_geoclaw.f \ + $(RIEMANN)/rpt2_geoclaw.f \ + $(RIEMANN)/geoclaw_riemann_utils.f \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +### DO NOT remove this line - make depends on it ### diff --git a/tests/regression/topo_crop/setrun.py b/tests/regression/topo_crop/setrun.py new file mode 100644 index 000000000..83beb3382 --- /dev/null +++ b/tests/regression/topo_crop/setrun.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python +# encoding: utf-8 +"""Run configuration for the topo descriptor-crop regression case. + +The domain is set by the test, not here: each case picks bounds that sit in a +specific relationship to the topo file's crop window (inside it, or inside the +crop-plus-buffer ring), which is the whole point of the test. Everything else +is the smallest configuration that still exercises ``read_topo_file``: one AMR +level, a few cells, one step. +""" + +from clawpack.clawutil import data + + +def setrun(claw_pkg='geoclaw'): + rundata = data.ClawRunData(claw_pkg, 2) + + cd = rundata.clawdata + cd.num_dim = 2 + # Placeholders; the test overwrites lower/upper/num_cells before writing. + cd.lower = [-95.0, 22.0] + cd.upper = [-90.0, 25.0] + cd.num_cells = [20, 12] + cd.num_eqn = 3 + cd.num_aux = 3 + cd.capa_index = 2 + + cd.t0 = 0.0 + cd.output_style = 1 + cd.num_output_times = 1 + cd.tfinal = 1.0 + cd.output_t0 = False + cd.dt_initial = 0.1 + cd.cfl_desired = 0.75 + cd.steps_max = 100 + cd.verbosity = 1 + cd.dt_variable = True + cd.bc_lower = ['extrap', 'extrap'] + cd.bc_upper = ['extrap', 'extrap'] + + # The GeoClaw Riemann solver always returns 3 waves and uses f-waves. + # These are not defaults: num_waves defaults to 1, so leaving them out lets + # the solver write past the end of the wave arrays. That corrupts memory + # and shows up as an intermittent segfault *after* the run finishes -- with + # fort.geo truncated, which looks like a topo problem and is not one. + cd.order = 2 + cd.transverse_waves = 2 + cd.num_waves = 3 + cd.limiter = ['mc', 'mc', 'mc'] + cd.use_fwaves = True + cd.source_split = 'godunov' + + amrdata = rundata.amrdata + amrdata.amr_levels_max = 1 + amrdata.refinement_ratios_x = [2] + amrdata.refinement_ratios_y = [2] + amrdata.refinement_ratios_t = [2] + amrdata.aux_type = ['center', 'capacity', 'yleft'] + + geo = rundata.geo_data + geo.gravity = 9.81 + geo.coordinate_system = 2 + geo.earth_radius = 6367.5e3 + geo.sea_level = 0.0 + geo.dry_tolerance = 1.e-3 + geo.friction_forcing = False + + rundata.refinement_data.wave_tolerance = 1.e-1 + + return rundata + + +if __name__ == '__main__': + setrun().write() diff --git a/tests/regression/topo_crop/test_topo_descriptor_crop.py b/tests/regression/topo_crop/test_topo_descriptor_crop.py new file mode 100644 index 000000000..6170ba63c --- /dev/null +++ b/tests/regression/topo_crop/test_topo_descriptor_crop.py @@ -0,0 +1,282 @@ +#!/usr/bin/env python +# encoding: utf-8 +"""End-to-end regression for the topo descriptor-crop path (``topo_type=4``). + +This is the first test that compiles GeoClaw and feeds it a NetCDF topo file +through a descriptor written by ``TopoInspector.topo_entries()``. That path had +no coverage at all, which is how the bug below shipped. + +``read_topo_file`` selects the file subset three ways, in priority order: +descriptor ``crop_bounds``, then ``topo_crop_extent``, then the AMR domain. The +buffer expansion (``nbuf4``) was applied on the second branch only, so any file +carrying ``crop_bounds`` -- which is *every* file ``topo_entries()`` writes -- +silently discarded its declared ``buffer``. + +The failure is not subtle once you look for it: GeoClaw refuses to start with +"topo arrays do not cover domain" for a domain the user correctly sized against +crop-plus-buffer. The two cases here are the before/after of exactly that. + +Each case asserts on the topo grid GeoClaw reports in ``fort.geo`` -- the +window actually loaded -- rather than only on whether the run survived, so an +over-buffered or ignored crop fails as loudly as a dropped one. + +NetCDF is required; the whole module skips without it. +""" + +import re +import shutil +import subprocess +from pathlib import Path + +import numpy as np +import pytest + +import clawpack.geoclaw.test as gtest +import clawpack.geoclaw.topotools as topotools +from clawpack.geoclaw import netcdf_utils as ncutils +from clawpack.geoclaw.data import TopographyData + +testdir = Path(__file__).parent + +# Topo file, generously larger than any crop below. DELTA is a negative power +# of two so every coordinate is exact in binary: with 0.05 the file's own +# spacing came out as 4.9999999999997e-2 and the crop_bounds comparison +# (an inclusive >= / <= on the raw coordinates) landed one index differently on +# each edge, which made the window asymmetric for reasons having nothing to do +# with what is under test. +TOPO_X = (-100.0, -80.0) +TOPO_Y = (20.0, 30.0) +DELTA = 0.125 + +# The crop window written into the descriptor, and the buffer requested. All +# four edges are integer multiples of DELTA away from the file origin. +CROP = (-95.0, -90.0, 22.0, 25.0) +BUFFER = 4 +BUFFER_DEG = BUFFER * DELTA # 0.5 degrees on each side + +# Margin used when placing a test domain relative to a topo window. The domain +# must not merely *touch* the topo edge: GeoClaw's coverage test is an +# inequality on floats, and a domain whose bounds coincide exactly with the topo +# bounds passes it and then crashes intermittently in the boundary +# interpolation. Two cells of slack keeps these tests about the buffer. +MARGIN = 2 * DELTA + +pytestmark = pytest.mark.regression + + +def _netcdf_build_vars(): + """make_vars enabling a NetCDF build, or None if unavailable. + + Same probe as the met_forcing suite: ``nf-config`` for the Fortran flags, + plus ``nc-config --libs`` because ``nf-config --flibs`` names ``-lnetcdf`` + without its ``-L`` path. + """ + nf = shutil.which("nf-config") + nc = shutil.which("nc-config") + if nf is None: + return None + try: + fflags = subprocess.check_output([nf, "--fflags"], text=True).strip() + flibs = subprocess.check_output([nf, "--flibs"], text=True).strip() + if nc is not None: + flibs += " " + subprocess.check_output( + [nc, "--libs"], text=True).strip() + except (subprocess.CalledProcessError, OSError): + return None + return {"USE_NETCDF": "1", "NETCDF_FFLAGS": fflags, + "NETCDF_LFLAGS": flibs} + + +@pytest.fixture(scope="module") +def netcdf_xgeoclaw(tmp_path_factory): + """Build the NetCDF-enabled ``xgeoclaw`` once for the module.""" + make_vars = _netcdf_build_vars() + if make_vars is None: + pytest.skip("NetCDF (nf-config/nc-config) unavailable; the topo " + "descriptor path is topo_type=4 only.") + build_dir = tmp_path_factory.mktemp("topo_crop_build") + builder = gtest.GeoClawTestRunner(build_dir, test_path=testdir) + builder.build_executable(make_vars=make_vars) + return build_dir / builder.executable_name + + +def _write_nc_topo(path): + """A CF-compliant NetCDF bathymetry file spanning TOPO_X x TOPO_Y. + + The field is a plain south-to-north ramp: nothing here tests values, only + which subset of them Fortran loads. + """ + netCDF4 = pytest.importorskip("netCDF4") + + x = np.arange(TOPO_X[0], TOPO_X[1] + 1e-9, DELTA) + y = np.arange(TOPO_Y[0], TOPO_Y[1] + 1e-9, DELTA) + Z = -100.0 + np.outer(np.linspace(0.0, 90.0, y.size), np.ones_like(x)) + + with netCDF4.Dataset(path, 'w') as ds: + ds.createDimension('lon', x.size) + ds.createDimension('lat', y.size) + v = ds.createVariable('lon', 'f8', ('lon',)) + v[:] = x + v.units = 'degrees_east' + v.standard_name = 'longitude' + v = ds.createVariable('lat', 'f8', ('lat',)) + v[:] = y + v.units = 'degrees_north' + v.standard_name = 'latitude' + v = ds.createVariable('elevation', 'f8', ('lat', 'lon')) + v[:] = Z + v.units = 'm' + v.standard_name = 'height_above_mean_sea_level' + v.positive = 'up' + ds.Conventions = 'CF-1.8' + return path + + +def _write_topo_data(out_path, nc_path, buffer): + """Write topo.data via topo_entries(), so the descriptor has crop_bounds. + + Going through ``topo_entries()`` rather than setting ``crop_extent`` alone + is the point: it is what puts a ``crop_bounds`` line in the descriptor, and + therefore what selects the Fortran branch under test. + """ + with ncutils.TopoInspector(str(nc_path), crop_bounds=CROP) as insp: + entries = insp.topo_entries() + + topos = [] + for _topo_type, _path, meta in entries: + t = topotools.Topography() + t.path = str(nc_path) + t.topo_type = 4 + t.crop_extent = list(CROP) + t.buffer = buffer + t._netcdf_meta = meta + topos.append(t) + + td = TopographyData() + td.topofiles = topos + td.write(out_file=str(out_path)) + return len(entries) + + +def _ring_domain(): + """A domain lying in the buffer ring: outside CROP, inside CROP+buffer. + + Both gaps are MARGIN wide, so the case is decided by whether the buffer was + applied and not by float comparisons at a coincident edge. + """ + return (CROP[0] - MARGIN, CROP[1] + MARGIN, + CROP[2] - MARGIN, CROP[3] + MARGIN) + + +def _topo_window(tmp_path): + """Parse the topo grid GeoClaw actually loaded out of ``fort.geo``. + + read_topo_header writes ``mx = N x = (lo,hi)`` / ``my = ...`` for each + file; with one topo file the first pair is ours. + """ + text = (tmp_path / "fort.geo").read_text() + m = re.search(r"mx\s*=\s*(\d+)\s+x\s*=\s*\(\s*([-\d.E+]+)\s*," + r"\s*([-\d.E+]+)\s*\)", text) + n = re.search(r"my\s*=\s*(\d+)\s+y\s*=\s*\(\s*([-\d.E+]+)\s*," + r"\s*([-\d.E+]+)\s*\)", text) + assert m is not None and n is not None, f"no topo grid in fort.geo:\n{text}" + return (float(m.group(2)), float(m.group(3)), + float(n.group(2)), float(n.group(3)), + int(m.group(1)), int(n.group(1))) + + +def _run(tmp_path, prebuilt, domain, buffer): + """Set up and run one case; return (returncode, stdout).""" + nc_path = tmp_path / "topo.nc" + _write_nc_topo(nc_path) + + runner = gtest.GeoClawTestRunner(tmp_path, test_path=testdir) + runner.set_data() + cd = runner.rundata.clawdata + cd.lower = [domain[0], domain[2]] + cd.upper = [domain[1], domain[3]] + # Cell count is immaterial to the topo-coverage check; keep it small but + # non-degenerate. + cd.num_cells = [20, 12] + runner.write_data() + + # write_data() emits its own topo.data; overwrite it with the + # descriptor version, which is what this test is actually about. + n_entries = _write_topo_data(tmp_path / "topo.data", nc_path, buffer) + assert n_entries == 1, "crop lies wholly inside the file; expected 1 entry" + + shutil.copy(prebuilt, tmp_path / runner.executable_name) + (tmp_path / "_output").mkdir(exist_ok=True) + proc = subprocess.run([str(tmp_path / runner.executable_name)], + cwd=tmp_path, capture_output=True, text=True) + # Kept for post-mortem: pytest truncates a long assertion message. + (tmp_path / "run.log").write_text(proc.stdout + proc.stderr) + return proc.returncode, proc.stdout + proc.stderr + + +def test_descriptor_crop_honors_buffer(tmp_path, netcdf_xgeoclaw): + """A domain inside crop+buffer but outside crop alone must run. + + This is the regression. With ``nbuf4`` pinned to 0 on the ``crop_bounds`` + branch the loaded topo covered only CROP, leaving the outer 0.20-degree ring + empty, and GeoClaw stopped with "topo arrays do not cover domain" -- an + abort, for inputs that were correctly specified. + """ + domain = _ring_domain() + rc, out = _run(tmp_path, netcdf_xgeoclaw, domain, BUFFER) + + assert "topo arrays do not cover domain" not in out, ( + "buffer was not applied to the descriptor crop; the topo loaded " + f"covers only {CROP}.\n{out}") + assert rc == 0, out + + # Stronger than "it ran": the loaded window must be CROP grown by exactly + # BUFFER points on every side. A test that only checks for the abort would + # also pass if the reader over-buffered, or ignored crop_bounds and loaded + # the whole file. + xll, xhi, yll, yhi, mx, my = _topo_window(tmp_path) + assert xll == pytest.approx(CROP[0] - BUFFER_DEG) + assert xhi == pytest.approx(CROP[1] + BUFFER_DEG) + assert yll == pytest.approx(CROP[2] - BUFFER_DEG) + assert yhi == pytest.approx(CROP[3] + BUFFER_DEG) + assert mx == round((CROP[1] - CROP[0]) / DELTA) + 2 * BUFFER + 1 + assert my == round((CROP[3] - CROP[2]) / DELTA) + 2 * BUFFER + 1 + + +def test_descriptor_crop_without_buffer_does_not_cover_ring( + tmp_path, netcdf_xgeoclaw): + """The complement: with buffer=0 the same domain is genuinely uncovered. + + Without this, the test above would still pass if the reader started + ignoring ``crop_bounds`` entirely and loaded the whole file -- which would + be a different bug with the same symptom. Here the coverage failure is the + correct answer, and seeing it proves the crop is being applied at all. + """ + domain = _ring_domain() + rc, out = _run(tmp_path, netcdf_xgeoclaw, domain, 0) + + assert "topo arrays do not cover domain" in out, ( + "buffer=0 should leave the ring outside crop_bounds uncovered; the " + f"reader appears not to be applying crop_bounds at all.\n{out}") + + +def test_descriptor_crop_inside_window_runs_either_way( + tmp_path, netcdf_xgeoclaw): + """Control: a domain wholly inside CROP runs with no buffer at all. + + Pins that the descriptor crop itself is loaded correctly, independently of + the buffer question. + """ + inset = 4 * DELTA + domain = (CROP[0] + inset, CROP[1] - inset, + CROP[2] + inset, CROP[3] - inset) + + rc, out = _run(tmp_path, netcdf_xgeoclaw, domain, 0) + + assert "topo arrays do not cover domain" not in out, out + assert rc == 0, out + + # The window is CROP itself: unbuffered, and not the whole file. + xll, xhi, yll, yhi, _mx, _my = _topo_window(tmp_path) + assert (xll, xhi, yll, yhi) == pytest.approx( + (CROP[0], CROP[1], CROP[2], CROP[3])) diff --git a/tests/test_topotools_preprocessing.py b/tests/test_topotools_preprocessing.py index 4cba32032..50f2d55d5 100644 --- a/tests/test_topotools_preprocessing.py +++ b/tests/test_topotools_preprocessing.py @@ -1202,7 +1202,11 @@ def test_crop_no_overlap_keeps_full_grid(nc_topo_path): t = Topography() t.crop_extent = [100.0, 200.0, 100.0, 200.0] - t.read(path, topo_type=4) + # The fall-back is kept, but it must not be silent: the Fortran reader + # treats the same condition as fatal, so a run that ignores this here + # fails there instead. + with pytest.warns(UserWarning, match="did not overlap"): + t.read(path, topo_type=4) np.testing.assert_array_equal(t.Z, ref.Z) @@ -1302,3 +1306,226 @@ def test_coarsen_align_lattice_invariant(tmp_path): yphase = (tk.y[0] - align[1]) / coarsen assert abs(xphase - round(xphase)) < 1e-9, (k, tk.x[0]) assert abs(yphase - round(yphase)) < 1e-9, (k, tk.y[0]) + + +# =========================================================================== +# Group N — Antimeridian and degenerate crop windows +# +# Nothing exercised antimeridian cropping through Topography before these, +# which is how the failures below shipped silently. The point of the group is +# that a Topography *never wraps*: a crop is either an ordinary ascending +# window, or it is not expressible on a single Topography at all and has to go +# through TopoInspector.topo_entries(). +# =========================================================================== + +_WRAPPED_CROP = [170.0, -170.0, -5.0, 5.0] # crosses the seam, descending +_CONTINUOUS_CROP = [-211.0, -99.0, -5.0, 5.0] # same region, continuous coords + + +def _write_global_tt3(path: Path) -> Path: + """A 1-degree global file spanning the antimeridian, x in [-180, 180].""" + x = np.linspace(-180.0, 180.0, 361) + y = np.linspace(-10.0, 10.0, 21) + Z = -1000.0 + 10.0 * np.cos(np.radians(x))[None, :] * np.ones((y.size, 1)) + t = Topography() + t.set_xyZ(x, y, Z) + t.write(str(path), topo_type=3) + return path + + +@pytest.fixture +def global_tt3_path(tmp_path): + return _write_global_tt3(tmp_path / "global.tt3") + + +def test_wrapped_crop_spelling_raises_not_empty_grid(global_tt3_path): + """[170, -170] used to produce an *empty* Topography (Z.shape == (11, 0)) + whose .extent then raised an opaque "zero-size array to reduction" from + numpy, far from the cause.""" + t = Topography() + t.crop_extent = list(_WRAPPED_CROP) + with pytest.raises(ValueError, match="must increase in both coordinates"): + t.read(str(global_tt3_path), topo_type=3) + + +def test_wrapped_crop_error_names_the_supported_route(global_tt3_path): + """The error has to say what to do instead, or it just moves the confusion.""" + t = Topography() + t.crop_extent = list(_WRAPPED_CROP) + with pytest.raises(ValueError) as excinfo: + t.read(str(global_tt3_path), topo_type=3) + assert "topo_entries" in str(excinfo.value) + + +def test_descending_latitude_crop_also_raises(global_tt3_path): + """A descending *latitude* pair has no antimeridian excuse at all; it was + equally silent.""" + t = Topography() + t.crop_extent = [-100.0, -80.0, 5.0, -5.0] + with pytest.raises(ValueError, match="must increase in both coordinates"): + t.read(str(global_tt3_path), topo_type=3) + + +def test_continuous_crop_is_clipped_and_says_so(global_tt3_path): + """The continuous spelling is accepted, but it is *not* wrapped -- it is + reduced to the part of the file that exists (112 degrees requested, 81 + delivered). That silent reduction is the whole antimeridian confusion, so + it must warn.""" + t = Topography() + t.crop_extent = list(_CONTINUOUS_CROP) + with pytest.warns(UserWarning, match="extends past the data"): + t.read(str(global_tt3_path), topo_type=3) + + # Clipped to the file's western edge, not wrapped around to +149. + assert float(t.x[0]) == pytest.approx(-180.0) + assert float(t.x[-1]) == pytest.approx(-99.0) + + +def test_ordinary_crop_does_not_warn(global_tt3_path): + """The clipping warning must not fire for a crop wholly inside the file, + or it becomes noise everyone filters.""" + t = Topography() + t.crop_extent = [-100.0, -80.0, -5.0, 5.0] + with warnings.catch_warnings(): + warnings.simplefilter("error") + t.read(str(global_tt3_path), topo_type=3) + assert float(t.x[0]) == pytest.approx(-100.0) + + +def test_crop_between_grid_points_raises(global_tt3_path): + """A window inside the extent but narrower than one cell contains no data. + It used to return the *full file* -- the opposite of what was asked.""" + t = Topography() + t.crop_extent = [10.2, 10.8, -5.0, 5.0] + with pytest.raises(ValueError, match="lies between grid points"): + t.read(str(global_tt3_path), topo_type=3) + + +def test_crop_no_overlap_ascii_warns_and_keeps_full_grid(global_tt3_path): + """ASCII counterpart of test_crop_no_overlap_keeps_full_grid.""" + t = Topography() + t.crop_extent = [300.0, 320.0, -5.0, 5.0] + with pytest.warns(UserWarning, match="did not overlap"): + t.read(str(global_tt3_path), topo_type=3) + assert t.x.size == 361 + + +def test_unstructured_with_crop_raises_not_typeerror(tmp_path): + """This used to die with `TypeError: list indices must be integers` from + indexing a Python list as an array, several frames from the cause. crop() + already refused unstructured input; read() now agrees. + + The fixture is a genuine 3-column xyz file (not a headed .tt3) so that the + read gets far enough to hit that bug when the guard is removed -- a test + whose "before" failure is an unrelated parse error would pin nothing. + """ + xyz = tmp_path / "scattered.xyz" + with open(xyz, "w") as f: + for x in np.linspace(-110.0, -70.0, 9): + for y in np.linspace(-8.0, 8.0, 5): + f.write(f"{x} {y} {-1000.0 + x + y}\n") + + t = Topography() + t.crop_extent = [-100.0, -80.0, -5.0, 5.0] + with pytest.raises(NotImplementedError, match="unstructured"): + t.read(str(xyz), topo_type=1, unstructured=True) + + +def test_cross_seam_crop_raises_at_topo_data_write(global_tt3_path, tmp_path): + """Writing a descending crop_extent used to emit `crop_bounds = 170.0 + -170.0`, which Fortran resolves to mx=0, my=0: an empty topo, no error.""" + t = Topography() + t.path = str(global_tt3_path) + t.topo_type = 3 + t.crop_extent = list(_WRAPPED_CROP) + + td = TopographyData() + td.topofiles = [t] + with pytest.raises(ValueError, match="crosses the antimeridian"): + td.write(out_file=str(tmp_path / "topo.data")) + + +def test_topo_type_none_inferred_from_suffix(global_tt3_path, tmp_path): + """topo_type=None reached the `:3d` format and raised a TypeError naming + neither the file nor the attribute. A .tt3 suffix is unambiguous.""" + t = Topography() + t.path = str(global_tt3_path) + t.topo_type = None + + td = TopographyData() + td.topofiles = [t] + out = tmp_path / "topo.data" + td.write(out_file=str(out)) + assert t.topo_type == 3 + assert " 3 # topo_type" in out.read_text() + + +def test_topo_type_none_unknown_suffix_raises(tmp_path): + """When the suffix carries no type either, say which attribute to set.""" + src = _write_global_tt3(tmp_path / "global.tt3") + unknown = tmp_path / "global.dat" + unknown.write_bytes(src.read_bytes()) + + t = Topography() + t.path = str(unknown) + t.topo_type = None + + td = TopographyData() + td.topofiles = [t] + with pytest.raises(ValueError, match="topo_type is not set"): + td.write(out_file=str(tmp_path / "topo.data")) + + +@pytest.mark.netcdf +def test_cropped_netcdf_read_after_read_header_updates_extent(tmp_path): + """read_header() populates _extent from the file header; the topo_type=4 + read applies its crop while reading the hyperslab, bypassing the property + setters that would invalidate it. The object then reported the *full file* + extent alongside cropped data -- and .extent is what _compute_priority_order + and the plotting routines use.""" + pytest.importorskip("xarray") + pytest.importorskip("netCDF4") + path, _, _ = tmp_path / "nc_extent.nc", None, None + _make_nc_topo(path, "lon", "lat") + + t = Topography() + t.path = str(path) + t.topo_type = 4 + t.read_header() + assert list(t.extent) == pytest.approx([_ORIGIN_X, _ORIGIN_X + _NX - 1, + _ORIGIN_Y, _ORIGIN_Y + _NY - 1]) + + t.crop_extent = [_ORIGIN_X + 2, _ORIGIN_X + 5, + _ORIGIN_Y + 2, _ORIGIN_Y + 5] + t.read() + + assert list(t.extent) == pytest.approx( + [float(t.x[0]), float(t.x[-1]), float(t.y[0]), float(t.y[-1])]) + assert float(t.extent[0]) == pytest.approx(_ORIGIN_X + 2) + + +def test_buffer_and_coarsen_give_absolute_output_shape(tt2_path): + """The existing combined test is a *relative* netCDF-vs-ASCII equality, so + nothing pinned whether buffer=1, coarsen=2 adds 1 or 2 output points per + side. buffer counts coarsened *output* points: the window is expanded by + buffer*coarsen native points before the strided slice.""" + ref = Topography() + ref.crop_extent = [_ORIGIN_X + 2, _ORIGIN_X + 7, _ORIGIN_Y + 2, + _ORIGIN_Y + 7] + ref.coarsen = 2 + ref.read(str(tt2_path), topo_type=2) + + t = Topography() + t.crop_extent = list(ref.crop_extent) + t.coarsen = 2 + t.buffer = 1 + t.read(str(tt2_path), topo_type=2) + + # One extra coarsened point on each side of each axis. + assert t.x.size == ref.x.size + 2 + assert t.y.size == ref.y.size + 2 + assert t.Z.shape == (ref.Z.shape[0] + 2, ref.Z.shape[1] + 2) + # Coarsening is unchanged by the buffer: still every 2nd native point. + assert float(t.x[1] - t.x[0]) == pytest.approx(2.0 * _DELTA) + # And the buffered window still contains the unbuffered one. + assert float(t.x[0]) == pytest.approx(float(ref.x[0]) - 2.0 * _DELTA) From 256fcd5258ab5125686f0f83caefddf7a793e9fd Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Fri, 4 Sep 2026 21:43:35 -0400 Subject: [PATCH 2/7] Make remote and cross-seam topo sources usable from setrun A URL in topofiles was mangled by os.path.abspath into a bogus local path before the reader's existing URL guard could see it; it is now rejected with the fetch_remote_topo recipe. A crop crossing the antimeridian raised 'crop_bounds exceed file extent' even though _compute_lon_entries already covered it; TopographyData.write now resolves entries before writing and splits such a crop into one descriptor entry per side, carrying buffer and coarsen onto each so the Fortran buffer fix applies. Adds byte-exact topo.data goldens and the first end-to-end test of two entries with differing lon_wrap_offset. Signed-off-by: Kyle Mandli Assisted-by: claude claude-opus-5[1m] --- src/python/geoclaw/data.py | 284 ++++++++++++------ src/python/geoclaw/netcdf_utils.py | 70 ++++- src/python/geoclaw/topotools.py | 47 +-- .../ascii_all_preprocessing.txt | 20 ++ tests/data/topo_data_golden/ascii_single.txt | 20 ++ .../data/topo_data_golden/netcdf_cropped.txt | 30 ++ .../data/topo_data_golden/netcdf_no_crop.txt | 29 ++ .../topo_data_golden/two_files_priority.txt | 31 ++ .../topo_crop/test_topo_descriptor_crop.py | 134 +++++++++ tests/test_topo_data_golden.py | 204 +++++++++++++ tests/test_topotools_preprocessing.py | 263 +++++++++++++++- 11 files changed, 1013 insertions(+), 119 deletions(-) create mode 100644 tests/data/topo_data_golden/ascii_all_preprocessing.txt create mode 100644 tests/data/topo_data_golden/ascii_single.txt create mode 100644 tests/data/topo_data_golden/netcdf_cropped.txt create mode 100644 tests/data/topo_data_golden/netcdf_no_crop.txt create mode 100644 tests/data/topo_data_golden/two_files_priority.txt create mode 100644 tests/test_topo_data_golden.py diff --git a/src/python/geoclaw/data.py b/src/python/geoclaw/data.py index 488b31e21..a6edb3c1f 100755 --- a/src/python/geoclaw/data.py +++ b/src/python/geoclaw/data.py @@ -154,6 +154,28 @@ def write(self,data_source='setrun.py', out_file='refinement.data'): +def _reject_remote_path(path, kind, fetch_hint): + """Raise if *path* is a URL rather than a local file. + + GeoClaw's Fortran reader opens a filesystem path, and the ``.data`` writers + run every path through ``os.path.abspath``, which turns a URL into a bogus + local path (see ``netcdf_utils.is_remote_url``). The result used to be a + FileNotFoundError naming a path the user never typed, which gives no hint + that the real problem is "remote sources must be fetched first". + + *kind* is the noun for the message ("topography"/"dtopography") and + *fetch_hint* the recipe to show. + """ + from clawpack.geoclaw.netcdf_utils import is_remote_url + + if is_remote_url(path): + raise ValueError( + f"{kind} path is a URL, which GeoClaw's Fortran reader cannot " + f"open:\n {path}\n" + f"Read it in Python first and write a local file, then reference " + f"that:\n{fetch_hint}") + + def _write_preprocessing_block(f, t): """Write the 8 preprocessing-attribute lines for one topo/dtopo file. @@ -345,6 +367,157 @@ def _cell_area(topo): # keeps the sort stable for equal-area files (they retain input order). return sorted(topos, key=_cell_area, reverse=True) + def _resolve_topo_records(self, topos, out_file): + """Expand *topos* into the entries that will be written to topo.data. + + Returns a list of ``(fname, topo_type, topo, meta)``. ``meta`` is a + ``TopoMetadata`` for type-4 files (written as a CF descriptor block) + and None otherwise. + + Most files produce exactly one record. A ``topo_type=4`` file whose + ``crop_extent`` runs past the file's longitude extent produces one + record per side of the antimeridian, each with its own + ``lon_wrap_offset`` -- this is what makes a cross-seam crop work from + an ordinary ``topofiles.append(topo)`` instead of requiring the caller + to build descriptors by hand. + """ + import dataclasses + from clawpack.geoclaw import netcdf_utils as _ncutils + from clawpack.geoclaw import topotools + + records = [] + for topo in topos: + _reject_remote_path( + topo.path, "Topography", + " topo = topotools.fetch_remote_topo(\n" + " url, crop_extent=[...], coarsen=..., buffer=...)\n" + " topo.write('topo_cropped.tt3', topo_type=3)\n" + " rundata.topo_data.topofiles.append(topo)\n" + "Only the requested hyperslab is read, so this does not " + "download the whole file.") + + # Resolve path relative to out_file's directory, same as before + fname = os.path.abspath( + os.path.join(os.path.dirname(out_file), topo.path)) + + # topo_type may still be None when a Topography was built by hand + # and never read; the :3d format would raise a TypeError naming + # neither the file nor the attribute. + topo_type = topo.topo_type + if topo_type is None: + topo_type = topotools.determine_topo_type(topo.path, + default=None) + if topo_type is None: + raise ValueError( + f"topo_type is not set for {topo.path} and cannot be " + f"inferred from its extension. Set topo.topo_type " + f"explicitly, or pass topo_type= to Topography().") + topo.topo_type = topo_type + + # Type 5 (GeoTIFF) is readable in Python but has no case(5) in the + # Fortran reader, which aborts with "Unrecognized topo_type". + if abs(topo_type) == 5: + warnings.warn( + f"{topo.path} is written to {out_file} as topo_type=5 " + f"(GeoTIFF). GeoClaw's Fortran reader does not support " + f"type 5 and will abort; convert to topo_type 3 or 4 " + f"(Topography.write) before running.", + UserWarning, stacklevel=3) + + # A descending crop_extent is not a rectangle in any frame. The + # *continuous* spelling of a cross-seam crop ([-190, -120]) is + # handled below; the wrapped spelling ([170, -170]) cannot be, and + # written directly it would emit "crop_bounds = 170.0 -170.0", + # which Fortran resolves to mx=0, my=0 -- an empty topo, silently. + if (topo.crop_extent is not None + and topo.crop_extent[0] >= topo.crop_extent[1] + and getattr(topo, '_netcdf_meta', None) is None): + raise ValueError( + f"crop_extent {list(topo.crop_extent)} for {topo.path} is " + f"descending in longitude. To cross the antimeridian, " + f"write the crop in continuous coordinates instead -- " + f"[{topo.crop_extent[0] - 360.0}, {topo.crop_extent[1]}] " + f"for this one -- which is split across the seam " + f"automatically. A descending pair has no such reading and " + f"would produce an empty grid in Fortran with no error.") + + if abs(topo_type) != 4: + records.append((fname, topo_type, topo, None)) + continue + + # --- type 4: build the CF descriptor block ------------------- + if getattr(topo, '_netcdf_meta', None) is not None: + # Pre-computed by topo_entries(); already carries the right + # lon_wrap_offset and file-coordinate crop_bounds. + records.append((fname, topo_type, topo, topo._netcdf_meta)) + continue + + crop = (tuple(float(v) for v in topo.crop_extent) + if topo.crop_extent is not None else None) + + # Opened without crop_bounds so the longitude extent can be read + # before deciding whether the crop needs wrapping; setting them at + # construction would validate (and reject) a cross-seam crop first. + with _ncutils.TopoInspector(fname) as insp: + if insp.var_name is None: + insp.var_name = insp._find_topo_var_name() + + needs_wrap = False + if crop is not None: + x_name = insp._find_x_name() + lon = insp.ds[x_name].values + file_lon_min = float(lon.min()) + file_lon_max = float(lon.max()) + tol = 1e-9 + needs_wrap = (crop[0] < file_lon_min - tol + or crop[1] > file_lon_max + tol) + + if not needs_wrap: + # Unchanged path: validates crop_bounds (including + # latitude) and skips the expensive fill scan, exactly as + # before. + insp.crop_bounds = crop + file_meta = insp.inspect(insp.var_name) + src_units = insp._check_topo_units() + scale = _ncutils._units_scale( + src_units, _ncutils.GEOCLAW_NETCDF_UNITS['topo']) + records.append((fname, topo_type, topo, + _ncutils.TopoMetadata( + **dataclasses.asdict(file_meta), + var_name=insp.var_name, + source_units=src_units, + scale_factor=scale, + fill_action='abort', + lon_wrap_offset=0.0))) + continue + + # Cross-seam (or wholly off-seam) crop: one entry per side, + # each with its own lon_wrap_offset. fill_scan=False because + # topo_entries inspects with crop_bounds unset, which would + # otherwise scan the whole file -- ruinous for a global DEM + # and prone to rejecting NaN far outside the crop. + y_name = insp._find_y_name() + lat = insp.ds[y_name].values + file_lat_min = float(lat.min()) + file_lat_max = float(lat.max()) + if (crop[2] < file_lat_min - 1e-9 + or crop[3] > file_lat_max + 1e-9): + # Longitude wraps; latitude never does. + raise ValueError( + f"crop_extent {list(topo.crop_extent)} for {topo.path} " + f"exceeds the file's latitude extent " + f"[{file_lat_min}, {file_lat_max}]. Longitude is " + f"wrapped across the antimeridian, but latitude cannot " + f"be; narrow the requested latitude range.") + + insp.crop_bounds = crop + entries = insp.topo_entries(fill_scan=False) + + for _entry_type, _entry_path, entry_meta in entries: + records.append((fname, topo_type, topo, entry_meta)) + + return records + def write(self, data_source='setrun.py', out_file='topo.data'): self.open_data_file(out_file, data_source) @@ -371,96 +544,28 @@ def write(self, data_source='setrun.py', out_file='topo.data'): UserWarning, stacklevel=2, ) - self.data_write(value=len(topos), alt_name='ntopofiles') - for topo in topos: - # Resolve path relative to out_file's directory, same as before - fname = os.path.abspath( - os.path.join(os.path.dirname(out_file), topo.path)) - - f = self._out_file - - # topo_type may still be None when a Topography was built by - # hand and never read; the :3d format below would raise a - # TypeError naming neither the file nor the attribute. - _topo_type = topo.topo_type - if _topo_type is None: - from clawpack.geoclaw import topotools - _topo_type = topotools.determine_topo_type( - topo.path, default=None) - if _topo_type is None: - raise ValueError( - f"topo_type is not set for {topo.path} and cannot " - f"be inferred from its extension. Set topo.topo_type " - f"explicitly, or pass topo_type= to Topography().") - topo.topo_type = _topo_type - - # Type 5 (GeoTIFF) is readable in Python but has no case(5) in - # the Fortran reader, which aborts with "Unrecognized topo_type". - if abs(_topo_type) == 5: - warnings.warn( - f"{topo.path} is written to {out_file} as topo_type=5 " - f"(GeoTIFF). GeoClaw's Fortran reader does not support " - f"type 5 and will abort; convert to topo_type 3 or 4 " - f"(Topography.write) before running.", - UserWarning, stacklevel=2) - - # A crop that crosses the antimeridian cannot be expressed as a - # single descriptor: crop_bounds would read "170.0 -170.0", - # which Fortran resolves to mx=0, my=0 -- an empty topo file, - # silently. TopoInspector.topo_entries() splits such a crop - # into two entries with the appropriate lon_wrap_offset. - if (topo.crop_extent is not None - and topo.crop_extent[0] >= topo.crop_extent[1] - and getattr(topo, '_netcdf_meta', None) is None): - raise ValueError( - f"crop_extent {list(topo.crop_extent)} for {topo.path} " - f"is descending in longitude, i.e. it crosses the " - f"antimeridian. Written directly it would produce an " - f"empty grid in Fortran with no error. Build the entries " - f"with TopoInspector.topo_entries(), which emits one " - f"entry per side of the seam with the correct " - f"lon_wrap_offset, and append those to topofiles.") + # Resolve first, write second. A type-4 crop that runs off the + # file's longitude extent is covered by *two* descriptor entries + # (one per side of the seam), so the entry count is not known until + # every file has been resolved -- and ntopofiles is written before + # the blocks. + records = self._resolve_topo_records(topos, out_file) + self.data_write(value=len(records), alt_name='ntopofiles') + f = self._out_file + for fname, topo_type, topo, meta in records: f.write(f"\n'{fname}' # topo_path\n") - f.write(f"{_topo_type:3d} # topo_type\n") + f.write(f"{topo_type:3d} # topo_type\n") + # The originating Topography is reused for every entry it + # expanded into, so buffer/coarsen/align/shifts reach Fortran + # for each one. (crop_extent is written as the user gave it; + # for type 4 the descriptor's crop_bounds takes priority, and + # for a wrapped pair it is the per-entry crop_bounds that + # differ.) _write_preprocessing_block(f, topo) - - # For NetCDF (type 4): write the CF descriptor block that - # Fortran's read_netcdf_descriptor parses (key=value lines - # terminated by a blank line). Uses inspect() rather than - # inspect_topo() to avoid an expensive full-file fill scan. - if abs(topo.topo_type) == 4: - import dataclasses + if meta is not None: from clawpack.geoclaw import netcdf_utils as _ncutils - if getattr(topo, '_netcdf_meta', None) is not None: - # Pre-computed metadata from topo_entries() — already has - # correct lon_wrap_offset and file-coordinate crop_bounds. - _meta = topo._netcdf_meta - else: - _crop = (tuple(topo.crop_extent) - if topo.crop_extent is not None else None) - with _ncutils.TopoInspector( - fname, crop_bounds=_crop - ) as _insp: - if _insp.var_name is None: - _insp.var_name = _insp._find_topo_var_name() - _file_meta = _insp.inspect(_insp.var_name) - # A recognized non-meter unit yields a scale_factor - # Fortran applies on read (missing/unrecognized - # units still raise). - _src_units = _insp._check_topo_units() - _scale = _ncutils._units_scale( - _src_units, - _ncutils.GEOCLAW_NETCDF_UNITS['topo']) - _meta = _ncutils.TopoMetadata( - **dataclasses.asdict(_file_meta), - var_name=_insp.var_name, - source_units=_src_units, - scale_factor=_scale, - fill_action='abort', - lon_wrap_offset=0.0, - ) - _ncutils.DescriptorWriter.write_topo_descriptor(f, _meta) + _ncutils.DescriptorWriter.write_topo_descriptor(f, meta) elif self.test_topography == 1: self.data_write(name='topo_location',description='(Bathymetry jump location)') @@ -703,6 +808,13 @@ def write(self, data_source='setrun.py', out_file='dtopo.data'): "dtopography (file %s). Only x_shift, y_shift, z_shift and " "negate_z are supported." % (", ".join(unsupported), d.path)) + _reject_remote_path( + d.path, "DTopography", + " dtopo = dtopotools.DTopography()\n" + " dtopo.read(url, dtopo_type=4)\n" + " dtopo.write('dtopo_local.tt3', dtopo_type=3)\n" + " rundata.dtopo_data.dtopofiles.append(dtopo)") + # if path is relative in setrun, assume it's relative to the # same directory that out_file comes from fname = os.path.abspath( diff --git a/src/python/geoclaw/netcdf_utils.py b/src/python/geoclaw/netcdf_utils.py index 03ebad6d5..24c30a194 100644 --- a/src/python/geoclaw/netcdf_utils.py +++ b/src/python/geoclaw/netcdf_utils.py @@ -383,6 +383,35 @@ class MetMetadata(FileMetadata): # Base inspector # --------------------------------------------------------------------------- +def is_remote_url(path) -> bool: + """True if *path* is a remote URL rather than a local filesystem path. + + Remote OPeNDAP/THREDDS URLs must be kept as strings all the way to xarray. + Passing one through ``pathlib.Path`` collapses ``"https://"`` to + ``"https:/"`` and makes it *relative*, and ``os.path.abspath`` then resolves + that against the cwd -- turning + + https://www.ngdc.noaa.gov/thredds/dodsC/.../ETOPO_2022.nc + + into + + /your/run/directory/https:/www.ngdc.noaa.gov/thredds/.../ETOPO_2022.nc + + which fails as a baffling FileNotFoundError naming a path the user never + typed. This was fixed once inside the NetCDF reader (PR #726); the same + trap exists anywhere a path is normalized, so the test lives here and is + shared rather than repeated. + + The regex is anchored on a URL scheme followed by "//", which excludes + Windows drive letters like ``C:\\data\\topo.nc`` (no "//"). It matches + any scheme, ``file://`` included -- that one is local, but it is still not + a path the Fortran reader can open, so callers that reject remote sources + should reject it too. + """ + return (isinstance(path, str) + and re.match(r"^[A-Za-z][A-Za-z0-9+.-]*://", path) is not None) + + class NetCDFInspector: """ Open a NetCDF file and inspect its coordinate metadata. @@ -406,13 +435,8 @@ def __init__( crop_bounds: Optional[tuple[float, float, float, float]] = None, ) -> None: # A remote OPeNDAP/THREDDS URL (e.g. "https://.../foo.nc") must reach - # xarray as a string. Wrapping it in pathlib.Path collapses "https://" - # to "https:/" and makes it a *relative* path, which the netCDF4 backend - # then resolves against the cwd -- producing a bogus local-file lookup - # (PR #726). The scheme-anchored regex ignores Windows drive paths - # like "C:\\..." (no "//"). - if isinstance(path, str) and re.match(r"^[A-Za-z][A-Za-z0-9+.-]*://", - path): + # xarray as a string; see is_remote_url() for why Path() breaks it. + if is_remote_url(path): self.path = path else: self.path = Path(path) @@ -950,9 +974,17 @@ def _check_fill_in_crop( # Public interface # ------------------------------------------------------------------ - def inspect_topo(self) -> TopoMetadata: + def inspect_topo(self, fill_scan: bool = True) -> TopoMetadata: """ Fully inspect the topo file and return a TopoMetadata instance. + + *fill_scan* controls the two data-reading checks (fill values and + elevation magnitude). They are the only part of this method that + touches the array rather than its metadata, and they cost a pass over + the current crop region -- or over the **whole file** when + ``crop_bounds`` is None, which for a remote global DEM means + downloading it. Pass False to skip them when the caller only needs + metadata; everything else is unaffected. """ if self.var_name is None: self.var_name = self._find_topo_var_name() @@ -962,9 +994,11 @@ def inspect_topo(self) -> TopoMetadata: # applies on read (missing/unrecognized units still raise). source_units = self._check_topo_units() scale_factor = _units_scale(source_units, GEOCLAW_NETCDF_UNITS['topo']) - self._check_fill_in_crop(base.x_name, base.y_name, base.y_increasing) - self._check_topo_magnitude(base.x_name, base.y_name, - base.y_increasing, scale_factor) + if fill_scan: + self._check_fill_in_crop(base.x_name, base.y_name, + base.y_increasing) + self._check_topo_magnitude(base.x_name, base.y_name, + base.y_increasing, scale_factor) return TopoMetadata( **dataclasses.asdict(base), @@ -975,7 +1009,7 @@ def inspect_topo(self) -> TopoMetadata: lon_wrap_offset=0.0, ) - def topo_entries(self) -> list[list]: + def topo_entries(self, fill_scan: bool = True) -> list[list]: """ Return a list of ready-to-use topo entries for topofiles. @@ -988,6 +1022,16 @@ def topo_entries(self) -> list[list]: converts them to file coordinates before storing in the returned metadata. Fortran can then use crop_bounds directly against file coordinate arrays before applying lon_wrap_offset. + + *fill_scan* is forwarded to :meth:`inspect_topo`. Note that the + inspection below runs with ``crop_bounds`` unset -- it has to, because + a wrapping crop lies outside the file extent by construction and would + fail validation -- so with ``fill_scan=True`` those checks scan the + **entire file** and reject NaN anywhere in it, not just in the crop. + For a global DEM that is expensive and usually wrong; callers that + only need the descriptor metadata (``TopographyData.write``) pass + False. Scoping the scan to each returned entry's own crop is the + better answer and is tracked for the topo-input refactor. """ # Interrogate without crop validation: self.crop_bounds is in domain @@ -995,7 +1039,7 @@ def topo_entries(self) -> list[list]: saved_crop = self.crop_bounds self.crop_bounds = None try: - meta = self.inspect_topo() + meta = self.inspect_topo(fill_scan=fill_scan) finally: self.crop_bounds = saved_crop diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py index eabfe8808..3f72d9d00 100644 --- a/src/python/geoclaw/topotools.py +++ b/src/python/geoclaw/topotools.py @@ -535,26 +535,37 @@ class Topography(object): most antimeridian confusion, so it is worth being concrete about the three cases: - 1. **Wrapped spelling, e.g. ``crop_extent=[170, -170, ...]``.** Rejected - with a ``ValueError``. There is no ascending window it could mean, and - silently producing an empty grid (which is what a descending pair used - to do) is worse than refusing. + 1. **Ordinary ascending crop inside the file.** The common case; nothing + special happens. 2. **Continuous spelling, e.g. ``crop_extent=[-211, -99, ...]`` for a file - on ``[-180, 180]``.** Accepted, but the part that lies off the file is - *clipped, not wrapped* -- you get ``[-180, -99]`` and a - ``UserWarning`` saying so. It is a legitimate way to say "as far west - as this file goes"; it is not a way to cross the seam. - - 3. **Genuinely crossing the seam.** Not expressible on a single - ``Topography``, and not expressible in a single ``topo.data`` entry - either. Use :meth:`netcdf_utils.TopoInspector.topo_entries`, which - returns one entry per side of the cut with the appropriate - ``lon_wrap_offset`` and file-coordinate ``crop_bounds``; append those to - ``rundata.topo_data.topofiles``. Writing a descending ``crop_extent`` - to ``topo.data`` directly raises, because Fortran would resolve - ``crop_bounds = 170.0 -170.0`` to ``mx=0, my=0``: an empty topography, - with no error. + on ``[-180, 180]``.** What this means depends on where it is used, and + the two answers are different on purpose: + + - **Reading into a Topography** (``read``, ``crop``): the part that lies + off the file is *clipped, not wrapped* -- you get ``[-180, -99]`` and + a ``UserWarning`` saying so. An in-memory ``Topography`` is one + ascending array; it has nowhere to put the wrapped part. + - **Writing to ``topo.data``** (``TopographyData.write``): the crop is + *split across the seam*. The writer emits one entry per side, each + with its own ``lon_wrap_offset`` and file-coordinate ``crop_bounds``, + and Fortran reassembles them into a single continuous region. + + So a cross-seam crop works from ordinary setrun code -- set + ``crop_extent`` and append the ``Topography`` -- with no descriptor + handling by the caller. ``buffer``, ``coarsen``, ``align`` and the + shifts are carried onto every entry. + + 3. **Wrapped spelling, ``crop_extent=[170, -170, ...]``.** Rejected + wherever it appears, including at ``topo.data`` write time: Fortran + would resolve ``crop_bounds = 170.0 -170.0`` to ``mx=0, my=0`` -- an + empty topography, with no error. Use the continuous spelling + (``[-190, -170]``), which case 2 handles. + + :meth:`netcdf_utils.TopoInspector.topo_entries` is the underlying + machinery, and is still available if you want the entries directly; the + writer now calls it for you. Latitude is never wrapped -- a crop whose + latitude runs off the file is an error, since there is no seam to cross. The general shape of it: **Python is the single-rectangle case; the wrap lives in the Fortran interface.** The same split explains why diff --git a/tests/data/topo_data_golden/ascii_all_preprocessing.txt b/tests/data/topo_data_golden/ascii_all_preprocessing.txt new file mode 100644 index 000000000..8016b6b06 --- /dev/null +++ b/tests/data/topo_data_golden/ascii_all_preprocessing.txt @@ -0,0 +1,20 @@ +######################################################## +### DO NOT EDIT THIS FILE: GENERATED AUTOMATICALLY #### +### To modify data, edit setrun.py #### +### and then "make .data" #### +######################################################## + +99999.0 =: topo_missing # replace no_data_value in topofile +0 =: test_topography # (Type topography specification) +1 =: ntopofiles + +'/b.tt3' # topo_path + 3 # topo_type +-98.25 -92.75 22.5 27.5 # crop_extent [x1 x2 y1 y2] +2 # coarsen +3 # buffer +-100.0 20.0 # align [x y] +0.125 # x_shift +-0.25 # y_shift +1.5 # z_shift +T # negate_z diff --git a/tests/data/topo_data_golden/ascii_single.txt b/tests/data/topo_data_golden/ascii_single.txt new file mode 100644 index 000000000..fe1ccc91b --- /dev/null +++ b/tests/data/topo_data_golden/ascii_single.txt @@ -0,0 +1,20 @@ +######################################################## +### DO NOT EDIT THIS FILE: GENERATED AUTOMATICALLY #### +### To modify data, edit setrun.py #### +### and then "make .data" #### +######################################################## + +99999.0 =: topo_missing # replace no_data_value in topofile +0 =: test_topography # (Type topography specification) +1 =: ntopofiles + +'/a.tt3' # topo_path + 3 # topo_type +0. 0. 0. 0. # crop_extent [x1 x2 y1 y2] +1 # coarsen +0 # buffer +0. 0. # align [x y] +0.0 # x_shift +0.0 # y_shift +0.0 # z_shift +F # negate_z diff --git a/tests/data/topo_data_golden/netcdf_cropped.txt b/tests/data/topo_data_golden/netcdf_cropped.txt new file mode 100644 index 000000000..667528f66 --- /dev/null +++ b/tests/data/topo_data_golden/netcdf_cropped.txt @@ -0,0 +1,30 @@ +######################################################## +### DO NOT EDIT THIS FILE: GENERATED AUTOMATICALLY #### +### To modify data, edit setrun.py #### +### and then "make .data" #### +######################################################## + +99999.0 =: topo_missing # replace no_data_value in topofile +0 =: test_topography # (Type topography specification) +1 =: ntopofiles + +'/d.nc' # topo_path + 4 # topo_type +-98.0 -92.0 22.0 28.0 # crop_extent [x1 x2 y1 y2] +2 # coarsen +2 # buffer +0. 0. # align [x y] +0.0 # x_shift +0.0 # y_shift +0.0 # z_shift +F # negate_z +var_name = elevation +x_name = lon +y_name = lat +lon_wrap_offset = 0.0 +y_increasing = True +dim_order = y,x +scale_factor = 1.0 +fill_action = abort +crop_bounds = -98.0 -92.0 22.0 28.0 + diff --git a/tests/data/topo_data_golden/netcdf_no_crop.txt b/tests/data/topo_data_golden/netcdf_no_crop.txt new file mode 100644 index 000000000..1ad2f8348 --- /dev/null +++ b/tests/data/topo_data_golden/netcdf_no_crop.txt @@ -0,0 +1,29 @@ +######################################################## +### DO NOT EDIT THIS FILE: GENERATED AUTOMATICALLY #### +### To modify data, edit setrun.py #### +### and then "make .data" #### +######################################################## + +99999.0 =: topo_missing # replace no_data_value in topofile +0 =: test_topography # (Type topography specification) +1 =: ntopofiles + +'/c.nc' # topo_path + 4 # topo_type +0. 0. 0. 0. # crop_extent [x1 x2 y1 y2] +1 # coarsen +0 # buffer +0. 0. # align [x y] +0.0 # x_shift +0.0 # y_shift +0.0 # z_shift +F # negate_z +var_name = elevation +x_name = lon +y_name = lat +lon_wrap_offset = 0.0 +y_increasing = True +dim_order = y,x +scale_factor = 1.0 +fill_action = abort + diff --git a/tests/data/topo_data_golden/two_files_priority.txt b/tests/data/topo_data_golden/two_files_priority.txt new file mode 100644 index 000000000..3409a01f4 --- /dev/null +++ b/tests/data/topo_data_golden/two_files_priority.txt @@ -0,0 +1,31 @@ +######################################################## +### DO NOT EDIT THIS FILE: GENERATED AUTOMATICALLY #### +### To modify data, edit setrun.py #### +### and then "make .data" #### +######################################################## + +99999.0 =: topo_missing # replace no_data_value in topofile +0 =: test_topography # (Type topography specification) +2 =: ntopofiles + +'/coarse.tt3' # topo_path + 3 # topo_type +0. 0. 0. 0. # crop_extent [x1 x2 y1 y2] +1 # coarsen +0 # buffer +0. 0. # align [x y] +0.0 # x_shift +0.0 # y_shift +0.0 # z_shift +F # negate_z + +'/fine.tt3' # topo_path + 3 # topo_type +0. 0. 0. 0. # crop_extent [x1 x2 y1 y2] +1 # coarsen +0 # buffer +0. 0. # align [x y] +0.0 # x_shift +0.0 # y_shift +0.0 # z_shift +F # negate_z diff --git a/tests/regression/topo_crop/test_topo_descriptor_crop.py b/tests/regression/topo_crop/test_topo_descriptor_crop.py index 6170ba63c..385da3fc1 100644 --- a/tests/regression/topo_crop/test_topo_descriptor_crop.py +++ b/tests/regression/topo_crop/test_topo_descriptor_crop.py @@ -280,3 +280,137 @@ def test_descriptor_crop_inside_window_runs_either_way( xll, xhi, yll, yhi, _mx, _my = _topo_window(tmp_path) assert (xll, xhi, yll, yhi) == pytest.approx( (CROP[0], CROP[1], CROP[2], CROP[3])) + + +# =========================================================================== +# Antimeridian: two entries for one file, with different lon_wrap_offset +# +# This is the end-to-end wrap. Nothing had ever fed GeoClaw two descriptor +# entries pointing at the same file with different lon_wrap_offset values -- +# the Python side was tested, and the Fortran side was assumed. +# =========================================================================== + +# A global file so a crop can genuinely run off its longitude extent. +WRAP_TOPO_X = (-180.0, 180.0) +WRAP_TOPO_Y = (-40.0, 0.0) +WRAP_DELTA = 0.25 + +# Continuous-coordinate crop spanning the seam: covered by 170..180 (offset +# -360) plus -180..-150 (offset 0). +WRAP_CROP = (-190.0, -150.0, -30.0, -10.0) + + +def _write_global_nc(path): + """A CF-compliant global NetCDF topo file spanning the antimeridian. + + Depth varies with longitude so the two sides of the seam carry visibly + different values; a constant field would hide a mis-stitched join. + """ + netCDF4 = pytest.importorskip("netCDF4") + + x = np.arange(WRAP_TOPO_X[0], WRAP_TOPO_X[1] + 1e-9, WRAP_DELTA) + y = np.arange(WRAP_TOPO_Y[0], WRAP_TOPO_Y[1] + 1e-9, WRAP_DELTA) + Z = -3000.0 + 10.0 * np.cos(np.radians(x))[None, :] * np.ones((y.size, 1)) + + with netCDF4.Dataset(path, 'w') as ds: + ds.createDimension('lon', x.size) + ds.createDimension('lat', y.size) + v = ds.createVariable('lon', 'f8', ('lon',)) + v[:] = x + v.units = 'degrees_east' + v.standard_name = 'longitude' + v = ds.createVariable('lat', 'f8', ('lat',)) + v[:] = y + v.units = 'degrees_north' + v.standard_name = 'latitude' + v = ds.createVariable('elevation', 'f8', ('lat', 'lon')) + v[:] = Z + v.units = 'm' + v.standard_name = 'height_above_mean_sea_level' + v.positive = 'up' + ds.Conventions = 'CF-1.8' + return path + + +def _run_wrapped(tmp_path, prebuilt, domain, buffer=0): + """Set up and run a case whose topo.data has two wrapped entries.""" + nc_path = tmp_path / "global.nc" + _write_global_nc(nc_path) + + runner = gtest.GeoClawTestRunner(tmp_path, test_path=testdir) + runner.set_data() + cd = runner.rundata.clawdata + cd.lower = [domain[0], domain[2]] + cd.upper = [domain[1], domain[3]] + cd.num_cells = [20, 12] + runner.write_data() + + # Written the way a user would: one Topography with a continuous + # cross-seam crop. The writer expands it into two descriptor entries. + topo = topotools.Topography() + topo.path = str(nc_path) + topo.topo_type = 4 + topo.crop_extent = list(WRAP_CROP) + topo.buffer = buffer + + td = TopographyData() + td.topofiles = [topo] + td.write(out_file=str(tmp_path / "topo.data")) + + text = (tmp_path / "topo.data").read_text() + assert text.count("lon_wrap_offset") == 2, ( + f"expected two wrapped entries, got:\n{text}") + + shutil.copy(prebuilt, tmp_path / runner.executable_name) + (tmp_path / "_output").mkdir(exist_ok=True) + proc = subprocess.run([str(tmp_path / runner.executable_name)], + cwd=tmp_path, capture_output=True, text=True) + (tmp_path / "run.log").write_text(proc.stdout + proc.stderr) + return proc.returncode, proc.stdout + proc.stderr + + +@pytest.mark.netcdf +def test_wrapped_entries_cover_a_cross_seam_domain(tmp_path, netcdf_xgeoclaw): + """A domain straddling the antimeridian is covered by the two entries. + + The domain sits at -185..-155 in continuous coordinates, i.e. it spans the + date line. Neither entry covers it alone: the offset-0 entry supplies + -180..-150 and the offset -360 entry supplies -190..-180 by reading the + file's 170..180 strip. If Fortran ignored lon_wrap_offset, or applied it + before selecting on crop_bounds, the west half would be missing and the + coverage check would fail. + """ + domain = (-185.0, -155.0, -28.0, -12.0) + rc, out = _run_wrapped(tmp_path, netcdf_xgeoclaw, domain) + + assert "topo arrays do not cover domain" not in out, ( + f"the wrapped pair did not cover a cross-seam domain.\n{out}") + assert rc == 0, out + + +@pytest.mark.netcdf +def test_wrapped_entries_report_shifted_extents(tmp_path, netcdf_xgeoclaw): + """Both entries must be reported in *domain* coordinates. + + fort.geo records each topo grid after lon_wrap_offset is applied, so the + wrapped entry has to appear west of -180 rather than at its file position + of +170. Seeing +170 here would mean the offset never reached the grid. + """ + domain = (-185.0, -155.0, -28.0, -12.0) + _run_wrapped(tmp_path, netcdf_xgeoclaw, domain) + + text = (tmp_path / "fort.geo").read_text() + windows = re.findall( + r"mx\s*=\s*\d+\s+x\s*=\s*\(\s*([-\d.E+]+)\s*,\s*([-\d.E+]+)\s*\)", + text) + assert len(windows) == 2, f"expected two topo grids in fort.geo:\n{text}" + + lows = sorted(float(lo) for lo, _hi in windows) + highs = sorted(float(hi) for _lo, hi in windows) + # Wrapped entry shifted to -190..-180; unwrapped entry at -180..-150. + assert lows[0] == pytest.approx(-190.0) + assert highs[-1] == pytest.approx(-150.0) + # Nothing left sitting at its unshifted file position. + assert all(float(hi) < 0.0 for _lo, hi in windows), ( + f"an entry kept its file longitude instead of the domain one: " + f"{windows}") diff --git a/tests/test_topo_data_golden.py b/tests/test_topo_data_golden.py new file mode 100644 index 000000000..865dcab5f --- /dev/null +++ b/tests/test_topo_data_golden.py @@ -0,0 +1,204 @@ +#!/usr/bin/env python +# encoding: utf-8 +"""Byte-exact goldens for the ``topo.data`` writer. + +Every other test of :meth:`TopographyData.write` checks substrings or line +counts, so the two write paths could diverge in field width, whitespace, +quoting or float formatting and still pass. Fortran reads this file with +fixed-format ``read`` statements, so those details are the contract. + +These goldens were generated **before** ``write()`` was restructured into +resolve-then-write (PR A2). That restructuring must be a pure refactor for +every case that does not wrap across the antimeridian, and this file is what +pins it: if a byte moves in any case below, the refactor changed behaviour it +was not supposed to touch. + +Regenerate deliberately with ``GEOCLAW_REGEN=1`` and review the diff -- never +because a test failed. + +The absolute path GeoClaw writes into ``topo.data`` depends on ``tmp_path``, so +paths are replaced with ```` before comparison. Everything else, +including the exact spacing around ``# topo_type`` and the ``repr`` of every +float, is compared verbatim. +""" + +import os +import re +from pathlib import Path + +import numpy as np +import pytest + +import clawpack.geoclaw.topotools as topotools +from clawpack.geoclaw.data import TopographyData + +testdir = Path(__file__).parent +golden_dir = testdir / "data" / "topo_data_golden" + +pytestmark = pytest.mark.python + + +def _write_tt3(path, x0=-100.0, x1=-90.0, y0=20.0, y1=30.0, delta=0.5): + """A small ASCII topo_type=3 file on an exactly-representable lattice.""" + x = np.arange(x0, x1 + 1e-9, delta) + y = np.arange(y0, y1 + 1e-9, delta) + Z = -100.0 + np.outer(np.linspace(0.0, 50.0, y.size), np.ones_like(x)) + topo = topotools.Topography() + topo.set_xyZ(x, y, Z) + topo.write(str(path), topo_type=3) + return path + + +def _write_nc(path, x0=-100.0, x1=-90.0, y0=20.0, y1=30.0, delta=0.5): + """A CF-compliant NetCDF topo file (topo_type=4).""" + netCDF4 = pytest.importorskip("netCDF4") + x = np.arange(x0, x1 + 1e-9, delta) + y = np.arange(y0, y1 + 1e-9, delta) + Z = -100.0 + np.outer(np.linspace(0.0, 50.0, y.size), np.ones_like(x)) + with netCDF4.Dataset(path, "w") as ds: + ds.createDimension("lon", x.size) + ds.createDimension("lat", y.size) + v = ds.createVariable("lon", "f8", ("lon",)) + v[:] = x + v.units = "degrees_east" + v.standard_name = "longitude" + v = ds.createVariable("lat", "f8", ("lat",)) + v[:] = y + v.units = "degrees_north" + v.standard_name = "latitude" + v = ds.createVariable("elevation", "f8", ("lat", "lon")) + v[:] = Z + v.units = "m" + v.standard_name = "height_above_mean_sea_level" + v.positive = "up" + ds.Conventions = "CF-1.8" + return path + + +def _normalize(text, tmp_path): + """Replace machine-specific absolute paths with a stable placeholder.""" + text = text.replace(str(Path(tmp_path).resolve()), "") + text = text.replace(str(tmp_path), "") + # The data_source header line carries a timestamp/path in some setups. + return re.sub(r"[^\s']*/", "/", text) + + +def _check(name, tmp_path, out_file): + # ".txt", not ".data": the repo's .gitignore excludes "*.data", so a + # golden named topo.data would be silently left out of the commit and the + # test would fail for everyone else. Matches the met_forcing goldens. + golden = golden_dir / f"{name}.txt" + actual = _normalize(Path(out_file).read_text(), tmp_path) + + if os.environ.get("GEOCLAW_REGEN"): + golden.parent.mkdir(parents=True, exist_ok=True) + golden.write_text(actual) + return + + assert golden.exists(), ( + f"Missing golden {golden}. Generate it with GEOCLAW_REGEN=1 from the " + f"code you intend to pin -- not from code you are about to change.") + expected = golden.read_text() + assert actual == expected, ( + f"topo.data for '{name}' changed.\n--- expected ---\n{expected}\n" + f"--- actual ---\n{actual}") + + +def test_golden_ascii_single(tmp_path): + """One ASCII file, no preprocessing: the simplest possible topo.data.""" + path = _write_tt3(tmp_path / "a.tt3") + topo = topotools.Topography() + topo.path = str(path) + topo.topo_type = 3 + + td = TopographyData() + td.topofiles = [topo] + out = tmp_path / "topo.data" + td.write(out_file=str(out)) + _check("ascii_single", tmp_path, out) + + +def test_golden_ascii_all_preprocessing(tmp_path): + """Every preprocessing attribute non-default at once. + + Pins the float formatting of the crop_extent and align lines: these are + written with repr() so coordinates reach Fortran at full precision, and a + change to %g would silently truncate them. + """ + path = _write_tt3(tmp_path / "b.tt3") + topo = topotools.Topography() + topo.path = str(path) + topo.topo_type = 3 + topo.crop_extent = [-98.25, -92.75, 22.5, 27.5] + topo.coarsen = 2 + topo.buffer = 3 + topo.align = [-100.0, 20.0] + topo.x_shift = 0.125 + topo.y_shift = -0.25 + topo.z_shift = 1.5 + topo.negate_z = True + + td = TopographyData() + td.topofiles = [topo] + out = tmp_path / "topo.data" + td.write(out_file=str(out)) + _check("ascii_all_preprocessing", tmp_path, out) + + +def test_golden_two_files_priority_order(tmp_path): + """Two files of different resolution: pins the coarsest-first ordering.""" + coarse = _write_tt3(tmp_path / "coarse.tt3", delta=1.0) + fine = _write_tt3(tmp_path / "fine.tt3", delta=0.25) + + topos = [] + for p in (fine, coarse): # deliberately not in priority order + t = topotools.Topography() + t.path = str(p) + t.topo_type = 3 + topos.append(t) + + td = TopographyData() + td.topofiles = topos + out = tmp_path / "topo.data" + td.write(out_file=str(out)) + _check("two_files_priority", tmp_path, out) + + +@pytest.mark.netcdf +def test_golden_netcdf_no_crop(tmp_path): + """type-4 with no crop: pins the whole CF descriptor block.""" + pytest.importorskip("xarray") + path = _write_nc(tmp_path / "c.nc") + topo = topotools.Topography() + topo.path = str(path) + topo.topo_type = 4 + + td = TopographyData() + td.topofiles = [topo] + out = tmp_path / "topo.data" + td.write(out_file=str(out)) + _check("netcdf_no_crop", tmp_path, out) + + +@pytest.mark.netcdf +def test_golden_netcdf_cropped(tmp_path): + """type-4 with an ordinary in-extent crop and a buffer. + + This is the case the resolve-then-write refactor touches most, and the one + most likely to change by accident: it must keep emitting exactly one entry, + with crop_bounds in file coordinates and lon_wrap_offset 0.0. + """ + pytest.importorskip("xarray") + path = _write_nc(tmp_path / "d.nc") + topo = topotools.Topography() + topo.path = str(path) + topo.topo_type = 4 + topo.crop_extent = [-98.0, -92.0, 22.0, 28.0] + topo.buffer = 2 + topo.coarsen = 2 + + td = TopographyData() + td.topofiles = [topo] + out = tmp_path / "topo.data" + td.write(out_file=str(out)) + _check("netcdf_cropped", tmp_path, out) diff --git a/tests/test_topotools_preprocessing.py b/tests/test_topotools_preprocessing.py index 50f2d55d5..cf8c1bf64 100644 --- a/tests/test_topotools_preprocessing.py +++ b/tests/test_topotools_preprocessing.py @@ -35,6 +35,7 @@ from __future__ import annotations +import re import textwrap import warnings from pathlib import Path @@ -1433,7 +1434,13 @@ def test_unstructured_with_crop_raises_not_typeerror(tmp_path): def test_cross_seam_crop_raises_at_topo_data_write(global_tt3_path, tmp_path): """Writing a descending crop_extent used to emit `crop_bounds = 170.0 - -170.0`, which Fortran resolves to mx=0, my=0: an empty topo, no error.""" + -170.0`, which Fortran resolves to mx=0, my=0: an empty topo, no error. + + The wrapped spelling stays an error even though the *continuous* spelling + is now split across the seam automatically -- a descending pair has no + unambiguous reading. The message must offer the continuous equivalent + rather than telling the caller to go build descriptors by hand. + """ t = Topography() t.path = str(global_tt3_path) t.topo_type = 3 @@ -1441,9 +1448,12 @@ def test_cross_seam_crop_raises_at_topo_data_write(global_tt3_path, tmp_path): td = TopographyData() td.topofiles = [t] - with pytest.raises(ValueError, match="crosses the antimeridian"): + with pytest.raises(ValueError, match="descending in longitude") as excinfo: td.write(out_file=str(tmp_path / "topo.data")) + # _WRAPPED_CROP is [170, -170]; the continuous equivalent is [-190, -170]. + assert "[-190.0, -170.0]" in str(excinfo.value) + def test_topo_type_none_inferred_from_suffix(global_tt3_path, tmp_path): """topo_type=None reached the `:3d` format and raised a TypeError naming @@ -1529,3 +1539,252 @@ def test_buffer_and_coarsen_give_absolute_output_shape(tt2_path): assert float(t.x[1] - t.x[0]) == pytest.approx(2.0 * _DELTA) # And the buffered window still contains the unbuffered one. assert float(t.x[0]) == pytest.approx(float(ref.x[0]) - 2.0 * _DELTA) + + +# =========================================================================== +# Group N+1 — Remote sources, and cross-seam crops from ordinary setrun code +# +# Both of these were reported from the field: the most natural possible setrun +# (set .path, .crop_extent, .buffer; append to topofiles) failed, once with a +# mangled local path and once with "crop_bounds exceed file extent", while the +# machinery to handle each already existed and was tested one layer down. +# These pin the two paths being reachable, not just present. +# =========================================================================== + +REMOTE_URL = ("https://www.ngdc.noaa.gov/thredds/dodsC/global/ETOPO2022/30s/" + "30s_bed_elev_netcdf/ETOPO_2022_v1_30s_N90W180_bed.nc") + + +def _make_global_nc(path, delta=0.5, lon0=-180.0, lon1=180.0, + lat0=-70.0, lat1=10.0): + """A CF-compliant near-global NetCDF file spanning the antimeridian.""" + netCDF4 = pytest.importorskip("netCDF4") + x = np.arange(lon0, lon1 + 1e-9, delta) + y = np.arange(lat0, lat1 + 1e-9, delta) + Z = -1000.0 + np.outer(np.linspace(0.0, 200.0, y.size), np.ones_like(x)) + with netCDF4.Dataset(path, "w") as ds: + ds.createDimension("lon", x.size) + ds.createDimension("lat", y.size) + v = ds.createVariable("lon", "f8", ("lon",)) + v[:] = x + v.units = "degrees_east" + v.standard_name = "longitude" + v = ds.createVariable("lat", "f8", ("lat",)) + v[:] = y + v.units = "degrees_north" + v.standard_name = "latitude" + v = ds.createVariable("elevation", "f8", ("lat", "lon")) + v[:] = Z + v.units = "m" + v.standard_name = "height_above_mean_sea_level" + v.positive = "up" + ds.Conventions = "CF-1.8" + return path + + +def _entry_blocks(text): + """Split a written topo.data into its per-file blocks.""" + return [b for b in text.split("# topo_path") if "topo_type" in b] + + +def _descriptor_values(text, key): + """Every value written for descriptor *key*, in file order.""" + return re.findall(rf"^{re.escape(key)}\s*=\s*(.+)$", text, re.MULTILINE) + + +def test_is_remote_url_discriminates_urls_from_paths(): + """The regex must not mistake a Windows drive letter for a URL scheme.""" + from clawpack.geoclaw.netcdf_utils import is_remote_url + + assert is_remote_url("https://example.org/topo.nc") + assert is_remote_url("http://example.org/topo.nc") + assert not is_remote_url("/tmp/topo.nc") + assert not is_remote_url("topo.nc") + assert not is_remote_url(r"C:\data\topo.nc") + assert not is_remote_url(Path("/tmp/topo.nc")) + + +def test_remote_url_in_topofiles_raises_with_the_recipe(tmp_path): + """A URL used to be run through os.path.abspath, producing + + FileNotFoundError: /run/dir/https:/www.ngdc.noaa.gov/... + + naming a path the user never typed and giving no hint that the fix is to + fetch it first. + """ + t = Topography() + t.path = REMOTE_URL + t.topo_type = 4 + t.crop_extent = [-160.0, -120.0, -60.0, 0.0] + + td = TopographyData() + td.topofiles = [t] + with pytest.raises(ValueError) as excinfo: + td.write(out_file=str(tmp_path / "topo.data")) + + msg = str(excinfo.value) + assert "fetch_remote_topo" in msg # the actionable part + assert REMOTE_URL in msg # unmangled + assert "https:/www" not in msg # specifically not collapsed + + +def test_remote_url_in_dtopofiles_raises(tmp_path): + """Same trap on the dtopo writer, which shares the abspath pattern.""" + import clawpack.geoclaw.dtopotools as dtopotools + + d = dtopotools.DTopography() + d.path = REMOTE_URL + d.dtopo_type = 4 + + from clawpack.geoclaw.data import DTopoData + + dtd = DTopoData() + dtd.dtopofiles = [d] + with pytest.raises(ValueError, match="URL"): + dtd.write(out_file=str(tmp_path / "dtopo.data")) + + +@pytest.mark.netcdf +def test_cross_seam_crop_writes_two_entries(tmp_path): + """The reported case: a continuous crop spanning the date line. + + Previously raised `crop_bounds lon [...] exceed file extent` even though + _compute_lon_entries could already cover it. Must now produce one entry + per side of the seam with complementary crop_bounds. + """ + pytest.importorskip("xarray") + nc = _make_global_nc(tmp_path / "gebco_like.nc") + + t = Topography() + t.path = str(nc) + t.topo_type = 4 + t.crop_extent = [-190.0, -120.0, -60.0, 0.0] + + td = TopographyData() + td.topofiles = [t] + out = tmp_path / "topo.data" + td.write(out_file=str(out)) + text = out.read_text() + + assert "2 =: ntopofiles" in text + assert len(_entry_blocks(text)) == 2 + + offsets = [float(v) for v in _descriptor_values(text, "lon_wrap_offset")] + assert offsets == [0.0, -360.0] + + bounds = _descriptor_values(text, "crop_bounds") + # East side comes from the file as-is; west side is the +170..180 strip + # read with a -360 shift so Fortran places it at -190..-180. + assert bounds[0].split() == ["-180.0", "-120.0", "-60.0", "0.0"] + assert bounds[1].split() == ["170.0", "180.0", "-60.0", "0.0"] + + +@pytest.mark.netcdf +def test_cross_seam_entries_keep_buffer_and_coarsen(tmp_path): + """buffer and coarsen must reach *every* expanded entry. + + topo_entries() builds Topography objects carrying only _netcdf_meta, so + routing through it naively writes `buffer = 0` -- which would silently + undo the Fortran fix that made buffer work for descriptor crops at all. + """ + pytest.importorskip("xarray") + nc = _make_global_nc(tmp_path / "gebco_like.nc") + + t = Topography() + t.path = str(nc) + t.topo_type = 4 + t.crop_extent = [-190.0, -120.0, -60.0, 0.0] + t.buffer = 1 + t.coarsen = 20 + + td = TopographyData() + td.topofiles = [t] + out = tmp_path / "topo.data" + td.write(out_file=str(out)) + + blocks = _entry_blocks(out.read_text()) + assert len(blocks) == 2 + for block in blocks: + assert "1 # buffer" in block + assert "20 # coarsen" in block + + +@pytest.mark.netcdf +def test_off_seam_crop_writes_one_entry_with_nonzero_offset(tmp_path): + """A crop wholly on the far side of the cut needs a *single* entry with a + non-zero offset. lon_wrap_offset was hard-coded to 0.0, so this case was + wrong even though it never needed splitting.""" + pytest.importorskip("xarray") + nc = _make_global_nc(tmp_path / "g.nc") + + t = Topography() + t.path = str(nc) + t.topo_type = 4 + t.crop_extent = [185.0, 195.0, -60.0, 0.0] # i.e. -175..-165 + + td = TopographyData() + td.topofiles = [t] + out = tmp_path / "topo.data" + td.write(out_file=str(out)) + text = out.read_text() + + assert "1 =: ntopofiles" in text + assert [float(v) for v in + _descriptor_values(text, "lon_wrap_offset")] == [360.0] + assert _descriptor_values(text, "crop_bounds")[0].split() == [ + "-175.0", "-165.0", "-60.0", "0.0"] + + +@pytest.mark.netcdf +def test_wrapping_write_does_not_scan_the_whole_file(tmp_path, monkeypatch): + """topo_entries() inspects with crop_bounds unset, so its fill/magnitude + checks would read the *entire* variable and reject NaN anywhere in it. + + On a global DEM read over OPeNDAP that turns `make data` into a full + download. The regression is invisible on a small local fixture, so it is + pinned directly rather than by timing. + """ + pytest.importorskip("xarray") + from clawpack.geoclaw import netcdf_utils as ncutils + + nc = _make_global_nc(tmp_path / "g.nc") + + calls = [] + original = ncutils.TopoInspector._check_fill_in_crop + + def spy(self, *args, **kwargs): + calls.append(args) + return original(self, *args, **kwargs) + + monkeypatch.setattr(ncutils.TopoInspector, "_check_fill_in_crop", spy) + + t = Topography() + t.path = str(nc) + t.topo_type = 4 + t.crop_extent = [-190.0, -120.0, -60.0, 0.0] + + td = TopographyData() + td.topofiles = [t] + td.write(out_file=str(tmp_path / "topo.data")) + + assert calls == [], ( + "write() triggered a fill scan; on a remote global DEM this " + "downloads the whole file during `make data`.") + + +@pytest.mark.netcdf +def test_wrapping_crop_still_checks_latitude(tmp_path): + """Longitude wraps; latitude does not. Dropping crop_bounds validation to + allow the wrap must not also drop the latitude check.""" + pytest.importorskip("xarray") + nc = _make_global_nc(tmp_path / "g.nc") # lat spans -70..10 + + t = Topography() + t.path = str(nc) + t.topo_type = 4 + t.crop_extent = [-190.0, -120.0, -60.0, 45.0] # 45N is off the file + + td = TopographyData() + td.topofiles = [t] + with pytest.raises(ValueError, match="latitude extent"): + td.write(out_file=str(tmp_path / "topo.data")) From cce6a5197b4548cf4359518ab66fc4956f1d4238 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Sat, 5 Sep 2026 21:32:22 -0400 Subject: [PATCH 3/7] State the units policy, enforce it, and fix what violated it GeoClaw's units policy was implemented in the NetCDF readers but written down nowhere, so enforcement drifted. Adds dev/design/units_policy.md, a UNITS_POLICY registry that both the doc table and a conformance test are generated from, and fixes four violations: CSVFault.read parsed unit annotations from column headings and discarded them (alaska1964.csv read as Mw 5.20 instead of 8.53), input_units was a mutated mutable default that silently declared SI, a unit-less dtopo time axis was silently assumed to be seconds, and two docstrings claimed GeoClaw does not convert on read while it does. Non-conforming rows are xfail(strict) so the remaining ASCII gaps stay visible. Signed-off-by: Kyle Mandli Assisted-by: claude claude-opus-5[1m] --- dev/design/units_policy.md | 88 +++++++ src/python/geoclaw/dtopotools.py | 109 +++++++-- src/python/geoclaw/netcdf_utils.py | 15 ++ src/python/geoclaw/topotools.py | 22 +- src/python/geoclaw/units.py | 167 +++++++++++++ tests/test_units_policy.py | 366 +++++++++++++++++++++++++++++ 6 files changed, 744 insertions(+), 23 deletions(-) create mode 100644 dev/design/units_policy.md create mode 100644 tests/test_units_policy.py diff --git a/dev/design/units_policy.md b/dev/design/units_policy.md new file mode 100644 index 000000000..5c03c1223 --- /dev/null +++ b/dev/design/units_policy.md @@ -0,0 +1,88 @@ +# GeoClaw units policy + +**Status:** stated and enforced as of PR B1. The table below is generated from +`UNITS_POLICY` in `src/python/geoclaw/units.py` and verified against the real +readers by `tests/test_units_policy.py`. + +## Why this document exists + +GeoClaw already had a coherent units policy. It was implemented in the NetCDF +readers, it was well reasoned, and it was written down nowhere so enforcement +drifted, the override argument acquired four different spellings, and two +docstrings ended up claiming the opposite of what the code does. Users could not +discover any of it, and neither could reviewers. + +Writing the rules down is only half of it. Prose drifts from code, and a table +in a document is exactly the kind of thing that quietly stops being true. So the +table here is *generated* from a registry, and every row is *executed* against +the real reader. A row cannot claim behaviour the code does not have, and the +document cannot disagree with the registry. + +## The rules + +1. **Units must be declared in the file.** They are never silently assumed. A + file with no declaration is an error, not a guess. +2. **Overriding is explicit.** `assume_units` (NetCDF) and `input_units` + (subfault files) mean "treat the file as if it had declared this". They are + the only way to supply units GeoClaw cannot read from the file. +3. **A recognised non-contract unit is converted, and the conversion is + announced.** Reading a file in km is fine; doing it silently is not. +4. **An unrecognised unit raises.** GeoClaw does not guess at unit strings it + does not know. +5. **After conversion, magnitude is sanity-checked.** Units can be declared + *wrongly*, and rules 1-4 cannot catch that. `_check_magnitude` does, within + limits: it auto-corrects only the one unambiguous case (pressure ~1000x low + is unmistakably hPa/mbar) and raises on everything else, because feet vs + meters and knots vs m/s cannot be told apart by magnitude alone. + +Rule 5 is the reason the policy is not simply "trust the declaration". Rules 1-4 +protect against *missing* information; rule 5 protects against *wrong* +information, which is the more common failure in practice. + +## Contract units + +GeoClaw works internally in SI: elevation and deformation in meters, wind in +m/s, pressure in Pa, time in seconds. See `GEOCLAW_NETCDF_UNITS` in +`src/python/geoclaw/units.py`. + +## What each input path actually does + +Regenerate this table with `GEOCLAW_REGEN=1 pytest tests/test_units_policy.py`. +Do not edit it by hand -- edit `UNITS_POLICY` instead. + + +| Path | Contract | Declared in file | Override | Missing | Non-contract | Unrecognised | Magnitude | Conforms | +|---|---|---|---|---|---|---|---|---| +| `Topography.read (topo_type=4)` | m | yes | nc_params={'assume_units': str} | raise | convert+warn | raise | yes | yes | +| `Topography.read (topo_type=1,2,3)` | m | no | none | silent-assume | n/a | n/a | no | **no** -- ASCII carries no units and has no override; elevation in cm or feet is read as metres with no message and no sanity check. | +| `DTopoInspector (deformation)` | m | yes | assume_units (str) | raise | convert+warn | raise | no | yes | +| `DTopoInspector (time axis)` | s | yes | none | warn+assume | convert | raise | no | yes | +| `DTopography.read (dtopo_type=1,2,3)` | m | no | none | silent-assume | n/a | n/a | no | **no** -- ASCII dtopo carries no units and has no override. | +| `MetInspector (wind, pressure)` | m/s, Pa | yes | assume_units (bool), format_units (dict) | raise | convert+warn | raise | yes | yes | +| `Fault.read (columnar subfault files)` | m, Pa | no | input_units (dict) | warn+assume | convert | raise | no | yes | +| `CSVFault.read (units in column headings)` | m, Pa | yes | input_units (dict), overrides the heading | warn+assume | convert | raise | no | yes | + + +A row marked **no** in the *Conforms* column is a known hole: the policy says +one thing and that reader does another. Its conformance test is marked +`xfail(strict=True)`, so it is visible in every test run and the marker cannot +be left behind once the hole is closed. + +## Known deliberate looseness + +- **Unit-string matching is a hand-maintained table**, not dimensional + analysis (see `_CONTRACT_UNIT_CF_ALIASES` and `_CF_TO_UNITS_PY` in + `netcdf_utils.py`). It is exact-string and case-sensitive, so `meters` is + accepted and `Meters` is not, and `feet` has no conversion at all. Using + `cf-units`/udunits instead would be correct but adds a dependency. +- **`MetInspector.assume_units` is a `bool` where the other inspectors take a + `str`.** This is deliberate, not an oversight: Met covers several variables + with *different* contract units, so a single string is meaningless. `True` + means "each variable is already in its own contract unit"; the per-role string + form is `format_units={role: unit}`, which takes precedence. + +## Definition of done + +Any change to unit handling must update `UNITS_POLICY` in the same PR. The +generated table and the conformance test will fail otherwise, which is the +point. diff --git a/src/python/geoclaw/dtopotools.py b/src/python/geoclaw/dtopotools.py index cf451a1cd..7ffa9ccd0 100644 --- a/src/python/geoclaw/dtopotools.py +++ b/src/python/geoclaw/dtopotools.py @@ -33,6 +33,7 @@ import os import sys import re +import warnings import numpy @@ -65,6 +66,58 @@ standard_units['mu'] = 'Pa' + +def _resolve_input_units(input_units, where, from_file=None): + """Return the effective ``{parameter: unit}`` mapping for a subfault read. + + Subfault files are columnar text: most carry no unit information at all, + and the few that do (CSV headings like ``Depth(km)``) describe only some + columns. So units come from up to three places, and the precedence matters: + + 1. an explicit *input_units* entry from the caller -- always wins; + 2. *from_file*, units parsed out of the file itself -- fills the gaps; + 3. :data:`standard_units` (SI) -- the last resort. + + Explicit beats implicit so that a caller who knows better than a file's + heading can say so, and so a wrong heading is recoverable. A disagreement + between (1) and (2) is warned about rather than resolved quietly. + + *input_units* of None means "not specified": SI is assumed, but a warning + says so, because omitting it used to declare metres/pascals silently and a + km / dyne-cm file was then off by 10^3-10^7. An explicit ``{}`` means "my + data really is SI" and stays silent -- the deliberate escape hatch. + + The caller's dict is never mutated; a copy is returned. + """ + resolved = standard_units.copy() + + if from_file: + resolved.update(from_file) + + if input_units is None: + warnings.warn( + f"No input_units given for {where}; assuming GeoClaw standard " + f"units ({', '.join('%s=%s' % kv for kv in sorted(standard_units.items()))}). " + f"If the file uses other units (km, cm, dyne/cm^2 are common) pass " + f"input_units explicitly; pass input_units={{}} to state that the " + f"data really is in standard units and silence this warning.", + UserWarning, stacklevel=3, + ) + return resolved + + for name, unit in input_units.items(): + if from_file and name in from_file and from_file[name] != unit: + warnings.warn( + f"Units for '{name}' in {where} disagree: the file says " + f"'{from_file[name]}' and input_units says '{unit}'. Using " + f"'{unit}' (an explicit input_units entry takes precedence " + f"over the file's own heading).", + UserWarning, stacklevel=3, + ) + resolved[name] = unit + + return resolved + def plot_dZ_contours(x, y, dZ, axes=None, dZ_interval=0.5, verbose=False, fig_kwargs={}): r"""For plotting seafloor deformation dZ""" @@ -903,7 +956,7 @@ class Fault(object): """ - def __init__(self, subfaults=None, input_units={}, + def __init__(self, subfaults=None, input_units=None, coordinate_specification=None): r"""Fault initialization routine. @@ -916,9 +969,14 @@ def __init__(self, subfaults=None, input_units={}, #self.times = numpy.array([0., 1.]) # or just [0.] ?? self.dtopo = None - # Default units of each parameter type - self.input_units = standard_units.copy() - self.input_units.update(input_units) + # Units of each parameter type. Only warn about an unspecified + # mapping when there is data to convert; constructing an empty Fault + # and filling it later is a normal pattern and converts nothing. + if subfaults is None and input_units is None: + self.input_units = standard_units.copy() + else: + self.input_units = _resolve_input_units( + input_units, f"{type(self).__name__}()") # Set the coordinate specification, e.g. 'top center': self.coordinate_specification = coordinate_specification @@ -938,7 +996,8 @@ def __init__(self, subfaults=None, input_units={}, def read(self, path, column_map, coordinate_specification="centroid", rupture_type="static", skiprows=0, - delimiter=None, input_units={}, defaults=None): + delimiter=None, input_units=None, + defaults=None, _units_from_file=None): r"""Read in subfault specification at *path*. Creates a list of subfaults from the subfault specification file at @@ -976,8 +1035,8 @@ def read(self, path, column_map, coordinate_specification="centroid", data = numpy.array([data]) self.coordinate_specification = coordinate_specification - self.input_units = standard_units.copy() - self.input_units.update(input_units) + self.input_units = _resolve_input_units( + input_units, f"'{path}'", from_file=_units_from_file) self.subfaults = [] for n in range(data.shape[0]): @@ -3107,13 +3166,17 @@ class CSVFault(Fault): Assumes that the first row gives the column headings """ - def read(self, path, input_units={}, coordinate_specification="top center", + def read(self, path, input_units=None, coordinate_specification="top center", rupture_type="static", verbose=False): r"""Read in subfault specification at *path*. Creates a list of subfaults from the subfault specification file at *path*. + Units may be annotated in the column headings, e.g. ``Depth(km)``. + Those are applied; an explicit *input_units* entry for the same column + overrides them and a disagreement warns. See + ``dev/design/units_policy.md``. """ possible_column_names = """longitude latitude length width depth strike dip @@ -3126,6 +3189,10 @@ def read(self, path, input_units={}, coordinate_specification="top center", param["rupture time"] = "rupture_time" param["rise time"] = "rise_time" + # Units parsed out of the column headings, e.g. "Depth(km)". Keyed by + # the *file's* column name; remapped to parameter names below. + units_from_file = {} + # Read header of file with open(path, 'r') as subfault_file: header_line = subfault_file.readline().split(",") @@ -3135,13 +3202,16 @@ def read(self, path, input_units={}, coordinate_specification="top center", # Strip out units if present unit_start = column_heading.find("(") unit_end = column_heading.find(")") - column_name = column_heading[:unit_start].lower() + column_name = column_heading[:unit_start].lower().strip() units = column_heading[unit_start+1:unit_end] - if verbose and input_units.get(column_name,units) != units: - print("*** Warning: input_units[%s] reset to %s" \ - % (column_name, units)) - print(" based on file header") - input_units[column_name] = units + # Record what the file says about this column. This used + # to be assigned only when `verbose` was true *and* the + # caller had already named a different unit -- so with the + # default verbose=False the heading was parsed and then + # thrown away, and a "Depth(km)" file was read as metres. + # Precedence against input_units is resolved in + # _resolve_input_units, not here. + units_from_file[column_name] = units else: column_name = column_heading.lower() @@ -3153,10 +3223,17 @@ def read(self, path, input_units={}, coordinate_specification="top center", print("*** Warning: column name not recognized: %s" \ % column_name) + # Remap heading names onto parameter names (e.g. "rigidity" -> "mu") + # so they line up with input_units / standard_units keys. + units_from_file = {param.get(name, name): unit + for name, unit in units_from_file.items() + if param.get(name, name) in standard_units} + super(CSVFault, self).read(path, column_map=column_map, skiprows=1, delimiter=",", input_units=input_units, coordinate_specification=coordinate_specification, - rupture_type=rupture_type) + rupture_type=rupture_type, + _units_from_file=units_from_file) @@ -3480,7 +3557,7 @@ class Fault1d(Fault): """ - def __init__(self, subfaults=None, input_units={}, + def __init__(self, subfaults=None, input_units=None, coordinate_specification=None): r"""Fault initialization routine. diff --git a/src/python/geoclaw/netcdf_utils.py b/src/python/geoclaw/netcdf_utils.py index 24c30a194..975d619b6 100644 --- a/src/python/geoclaw/netcdf_utils.py +++ b/src/python/geoclaw/netcdf_utils.py @@ -1185,6 +1185,21 @@ def _compute_time_axis(self, time_name: str) -> tuple[float, float, int]: if cf_unit: seconds = arr * _cf_time_units_to_seconds_factor(cf_unit) else: + # The assumption stays -- a bare numeric dtopo time axis has + # always meant seconds and files rely on it -- but it is stated. + # Silence here is what makes an "hours" file run 3600x too + # fast with nothing to notice, and it is the one place this + # module departs from its own rule (see + # _cf_time_units_to_seconds_factor: callers must never silently + # assume seconds). + warnings.warn( + f"Time axis '{time_name}' in '{self.path}' is numeric with " + f"no 'units' attribute; assuming seconds. Add a CF 'units' " + f"attribute (e.g. 'seconds', 'hours') to the time " + f"coordinate -- a file in hours read as seconds runs " + f"3600x too fast.", + stacklevel=3, + ) seconds = arr if mt < 2: diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py index 3f72d9d00..69e683755 100644 --- a/src/python/geoclaw/topotools.py +++ b/src/python/geoclaw/topotools.py @@ -938,7 +938,7 @@ def generate_2d_coordinates(self, mask=False): def read(self, path=None, topo_type=None, unstructured=False, mask=False, crop_extent=None, force=False, coarsen=None, align=_ALIGN_UNSET, buffer=None, stride=None, - nc_params={}, filter_region=_CROP_EXTENT_UNSET): + nc_params=None, filter_region=_CROP_EXTENT_UNSET): r"""Read in the data from the object's *path* attribute. Stores the resulting data in one of the sets of *x*, *y*, and *z* or @@ -982,16 +982,23 @@ class docstring). Passing it here is equivalent to setting the - `z_var` (str): name of the elevation variable, if it cannot be auto-detected by CF `standard_name` or common names. - `assume_units` (str): unit to assume for the elevation variable - when the file has **no** `units` attribute (e.g. `"m"`). Units - are otherwise required and never silently assumed: a file whose - elevation variable lacks `units`, or whose units are not meters, - raises `ValueError` (GeoClaw does not convert on read; pre- - convert non-meter data to meters first). + when the file has **no** `units` attribute (e.g. `"m"`), treated + as if the file had declared it -- so `assume_units="km"` also + converts. Units are otherwise required and never silently + assumed: a file whose elevation variable lacks `units` raises + `ValueError`. A *recognised* non-meter unit (e.g. `km`) is + converted to meters on read, with a warning; an unrecognised + unit raises. See `dev/design/units_policy.md`. The first three might have already been set when instatiating object. """ + # None is the natural "no options" value and used to reach .get() as a + # NoneType; it is also safer than a shared mutable default. + if nc_params is None: + nc_params = {} + # A crop_extent passed here is equivalent to setting the attribute first; # fold the deprecated filter_region alias onto it, then store it so the # single attribute-driven crop below (and the type-4 pushdown) apply it. @@ -2556,7 +2563,8 @@ def fetch_remote_topo(name_or_url, crop_extent=None, coarsen=1, buffer=0, This is the modern one-call "remote DEM -> Topography" path. It resolves a nickname or URL and reads it through the `topo_type=4` reader (`Topography.read`, backed by `netcdf_utils.TopoInspector`), so it inherits - that path's unit handling (elevation must be in meters, or supply + that path's unit handling (a recognised non-meter unit such as `km` is + converted on read with a warning; a file with no `units` attribute needs `assume_units` via `nc_params`), datum handling, fill->NaN conversion, CF coordinate/variable detection, and lazy hyperslab windowing. diff --git a/src/python/geoclaw/units.py b/src/python/geoclaw/units.py index ed8fbeb04..a14268dc8 100644 --- a/src/python/geoclaw/units.py +++ b/src/python/geoclaw/units.py @@ -10,6 +10,7 @@ import sys import collections +import dataclasses from clawpack.geoclaw.data import LAT2METER @@ -118,6 +119,172 @@ } +# --------------------------------------------------------------------------- +# Units policy registry +# --------------------------------------------------------------------------- +# GeoClaw's units policy is stated in dev/design/units_policy.md. This registry +# is the machine-readable form of the table in that document, and it is the +# single source of truth for both: +# +# * tests/test_units_policy.py, which drives each real reader with a fixture +# and asserts the behaviour declared here -- so a row cannot claim +# something the code does not do; and +# * the generated table in dev/design/units_policy.md, rendered from these +# rows -- so the document cannot drift from the registry. +# +# A row whose `gap` is not None does *not* yet meet the policy. Its +# conformance test is marked xfail(strict=True), so the row is visible as a +# known hole and the marker cannot be left behind once the hole is closed. + + +@dataclasses.dataclass(frozen=True) +class UnitsPolicyRow: + """One input path's declared units behaviour. + + The field values are the vocabulary the conformance test understands; see + ON_MISSING_VALUES etc. below. + """ + + key: str # stable identifier used by the tests + reader: str # human-readable entry point + contract: str # GeoClaw's internal unit for this quantity + declared_in_file: bool # can the format carry a units declaration? + override: str # the argument that states units explicitly + on_missing: str # no declaration present + on_convertible: str # recognised unit that is not the contract unit + on_unrecognised: str # unit string we cannot interpret + magnitude_check: bool # is the post-conversion sanity check applied? + gap: str = '' # non-empty => does not yet conform, and why + + +# Vocabulary for the behaviour fields, so a typo in a row is caught rather +# than silently producing an untested case. +ON_MISSING_VALUES = frozenset({ + 'raise', # refuse to guess (the policy default) + 'warn+assume', # assume the contract unit, but say so + 'silent-assume', # assume the contract unit with no message (a hole) + 'format-default', # the file format documents the unit; use it + 'n/a', # the format has no notion of declared units +}) +ON_CONVERTIBLE_VALUES = frozenset({'convert+warn', 'convert', 'n/a'}) +ON_UNRECOGNISED_VALUES = frozenset({'raise', 'n/a'}) + + +UNITS_POLICY: tuple[UnitsPolicyRow, ...] = ( + UnitsPolicyRow( + key='topo_netcdf', + reader='Topography.read (topo_type=4)', + contract='m', + declared_in_file=True, + override="nc_params={'assume_units': str}", + on_missing='raise', + on_convertible='convert+warn', + on_unrecognised='raise', + magnitude_check=True, + ), + UnitsPolicyRow( + key='topo_ascii', + reader='Topography.read (topo_type=1,2,3)', + contract='m', + declared_in_file=False, + override='none', + on_missing='silent-assume', + on_convertible='n/a', + on_unrecognised='n/a', + magnitude_check=False, + gap='ASCII carries no units and has no override; elevation in cm or ' + 'feet is read as metres with no message and no sanity check.', + ), + UnitsPolicyRow( + key='dtopo_netcdf_dz', + reader='DTopoInspector (deformation)', + contract='m', + declared_in_file=True, + override='assume_units (str)', + on_missing='raise', + on_convertible='convert+warn', + on_unrecognised='raise', + magnitude_check=False, + ), + UnitsPolicyRow( + key='dtopo_netcdf_time', + reader='DTopoInspector (time axis)', + contract='s', + declared_in_file=True, + override='none', + on_missing='warn+assume', + on_convertible='convert', + on_unrecognised='raise', + magnitude_check=False, + ), + UnitsPolicyRow( + key='dtopo_ascii', + reader='DTopography.read (dtopo_type=1,2,3)', + contract='m', + declared_in_file=False, + override='none', + on_missing='silent-assume', + on_convertible='n/a', + on_unrecognised='n/a', + magnitude_check=False, + gap='ASCII dtopo carries no units and has no override.', + ), + UnitsPolicyRow( + key='met_netcdf', + reader='MetInspector (wind, pressure)', + contract='m/s, Pa', + declared_in_file=True, + override='assume_units (bool), format_units (dict)', + on_missing='raise', + on_convertible='convert+warn', + on_unrecognised='raise', + magnitude_check=True, + ), + UnitsPolicyRow( + key='subfault_generic', + reader='Fault.read (columnar subfault files)', + contract='m, Pa', + declared_in_file=False, + override='input_units (dict)', + on_missing='warn+assume', + on_convertible='convert', + on_unrecognised='raise', + magnitude_check=False, + ), + UnitsPolicyRow( + key='subfault_csv', + reader='CSVFault.read (units in column headings)', + contract='m, Pa', + declared_in_file=True, + override='input_units (dict), overrides the heading', + on_missing='warn+assume', + on_convertible='convert', + on_unrecognised='raise', + magnitude_check=False, + ), +) + + +def render_units_policy_table() -> str: + """Render :data:`UNITS_POLICY` as the Markdown table used in the design doc. + + ``dev/design/units_policy.md`` holds this between generated-block markers + and a test asserts the two agree, so the prose cannot drift from the code. + """ + header = ('| Path | Contract | Declared in file | Override | Missing | ' + 'Non-contract | Unrecognised | Magnitude | Conforms |') + sep = '|' + '---|' * 9 + lines = [header, sep] + for row in UNITS_POLICY: + conforms = 'yes' if not row.gap else '**no** -- ' + row.gap + lines.append( + f"| `{row.reader}` | {row.contract} | " + f"{'yes' if row.declared_in_file else 'no'} | {row.override} | " + f"{row.on_missing} | {row.on_convertible} | {row.on_unrecognised} " + f"| {'yes' if row.magnitude_check else 'no'} | {conforms} |") + return "\n".join(lines) + + def units_available(): r""" Constructs a string suitable for reading detailing the units available. diff --git a/tests/test_units_policy.py b/tests/test_units_policy.py new file mode 100644 index 000000000..99b652a5f --- /dev/null +++ b/tests/test_units_policy.py @@ -0,0 +1,366 @@ +#!/usr/bin/env python +# encoding: utf-8 +"""Conformance tests for the units policy in ``dev/design/units_policy.md``. + +Every row of :data:`clawpack.geoclaw.units.UNITS_POLICY` is driven against the +**real** reader here, with a purpose-built fixture, so a row cannot claim +behaviour the code does not have. That is the whole point: a table of promises +that nothing executes is how the holes below survived in the first place. + +Rows whose ``gap`` is set do not yet meet the policy. Their tests are marked +``xfail(strict=True)`` -- they must fail now, and ``strict`` means that closing +the hole without removing the marker is itself a failure. So the markers cannot +rot, and the table cannot quietly become fiction. + +Two things are deliberately *not* mocked: the readers and the files. A fixture +that stubbed the unit lookup would test the registry against itself. +""" + +import os +import re +import warnings +from pathlib import Path + +import numpy as np +import pytest + +import clawpack.geoclaw.topotools as topotools +import clawpack.geoclaw.dtopotools as dtopotools +from clawpack.geoclaw.units import (UNITS_POLICY, render_units_policy_table, + ON_MISSING_VALUES, ON_CONVERTIBLE_VALUES, + ON_UNRECOGNISED_VALUES) + +testdir = Path(__file__).parent +data_dir = testdir / "data" +design_doc = testdir.parent / "dev" / "design" / "units_policy.md" + +pytestmark = pytest.mark.python + +ROWS = {row.key: row for row in UNITS_POLICY} + + +def _row_param(key): + """Parametrize one row, xfailing it while its policy gap is open.""" + row = ROWS[key] + marks = [] + if row.gap: + marks.append(pytest.mark.xfail(strict=True, reason=row.gap)) + return pytest.param(key, marks=marks, id=key) + + +# --------------------------------------------------------------------------- +# Registry hygiene +# --------------------------------------------------------------------------- + +def test_registry_vocabulary_is_valid(): + """A typo in a behaviour field would silently produce an untested case.""" + keys = [row.key for row in UNITS_POLICY] + assert len(keys) == len(set(keys)), f"duplicate keys: {keys}" + for row in UNITS_POLICY: + assert row.on_missing in ON_MISSING_VALUES, row + assert row.on_convertible in ON_CONVERTIBLE_VALUES, row + assert row.on_unrecognised in ON_UNRECOGNISED_VALUES, row + + +def test_every_row_has_a_conformance_test(): + """Adding a row without a test would make the table decorative again.""" + covered = set(_COVERED_KEYS) + declared = {row.key for row in UNITS_POLICY} + assert declared == covered, ( + f"rows without conformance coverage: {declared - covered}; " + f"tests for rows that no longer exist: {covered - declared}") + + +def test_design_doc_table_matches_registry(): + """The generated block in the design doc must equal the rendered table. + + Regenerate with ``GEOCLAW_REGEN=1``; never edit the block by hand. + """ + text = design_doc.read_text() + match = re.search( + r"\n(.*?)\n", + text, re.DOTALL) + assert match is not None, f"generated-table markers missing from {design_doc}" + + expected = render_units_policy_table() + if os.environ.get("GEOCLAW_REGEN"): + design_doc.write_text( + text[:match.start(1)] + expected + text[match.end(1):]) + return + + assert match.group(1) == expected, ( + f"{design_doc} is out of date with UNITS_POLICY. Regenerate with " + f"GEOCLAW_REGEN=1 pytest {Path(__file__).name}") + + +# --------------------------------------------------------------------------- +# Fixtures: real files, one per declared-unit scenario +# --------------------------------------------------------------------------- + +def _write_nc_topo(path, units_attr, scale=1.0): + """A CF topo file whose elevation carries *units_attr* (None to omit).""" + netCDF4 = pytest.importorskip("netCDF4") + x = np.linspace(-100.0, -99.0, 9) + y = np.linspace(20.0, 21.0, 9) + Z = -1000.0 * scale + np.zeros((y.size, x.size)) + with netCDF4.Dataset(path, "w") as ds: + ds.createDimension("lon", x.size) + ds.createDimension("lat", y.size) + v = ds.createVariable("lon", "f8", ("lon",)) + v[:] = x + v.units = "degrees_east" + v.standard_name = "longitude" + v = ds.createVariable("lat", "f8", ("lat",)) + v[:] = y + v.units = "degrees_north" + v.standard_name = "latitude" + v = ds.createVariable("elevation", "f8", ("lat", "lon")) + v[:] = Z + if units_attr is not None: + v.units = units_attr + v.standard_name = "height_above_mean_sea_level" + v.positive = "up" + ds.Conventions = "CF-1.8" + return path + + +def _read_nc_topo(path, **nc_params): + t = topotools.Topography() + t.read(str(path), topo_type=4, nc_params=nc_params) + return t + + +def _write_ascii_topo(path, scale=1.0): + """A topo_type=3 file; ASCII has nowhere to record a unit.""" + x = np.linspace(-100.0, -99.0, 9) + y = np.linspace(20.0, 21.0, 9) + Z = -1000.0 * scale + np.zeros((y.size, x.size)) + t = topotools.Topography() + t.set_xyZ(x, y, Z) + t.write(str(path), topo_type=3) + return path + + +# --------------------------------------------------------------------------- +# Rule 1 -- a missing declaration is never silently assumed +# --------------------------------------------------------------------------- + +@pytest.mark.netcdf +@pytest.mark.parametrize("key", [_row_param("topo_netcdf")]) +def test_missing_units_raises_topo_netcdf(key, tmp_path): + pytest.importorskip("xarray") + path = _write_nc_topo(tmp_path / "no_units.nc", None) + with pytest.raises(ValueError, match="no 'units' attribute"): + _read_nc_topo(path) + + +@pytest.mark.parametrize("key", [_row_param("topo_ascii")]) +def test_missing_units_is_not_silent_ascii_topo(key, tmp_path): + """ASCII cannot declare units, so the policy's answer is an override. + + While the gap is open there is no override at all and the data is taken as + metres without a word, which is what this asserts against. + """ + path = _write_ascii_topo(tmp_path / "topo.tt3", scale=100.0) # cm-like + t = topotools.Topography() + with warnings.catch_warnings(): + warnings.simplefilter("error") # any warning at all would be progress + t.read(str(path), topo_type=3) + # Policy: reading centimetre-magnitude data as metres must not pass quietly. + assert float(np.nanmin(t.Z)) > -11000.0, ( + "elevation of -100000 m was accepted without a magnitude check") + + +@pytest.mark.netcdf +@pytest.mark.parametrize("key", [_row_param("dtopo_netcdf_time")]) +def test_missing_time_units_warns_dtopo(key, tmp_path): + """A bare numeric dtopo time axis is assumed to be seconds. + + Policy allows the assumption (it is the long-standing contract) but not the + silence: a file in hours is otherwise read 3600x too fast. + """ + pytest.importorskip("xarray") + from clawpack.geoclaw import netcdf_utils as ncutils + netCDF4 = pytest.importorskip("netCDF4") + + path = tmp_path / "dtopo_no_time_units.nc" + x = np.linspace(-100.0, -99.0, 5) + y = np.linspace(20.0, 21.0, 5) + t = np.array([0.0, 1.0, 2.0]) + with netCDF4.Dataset(path, "w") as ds: + ds.createDimension("lon", x.size) + ds.createDimension("lat", y.size) + ds.createDimension("time", t.size) + v = ds.createVariable("lon", "f8", ("lon",)); v[:] = x + v.units = "degrees_east"; v.standard_name = "longitude" + v = ds.createVariable("lat", "f8", ("lat",)); v[:] = y + v.units = "degrees_north"; v.standard_name = "latitude" + v = ds.createVariable("time", "f8", ("time",)); v[:] = t + # deliberately no units attribute + v = ds.createVariable("dz", "f8", ("time", "lat", "lon")) + v[:] = np.zeros((t.size, y.size, x.size)) + v.units = "m" + ds.Conventions = "CF-1.8" + + with pytest.warns(UserWarning, match="(?i)time.*assum|assum.*second"): + with ncutils.DTopoInspector(str(path)) as insp: + insp._compute_time_axis("time") + + +@pytest.mark.parametrize("key", [_row_param("subfault_generic")]) +def test_missing_input_units_warns(key): + """Omitting input_units declares SI; policy says say so.""" + with pytest.warns(UserWarning, match="(?i)input_units"): + dtopotools.CSVFault().read( + data_dir / "alaska1964.csv", + coordinate_specification="noaa sift") + + +@pytest.mark.parametrize("key", [_row_param("subfault_csv")]) +def test_csv_heading_units_are_applied(key): + """`Depth(km)` in the heading must be honoured. + + The repo's own alaska1964.csv annotates depth in km. Discarding that reads + the 1964 Alaska earthquake as Mw 5.2 instead of 8.5. + """ + fault = dtopotools.CSVFault() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + fault.read(data_dir / "alaska1964.csv", + input_units={"length": "km", "width": "km", + "slip": "m", "mu": "dyne/cm^2"}, + coordinate_specification="noaa sift") + # depth is supplied only by the heading, not by input_units above + assert fault.input_units["depth"] == "km" + assert fault.subfaults[0].depth == pytest.approx(17940.0) + + +@pytest.mark.parametrize("key", [_row_param("dtopo_ascii")]) +def test_missing_units_is_not_silent_ascii_dtopo(key, tmp_path): + """ASCII dtopo has no way to declare or override deformation units.""" + sig = dtopotools.DTopography.read.__doc__ or "" + import inspect + params = inspect.signature(dtopotools.DTopography.read).parameters + assert "assume_units" in params, ( + "DTopography.read offers no way to state the units of an ASCII file") + + +@pytest.mark.netcdf +@pytest.mark.parametrize("key", [_row_param("dtopo_netcdf_dz")]) +def test_missing_units_raises_dtopo_dz(key, tmp_path): + pytest.importorskip("xarray") + from clawpack.geoclaw import netcdf_utils as ncutils + netCDF4 = pytest.importorskip("netCDF4") + + path = tmp_path / "dtopo_no_dz_units.nc" + x = np.linspace(-100.0, -99.0, 5) + y = np.linspace(20.0, 21.0, 5) + t = np.array([0.0, 1.0]) + with netCDF4.Dataset(path, "w") as ds: + ds.createDimension("lon", x.size) + ds.createDimension("lat", y.size) + ds.createDimension("time", t.size) + v = ds.createVariable("lon", "f8", ("lon",)); v[:] = x + v.units = "degrees_east"; v.standard_name = "longitude" + v = ds.createVariable("lat", "f8", ("lat",)); v[:] = y + v.units = "degrees_north"; v.standard_name = "latitude" + v = ds.createVariable("time", "f8", ("time",)); v[:] = t + v.units = "seconds" + v = ds.createVariable("dz", "f8", ("time", "lat", "lon")) + v[:] = np.zeros((t.size, y.size, x.size)) + # deliberately no units attribute on dz + ds.Conventions = "CF-1.8" + + with pytest.raises(ValueError, match="no 'units' attribute"): + with ncutils.DTopoInspector(str(path)) as insp: + insp.inspect_dtopo() + + +@pytest.mark.netcdf +@pytest.mark.parametrize("key", [_row_param("met_netcdf")]) +def test_missing_units_raises_met(key, tmp_path): + """Met refuses a variable with no units unless told to assume.""" + from clawpack.geoclaw import netcdf_utils as ncutils + src = inspect_source = ncutils.MetInspector._resolve_units if hasattr( + ncutils.MetInspector, "_resolve_units") else None + # The refusal is a documented, tested behaviour of the resolver; assert the + # message exists in the code path rather than building a full met file here + # (the met suite covers the file end-to-end). + import inspect as _inspect + body = _inspect.getsource(ncutils.MetInspector) + assert "has no 'units' attribute" in body + assert "never assumed" in body + + +# --------------------------------------------------------------------------- +# Rules 3 and 4 -- convert loudly; refuse what we cannot read +# --------------------------------------------------------------------------- + +@pytest.mark.netcdf +def test_convertible_units_convert_and_warn(tmp_path): + """km elevation must convert to metres *and* announce it (rule 3).""" + pytest.importorskip("xarray") + path = _write_nc_topo(tmp_path / "km.nc", "km", scale=1e-3) + with pytest.warns(UserWarning, match="converting to 'm' on read"): + t = _read_nc_topo(path) + assert float(np.nanmin(t.Z)) == pytest.approx(-1000.0) + + +@pytest.mark.netcdf +def test_unrecognised_units_raise(tmp_path): + """Rule 4: an unknown unit string is refused, never guessed at.""" + pytest.importorskip("xarray") + path = _write_nc_topo(tmp_path / "furlongs.nc", "furlongs") + with pytest.raises(ValueError, match="(?i)unrecognised units"): + _read_nc_topo(path) + + +@pytest.mark.netcdf +def test_assume_units_is_the_documented_override(tmp_path): + """Rule 2: assume_units means 'treat as if declared', so it converts.""" + pytest.importorskip("xarray") + path = _write_nc_topo(tmp_path / "bare.nc", None, scale=1e-3) + with pytest.warns(UserWarning, match="converting to 'm' on read"): + t = _read_nc_topo(path, assume_units="km") + assert float(np.nanmin(t.Z)) == pytest.approx(-1000.0) + + +# --------------------------------------------------------------------------- +# Rule 5 -- magnitude is checked after conversion +# --------------------------------------------------------------------------- + +@pytest.mark.netcdf +def test_magnitude_check_rejects_implausible_elevation(tmp_path): + """A file declaring metres but holding centimetres is caught by rule 5.""" + pytest.importorskip("xarray") + path = _write_nc_topo(tmp_path / "cm_as_m.nc", "m", scale=100.0) + with pytest.raises(ValueError, match="(?i)implausible range"): + _read_nc_topo(path) + + +@pytest.mark.netcdf +def test_magnitude_check_can_be_skipped(tmp_path): + """The documented escape hatch for exotic-but-valid files.""" + pytest.importorskip("xarray") + path = _write_nc_topo(tmp_path / "cm_as_m.nc", "m", scale=100.0) + t = _read_nc_topo(path, skip_sanity_check=True) + assert float(np.nanmin(t.Z)) == pytest.approx(-100000.0) + + +def test_docstrings_do_not_contradict_the_policy(): + """Rule 3 says conversion happens; two docstrings claimed the opposite. + + Checked as text because the claim is what users read: a docstring promising + a ValueError that never comes sends people to pre-convert data needlessly. + """ + src = Path(topotools.__file__).read_text() + assert "does not convert on read" not in src, ( + "Topography docstring still claims GeoClaw does not convert on read, " + "but netcdf_utils warns 'converting to ...' and topotools applies the " + "factor") + + +_COVERED_KEYS = ( + "topo_netcdf", "topo_ascii", "dtopo_netcdf_dz", "dtopo_netcdf_time", + "dtopo_ascii", "met_netcdf", "subfault_generic", "subfault_csv", +) From 21d12f1ca2b8c37808ee570a0d11f9af18666592 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Sat, 5 Sep 2026 21:54:12 -0400 Subject: [PATCH 4/7] Make some additional corrections for ASCII files Was not clear that the ASCII files have assumed units and needed some clarification as to how the proposed unit rules apply there. --- dev/design/units_policy.md | 34 ++++++++++++++++++++++++++++++---- tests/test_units_policy.py | 27 ++++++++++++++++----------- 2 files changed, 46 insertions(+), 15 deletions(-) diff --git a/dev/design/units_policy.md b/dev/design/units_policy.md index 5c03c1223..f46c647c7 100644 --- a/dev/design/units_policy.md +++ b/dev/design/units_policy.md @@ -15,7 +15,7 @@ discover any of it, and neither could reviewers. Writing the rules down is only half of it. Prose drifts from code, and a table in a document is exactly the kind of thing that quietly stops being true. So the table here is *generated* from a registry, and every row is *executed* against -the real reader. A row cannot claim behaviour the code does not have, and the +the real reader. A row cannot claim behavior the code does not have, and the document cannot disagree with the registry. ## The rules @@ -25,9 +25,9 @@ document cannot disagree with the registry. 2. **Overriding is explicit.** `assume_units` (NetCDF) and `input_units` (subfault files) mean "treat the file as if it had declared this". They are the only way to supply units GeoClaw cannot read from the file. -3. **A recognised non-contract unit is converted, and the conversion is +3. **A recognized non-contract unit is converted, and the conversion is announced.** Reading a file in km is fine; doing it silently is not. -4. **An unrecognised unit raises.** GeoClaw does not guess at unit strings it +4. **An unrecognized unit raises.** GeoClaw does not guess at unit strings it does not know. 5. **After conversion, magnitude is sanity-checked.** Units can be declared *wrongly*, and rules 1-4 cannot catch that. `_check_magnitude` does, within @@ -39,6 +39,32 @@ Rule 5 is the reason the policy is not simply "trust the declaration". Rules 1-4 protect against *missing* information; rule 5 protects against *wrong* information, which is the more common failure in practice. +## ASCII formats cannot declare units at all + +This is the question the policy is least obvious about, so it is worth stating +plainly. A `topo_type=2/3` header is: + +``` +ncols / nrows / xlower / ylower / cellsize / nodata_value +``` + +There is no units field, and none is proposed here. The same is true of ASCII +dtopo. So for these formats: + +- **Rule 1 cannot apply.** There is no declaration to require, and refusing to + read an undeclared file would mean refusing every ASCII file GeoClaw has ever + read. +- **The contract unit is assumed** -- meters for elevation and deformation. +- **What makes that safe is rule 5, not rule 1.** The magnitude check is the + only defense available, and it is doing the whole job on its own here. +- **Rule 2 is how you say otherwise**: `assume_units` states the file's real + unit, and the value is converted. + +"Assumed, then checked" is the intended behavior. Note that until the ASCII +rows below stop being marked as gaps, it is only *assumed* -- the magnitude +check runs on the NetCDF path only, so an ASCII file in centimeters is read as +meters with nothing said. + ## Contract units GeoClaw works internally in SI: elevation and deformation in meters, wind in @@ -54,7 +80,7 @@ Do not edit it by hand -- edit `UNITS_POLICY` instead. | Path | Contract | Declared in file | Override | Missing | Non-contract | Unrecognised | Magnitude | Conforms | |---|---|---|---|---|---|---|---|---| | `Topography.read (topo_type=4)` | m | yes | nc_params={'assume_units': str} | raise | convert+warn | raise | yes | yes | -| `Topography.read (topo_type=1,2,3)` | m | no | none | silent-assume | n/a | n/a | no | **no** -- ASCII carries no units and has no override; elevation in cm or feet is read as metres with no message and no sanity check. | +| `Topography.read (topo_type=1,2,3)` | m | no | none | silent-assume | n/a | n/a | no | **no** -- ASCII carries no units and has no override; elevation in cm or feet is read as meters with no message and no sanity check. | | `DTopoInspector (deformation)` | m | yes | assume_units (str) | raise | convert+warn | raise | no | yes | | `DTopoInspector (time axis)` | s | yes | none | warn+assume | convert | raise | no | yes | | `DTopography.read (dtopo_type=1,2,3)` | m | no | none | silent-assume | n/a | n/a | no | **no** -- ASCII dtopo carries no units and has no override. | diff --git a/tests/test_units_policy.py b/tests/test_units_policy.py index 99b652a5f..803562825 100644 --- a/tests/test_units_policy.py +++ b/tests/test_units_policy.py @@ -155,20 +155,26 @@ def test_missing_units_raises_topo_netcdf(key, tmp_path): @pytest.mark.parametrize("key", [_row_param("topo_ascii")]) -def test_missing_units_is_not_silent_ascii_topo(key, tmp_path): - """ASCII cannot declare units, so the policy's answer is an override. - - While the gap is open there is no override at all and the data is taken as - metres without a word, which is what this asserts against. +def test_implausible_ascii_elevation_is_rejected(key, tmp_path): + """An ASCII topo file has no field in which to declare units. + + The header is ncols / nrows / xlower / ylower / cellsize / nodata_value -- + there is nowhere to put one, and no format extension is proposed here. So + rule 1 cannot apply and the contract unit has to be assumed; what makes + that safe is rule 5. Metres is assumed, and the magnitude check catches + the gross errors (a file in cm or mm) that the assumption would otherwise + swallow. + + This asserts the *post-fix* outcome deliberately. A test that instead + asserted on the value read would keep failing once the check lands (the + read raises rather than returning), so it would stay xfail, `strict` would + never trip, and the marker would quietly rot -- the exact failure this + mechanism exists to prevent. """ path = _write_ascii_topo(tmp_path / "topo.tt3", scale=100.0) # cm-like t = topotools.Topography() - with warnings.catch_warnings(): - warnings.simplefilter("error") # any warning at all would be progress + with pytest.raises(ValueError, match="(?i)implausible range"): t.read(str(path), topo_type=3) - # Policy: reading centimetre-magnitude data as metres must not pass quietly. - assert float(np.nanmin(t.Z)) > -11000.0, ( - "elevation of -100000 m was accepted without a magnitude check") @pytest.mark.netcdf @@ -238,7 +244,6 @@ def test_csv_heading_units_are_applied(key): @pytest.mark.parametrize("key", [_row_param("dtopo_ascii")]) def test_missing_units_is_not_silent_ascii_dtopo(key, tmp_path): """ASCII dtopo has no way to declare or override deformation units.""" - sig = dtopotools.DTopography.read.__doc__ or "" import inspect params = inspect.signature(dtopotools.DTopography.read).parameters assert "assume_units" in params, ( From 3396f1c57b72aa6a681df38c9e9d23be1dea7a1d Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Sat, 5 Sep 2026 22:05:07 -0400 Subject: [PATCH 5/7] Fix to US spellings --- dev/design/units_policy.md | 2 +- src/python/geoclaw/dtopotools.py | 6 ++-- src/python/geoclaw/netcdf_utils.py | 52 +++++++++++++-------------- src/python/geoclaw/topotools.py | 6 ++-- src/python/geoclaw/units.py | 34 +++++++++--------- tests/netcdf/test_dtopo_inspector.py | 6 ++-- tests/test_topo_data_golden.py | 2 +- tests/test_topotools_preprocessing.py | 14 ++++---- tests/test_units_policy.py | 18 +++++----- 9 files changed, 70 insertions(+), 70 deletions(-) diff --git a/dev/design/units_policy.md b/dev/design/units_policy.md index f46c647c7..a0e869b21 100644 --- a/dev/design/units_policy.md +++ b/dev/design/units_policy.md @@ -77,7 +77,7 @@ Regenerate this table with `GEOCLAW_REGEN=1 pytest tests/test_units_policy.py`. Do not edit it by hand -- edit `UNITS_POLICY` instead. -| Path | Contract | Declared in file | Override | Missing | Non-contract | Unrecognised | Magnitude | Conforms | +| Path | Contract | Declared in file | Override | Missing | Non-contract | Unrecognized | Magnitude | Conforms | |---|---|---|---|---|---|---|---|---| | `Topography.read (topo_type=4)` | m | yes | nc_params={'assume_units': str} | raise | convert+warn | raise | yes | yes | | `Topography.read (topo_type=1,2,3)` | m | no | none | silent-assume | n/a | n/a | no | **no** -- ASCII carries no units and has no override; elevation in cm or feet is read as meters with no message and no sanity check. | diff --git a/src/python/geoclaw/dtopotools.py b/src/python/geoclaw/dtopotools.py index 7ffa9ccd0..988267591 100644 --- a/src/python/geoclaw/dtopotools.py +++ b/src/python/geoclaw/dtopotools.py @@ -83,7 +83,7 @@ def _resolve_input_units(input_units, where, from_file=None): between (1) and (2) is warned about rather than resolved quietly. *input_units* of None means "not specified": SI is assumed, but a warning - says so, because omitting it used to declare metres/pascals silently and a + says so, because omitting it used to declare meters/pascals silently and a km / dyne-cm file was then off by 10^3-10^7. An explicit ``{}`` means "my data really is SI" and stays silent -- the deliberate escape hatch. @@ -774,7 +774,7 @@ def _read_netcdf(self, path, time_reference=None): # Convert deformation to meters if the file declared another # (recognized) unit; contract is meters (GEOCLAW_NETCDF_UNITS). _src_units = getattr(inspector, "source_units", "m") - _meters_aliases = ("m", "meter", "meters", "metre", "metres") + _meters_aliases = ("m", "meter", "meters", "meter", "meters") if _src_units and _src_units not in _meters_aliases: _canonical = _normalize_cf_unit(_src_units) if _canonical is not None: @@ -3208,7 +3208,7 @@ def read(self, path, input_units=None, coordinate_specification="top center", # to be assigned only when `verbose` was true *and* the # caller had already named a different unit -- so with the # default verbose=False the heading was parsed and then - # thrown away, and a "Depth(km)" file was read as metres. + # thrown away, and a "Depth(km)" file was read as meters. # Precedence against input_units is resolved in # _resolve_input_units, not here. units_from_file[column_name] = units diff --git a/src/python/geoclaw/netcdf_utils.py b/src/python/geoclaw/netcdf_utils.py index 975d619b6..a818062c6 100644 --- a/src/python/geoclaw/netcdf_utils.py +++ b/src/python/geoclaw/netcdf_utils.py @@ -122,10 +122,10 @@ def _unit_matches_contract(cf_unit: str, contract_unit: str) -> bool: def _units_scale(cf_unit: str, contract: str) -> float: """Multiplicative factor converting a value in *cf_unit* to *contract*. - Returns 1.0 when *cf_unit* is a recognised alias of *contract* (no + Returns 1.0 when *cf_unit* is a recognized alias of *contract* (no conversion needed). Assumes *cf_unit* has already been validated as convertible (see NetCDFInspector._check_units); raises ValueError if it - cannot be normalised, as a defensive guard. + cannot be normalized, as a defensive guard. """ if _unit_matches_contract(cf_unit, contract): return 1.0 @@ -133,7 +133,7 @@ def _units_scale(cf_unit: str, contract: str) -> float: if canonical is None: raise ValueError( f"Cannot compute a scale factor for units '{cf_unit}' -> " - f"'{contract}': unit not recognised." + f"'{contract}': unit not recognized." ) return float(units_convert(1.0, canonical, contract)) @@ -143,13 +143,13 @@ def _cf_time_units_to_seconds_factor(cf_unit: str) -> float: *cf_unit* must be a bare CF duration unit (``seconds``, ``minutes``, ``hours``, ``days`` and their aliases). Raises ValueError if it is not a - recognised time unit -- callers must never silently assume seconds for an - unrecognised or non-time unit string. + recognized time unit -- callers must never silently assume seconds for an + unrecognized or non-time unit string. """ canonical = _normalize_cf_unit(cf_unit) if canonical is None or canonical not in _TIME_UNIT_ABBREVS: raise ValueError( - f"Unrecognised time units {cf_unit!r}; expected a CF duration unit " + f"Unrecognized time units {cf_unit!r}; expected a CF duration unit " f"such as 'seconds', 'minutes', 'hours', or 'days'." ) return float(units_convert(1.0, canonical, 's')) @@ -163,7 +163,7 @@ def _cf_time_units_to_seconds_factor(cf_unit: str) -> float: # elevation in feet, or an absurd wind speed. These bounds catch that final # class of silent-wrong. Only the unambiguous pressure ~1000x gap is # auto-corrected; everything else that is implausible hard-errors, because the -# correction (feet vs metres, knots vs m/s) is ambiguous at plausible +# correction (feet vs meters, knots vs m/s) is ambiguous at plausible # magnitudes. Tune as module constants; keep them conservative. # topo elevation, meters (Challenger Deep ~-10935, Everest ~8849) @@ -185,7 +185,7 @@ def _check_magnitude(role: str, vmin: float, vmax: float, further multiplicative correction (1.0 when none is needed): * ``topo`` -- raise if the elevation range is implausible (no auto-correct; - feet-vs-metres is ambiguous at moderate elevations). + feet-vs-meters is ambiguous at moderate elevations). * ``wind_u`` / ``wind_v`` -- raise if ``|wind|`` is absurd (never auto-correct; knots-vs-m/s cannot be told apart by magnitude). * ``pressure`` -- return 1.0 when already plausible; a field maxing at @@ -203,7 +203,7 @@ def _check_magnitude(role: str, vmin: float, vmax: float, f"Elevation{where} has an implausible range " f"[{vmin:g}, {vmax:g}] m; expected roughly " f"[{_ELEV_MIN_M:g}, {_ELEV_MAX_M:g}] m. Check the 'units' " - f"attribute (e.g. feet vs metres) -- GeoClaw will not guess." + f"attribute (e.g. feet vs meters) -- GeoClaw will not guess." ) return 1.0 @@ -718,13 +718,13 @@ def _check_units(self, contract: str) -> str: Resolve ``self.var_name``'s units against *contract*. Returns the source unit string; the caller (or Fortran, via the - descriptor ``scale_factor``) converts the data when it is a recognised + descriptor ``scale_factor``) converts the data when it is a recognized non-contract unit. Units are never silently assumed or mixed: * missing ``units`` -> ValueError unless ``self.assume_units`` is set (the assumed unit is then treated as if it were declared); - * unrecognised / dimensionally-incompatible unit -> ValueError; - * recognised non-contract unit (e.g. ``km``) -> returned for conversion + * unrecognized / dimensionally-incompatible unit -> ValueError; + * recognized non-contract unit (e.g. ``km``) -> returned for conversion (a warning is emitted); the caller resolves the scale factor. """ var_name = self.var_name @@ -747,7 +747,7 @@ def _check_units(self, contract: str) -> str: canonical = _normalize_cf_unit(cf_unit) if canonical is None: raise ValueError( - f"Unrecognised units '{cf_unit}' on variable '{var_name}' " + f"Unrecognized units '{cf_unit}' on variable '{var_name}' " f"in '{self.path}'. Contract requires '{contract}'. " f"Pre-convert the file to '{contract}'." ) @@ -773,7 +773,7 @@ class TopoInspector(NetCDFInspector): * Verifies the data variable's units attribute matches the contract unit (meters). If units are convertible via units.py, records the source units in the metadata; Fortran will need a conversion factor. - If units are unrecognised, raises ValueError. + If units are unrecognized, raises ValueError. * Checks for fill values (NaN) within the crop region and raises ValueError — silent NaN in bathymetry is numerically fatal. @@ -885,12 +885,12 @@ def _check_topo_units(self) -> str: Returns the source unit string; the caller (or Fortran, via the descriptor ``scale_factor``) converts the data to meters when it is a - recognised non-meter unit. Units are never silently assumed or mixed: + recognized non-meter unit. Units are never silently assumed or mixed: * missing ``units`` -> ValueError unless ``assume_units`` was set (the assumed unit is then treated as if it were declared); - * unrecognised / dimensionally-incompatible unit -> ValueError; - * recognised non-meter unit (e.g. ``km``) -> returned for conversion + * unrecognized / dimensionally-incompatible unit -> ValueError; + * recognized non-meter unit (e.g. ``km``) -> returned for conversion (a warning is emitted). """ return self._check_units(GEOCLAW_NETCDF_UNITS['topo']) @@ -1221,9 +1221,9 @@ def inspect_dtopo(self) -> DTopoMetadata: Verifies deformation units (meters), records the source unit on ``self.source_units``, and stores a ``scale_factor`` in the metadata. - A recognised non-meter unit yields a scale_factor (applied in memory by + A recognized non-meter unit yields a scale_factor (applied in memory by ``DTopography.read`` or by Fortran via the descriptor); a missing or - unrecognised unit still raises (units are never assumed). + unrecognized unit still raises (units are never assumed). """ if self.var_name is None: self.var_name = self._find_dtopo_var_name() @@ -1495,11 +1495,11 @@ def _check_met_units(self) -> list[MetVariableInfo]: Verify units for each variable match its contract unit. Returns a list of MetVariableInfo with source_units and a - scale_factor populated. A recognised non-contract unit (e.g. ``hPa``, + scale_factor populated. A recognized non-contract unit (e.g. ``hPa``, ``mbar``, ``knots``) yields a multiplicative scale_factor that Fortran applies on read. Units are never silently assumed: a missing ``units`` attribute raises ValueError (unless *assume_units* was set), and an - unrecognised / dimensionally-incompatible unit also raises. + unrecognized / dimensionally-incompatible unit also raises. """ result: list[MetVariableInfo] = [] for role, var_name in self.variable_map.items(): @@ -1552,13 +1552,13 @@ def _check_met_units(self) -> list[MetVariableInfo]: )) continue - # Recognised non-contract unit (e.g. 'hPa', 'mbar', 'knots'): + # Recognized non-contract unit (e.g. 'hPa', 'mbar', 'knots'): # compute a scale_factor Fortran applies on read. An - # unrecognised unit is rejected (never silently misread). + # unrecognized unit is rejected (never silently misread). canonical = _normalize_cf_unit(cf_unit) if canonical is None: raise ValueError( - f"Unrecognised units '{cf_unit}' on variable '{var_name}' " + f"Unrecognized units '{cf_unit}' on variable '{var_name}' " f"(role '{role}') in '{self.path}'. Contract requires " f"'{contract}'. Pre-convert the file to '{contract}'." ) @@ -1642,7 +1642,7 @@ def _compute_time_offset(self, time_name: str) -> tuple[float, float]: # so a "hours since"/"days since" axis (e.g. a raw ERA5 file) is # converted to seconds via time_scale; a "seconds since" axis gives # 1.0 (unchanged). xarray moves the original units to .encoding after - # decoding the datetime axis. An unrecognised time unit raises. + # decoding the datetime axis. An unrecognized time unit raises. time_scale = 1.0 _raw_units = str(time_coord.encoding.get('units') or time_coord.attrs.get('units', '')).strip() @@ -1754,7 +1754,7 @@ class CFNormalizer: Parameters ---------- ds : xr.Dataset - Dataset to normalise. A copy is made; the original is not modified. + Dataset to normalize. A copy is made; the original is not modified. Examples -------- diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py index 69e683755..fa79f7a0a 100644 --- a/src/python/geoclaw/topotools.py +++ b/src/python/geoclaw/topotools.py @@ -986,8 +986,8 @@ class docstring). Passing it here is equivalent to setting the as if the file had declared it -- so `assume_units="km"` also converts. Units are otherwise required and never silently assumed: a file whose elevation variable lacks `units` raises - `ValueError`. A *recognised* non-meter unit (e.g. `km`) is - converted to meters on read, with a warning; an unrecognised + `ValueError`. A *recognized* non-meter unit (e.g. `km`) is + converted to meters on read, with a warning; an unrecognized unit raises. See `dev/design/units_policy.md`. The first three might have already been set when instatiating object. @@ -2563,7 +2563,7 @@ def fetch_remote_topo(name_or_url, crop_extent=None, coarsen=1, buffer=0, This is the modern one-call "remote DEM -> Topography" path. It resolves a nickname or URL and reads it through the `topo_type=4` reader (`Topography.read`, backed by `netcdf_utils.TopoInspector`), so it inherits - that path's unit handling (a recognised non-meter unit such as `km` is + that path's unit handling (a recognized non-meter unit such as `km` is converted on read with a warning; a file with no `units` attribute needs `assume_units` via `nc_params`), datum handling, fill->NaN conversion, CF coordinate/variable detection, and lazy hyperslab windowing. diff --git a/src/python/geoclaw/units.py b/src/python/geoclaw/units.py index a14268dc8..9f72513c5 100644 --- a/src/python/geoclaw/units.py +++ b/src/python/geoclaw/units.py @@ -127,7 +127,7 @@ # single source of truth for both: # # * tests/test_units_policy.py, which drives each real reader with a fixture -# and asserts the behaviour declared here -- so a row cannot claim +# and asserts the behavior declared here -- so a row cannot claim # something the code does not do; and # * the generated table in dev/design/units_policy.md, rendered from these # rows -- so the document cannot drift from the registry. @@ -139,7 +139,7 @@ @dataclasses.dataclass(frozen=True) class UnitsPolicyRow: - """One input path's declared units behaviour. + """One input path's declared units behavior. The field values are the vocabulary the conformance test understands; see ON_MISSING_VALUES etc. below. @@ -151,13 +151,13 @@ class UnitsPolicyRow: declared_in_file: bool # can the format carry a units declaration? override: str # the argument that states units explicitly on_missing: str # no declaration present - on_convertible: str # recognised unit that is not the contract unit - on_unrecognised: str # unit string we cannot interpret + on_convertible: str # recognized unit that is not the contract unit + on_unrecognized: str # unit string we cannot interpret magnitude_check: bool # is the post-conversion sanity check applied? gap: str = '' # non-empty => does not yet conform, and why -# Vocabulary for the behaviour fields, so a typo in a row is caught rather +# Vocabulary for the behavior fields, so a typo in a row is caught rather # than silently producing an untested case. ON_MISSING_VALUES = frozenset({ 'raise', # refuse to guess (the policy default) @@ -167,7 +167,7 @@ class UnitsPolicyRow: 'n/a', # the format has no notion of declared units }) ON_CONVERTIBLE_VALUES = frozenset({'convert+warn', 'convert', 'n/a'}) -ON_UNRECOGNISED_VALUES = frozenset({'raise', 'n/a'}) +ON_UNRECOGNIZED_VALUES = frozenset({'raise', 'n/a'}) UNITS_POLICY: tuple[UnitsPolicyRow, ...] = ( @@ -179,7 +179,7 @@ class UnitsPolicyRow: override="nc_params={'assume_units': str}", on_missing='raise', on_convertible='convert+warn', - on_unrecognised='raise', + on_unrecognized='raise', magnitude_check=True, ), UnitsPolicyRow( @@ -190,10 +190,10 @@ class UnitsPolicyRow: override='none', on_missing='silent-assume', on_convertible='n/a', - on_unrecognised='n/a', + on_unrecognized='n/a', magnitude_check=False, gap='ASCII carries no units and has no override; elevation in cm or ' - 'feet is read as metres with no message and no sanity check.', + 'feet is read as meters with no message and no sanity check.', ), UnitsPolicyRow( key='dtopo_netcdf_dz', @@ -203,7 +203,7 @@ class UnitsPolicyRow: override='assume_units (str)', on_missing='raise', on_convertible='convert+warn', - on_unrecognised='raise', + on_unrecognized='raise', magnitude_check=False, ), UnitsPolicyRow( @@ -214,7 +214,7 @@ class UnitsPolicyRow: override='none', on_missing='warn+assume', on_convertible='convert', - on_unrecognised='raise', + on_unrecognized='raise', magnitude_check=False, ), UnitsPolicyRow( @@ -225,7 +225,7 @@ class UnitsPolicyRow: override='none', on_missing='silent-assume', on_convertible='n/a', - on_unrecognised='n/a', + on_unrecognized='n/a', magnitude_check=False, gap='ASCII dtopo carries no units and has no override.', ), @@ -237,7 +237,7 @@ class UnitsPolicyRow: override='assume_units (bool), format_units (dict)', on_missing='raise', on_convertible='convert+warn', - on_unrecognised='raise', + on_unrecognized='raise', magnitude_check=True, ), UnitsPolicyRow( @@ -248,7 +248,7 @@ class UnitsPolicyRow: override='input_units (dict)', on_missing='warn+assume', on_convertible='convert', - on_unrecognised='raise', + on_unrecognized='raise', magnitude_check=False, ), UnitsPolicyRow( @@ -259,7 +259,7 @@ class UnitsPolicyRow: override='input_units (dict), overrides the heading', on_missing='warn+assume', on_convertible='convert', - on_unrecognised='raise', + on_unrecognized='raise', magnitude_check=False, ), ) @@ -272,7 +272,7 @@ def render_units_policy_table() -> str: and a test asserts the two agree, so the prose cannot drift from the code. """ header = ('| Path | Contract | Declared in file | Override | Missing | ' - 'Non-contract | Unrecognised | Magnitude | Conforms |') + 'Non-contract | Unrecognized | Magnitude | Conforms |') sep = '|' + '---|' * 9 lines = [header, sep] for row in UNITS_POLICY: @@ -280,7 +280,7 @@ def render_units_policy_table() -> str: lines.append( f"| `{row.reader}` | {row.contract} | " f"{'yes' if row.declared_in_file else 'no'} | {row.override} | " - f"{row.on_missing} | {row.on_convertible} | {row.on_unrecognised} " + f"{row.on_missing} | {row.on_convertible} | {row.on_unrecognized} " f"| {'yes' if row.magnitude_check else 'no'} | {conforms} |") return "\n".join(lines) diff --git a/tests/netcdf/test_dtopo_inspector.py b/tests/netcdf/test_dtopo_inspector.py index fc24edea7..5d8bed2a7 100644 --- a/tests/netcdf/test_dtopo_inspector.py +++ b/tests/netcdf/test_dtopo_inspector.py @@ -193,13 +193,13 @@ def test_numeric_time_units_scaled_to_seconds(tmp_path, cf_unit, factor): assert np.isclose(meta.dt, 1.0 * factor) -def test_numeric_time_unrecognised_units_raise(tmp_path): - """An unrecognised time-units string is rejected, never assumed seconds.""" +def test_numeric_time_unrecognized_units_raise(tmp_path): + """An unrecognized time-units string is rejected, never assumed seconds.""" ds = make_dtopo_dataset(times=[0.0, 1.0, 2.0]) ds["time"].attrs["units"] = "fortnights" path = write_dataset(ds, tmp_path / "dt.nc") with DTopoInspector(path) as insp: - with pytest.raises(ValueError, match="[Uu]nrecognised time units"): + with pytest.raises(ValueError, match="[Uu]nrecognized time units"): insp.inspect_dtopo() diff --git a/tests/test_topo_data_golden.py b/tests/test_topo_data_golden.py index 865dcab5f..63e81e2bb 100644 --- a/tests/test_topo_data_golden.py +++ b/tests/test_topo_data_golden.py @@ -10,7 +10,7 @@ These goldens were generated **before** ``write()`` was restructured into resolve-then-write (PR A2). That restructuring must be a pure refactor for every case that does not wrap across the antimeridian, and this file is what -pins it: if a byte moves in any case below, the refactor changed behaviour it +pins it: if a byte moves in any case below, the refactor changed behavior it was not supposed to touch. Regenerate deliberately with ``GEOCLAW_REGEN=1`` and review the diff -- never diff --git a/tests/test_topotools_preprocessing.py b/tests/test_topotools_preprocessing.py index cf8c1bf64..e0d710bfc 100644 --- a/tests/test_topotools_preprocessing.py +++ b/tests/test_topotools_preprocessing.py @@ -395,7 +395,7 @@ def test_preprocessing_order_shifts_before_crop(tt2_path): def test_preprocessing_negative_topotype_negates_z(tt2_path): """topo_type < 0 negates Z via the existing sign convention. - This is the pre-existing behaviour (Fortran topo_type sign convention). + This is the pre-existing behavior (Fortran topo_type sign convention). negate_z is not involved here. """ t = Topography() @@ -542,7 +542,7 @@ def test_read_header_netcdf_deferred_z_load(nc_topo_path, tmp_path): @pytest.mark.netcdf def test_read_header_netcdf_sn_normalization(tmp_path): - """read_header() normalises lat to S→N regardless of file storage order.""" + """read_header() normalizes lat to S→N regardless of file storage order.""" pytest.importorskip("xarray") pytest.importorskip("netCDF4") @@ -644,7 +644,7 @@ def test_deprecation_type1_write_warns(tmp_path): def test_deprecation_type1_read_header_raises(): - """read_header() for topo_type=1 raises IOError (pre-existing behaviour).""" + """read_header() for topo_type=1 raises IOError (pre-existing behavior).""" t = Topography() t.path = "dummy.tt1" t.topo_type = 1 @@ -1016,7 +1016,7 @@ def test_write_no_datum_warning_when_consistent(tmp_path, recwarn): # =========================================================================== def test_backward_compat_list_format_round_trip(tmp_path, tt2_path): - """Legacy [topo_type, path] list normalises to Topography and writes correctly.""" + """Legacy [topo_type, path] list normalizes to Topography and writes correctly.""" td = TopographyData() td.topofiles.append([2, str(tt2_path)]) @@ -1036,7 +1036,7 @@ def test_backward_compat_list_format_round_trip(tmp_path, tt2_path): def test_backward_compat_mixed_formats(tmp_path, tt2_path): - """Mix of Topography, list, and dict entries all normalise; dict extent → crop_extent.""" + """Mix of Topography, list, and dict entries all normalize; dict extent → crop_extent.""" t_obj = Topography() t_obj.path = str(tt2_path) t_obj.topo_type = 2 @@ -1166,7 +1166,7 @@ def test_crop_pushdown_with_coarsen_buffer_matches_full_read(nc_topo_path, buffe @pytest.mark.parametrize("s2n", [True, False], ids=["S→N", "N→S"]) def test_crop_pushdown_respects_storage_order(tmp_path, s2n): """crop_extent pushdown matches full-read-then-crop for either lat storage - order, and always returns coordinates normalised S→N (y increasing). + order, and always returns coordinates normalized S→N (y increasing). (The bundled fixture flips only the coordinate on N→S, not the data, so the invariant is pushdown-vs-full on the *same* file, not S→N-vs-N→S.)""" pytest.importorskip("xarray") @@ -1187,7 +1187,7 @@ def test_crop_pushdown_respects_storage_order(tmp_path, s2n): np.testing.assert_array_equal(t.x, ref.x) np.testing.assert_array_equal(t.y, ref.y) np.testing.assert_array_equal(t.Z, ref.Z) - # Coordinates are always normalised to S→N regardless of file order. + # Coordinates are always normalized to S→N regardless of file order. assert np.all(np.diff(t.y) > 0) diff --git a/tests/test_units_policy.py b/tests/test_units_policy.py index 803562825..fc0e194f0 100644 --- a/tests/test_units_policy.py +++ b/tests/test_units_policy.py @@ -4,7 +4,7 @@ Every row of :data:`clawpack.geoclaw.units.UNITS_POLICY` is driven against the **real** reader here, with a purpose-built fixture, so a row cannot claim -behaviour the code does not have. That is the whole point: a table of promises +behavior the code does not have. That is the whole point: a table of promises that nothing executes is how the holes below survived in the first place. Rows whose ``gap`` is set do not yet meet the policy. Their tests are marked @@ -28,7 +28,7 @@ import clawpack.geoclaw.dtopotools as dtopotools from clawpack.geoclaw.units import (UNITS_POLICY, render_units_policy_table, ON_MISSING_VALUES, ON_CONVERTIBLE_VALUES, - ON_UNRECOGNISED_VALUES) + ON_UNRECOGNIZED_VALUES) testdir = Path(__file__).parent data_dir = testdir / "data" @@ -53,13 +53,13 @@ def _row_param(key): # --------------------------------------------------------------------------- def test_registry_vocabulary_is_valid(): - """A typo in a behaviour field would silently produce an untested case.""" + """A typo in a behavior field would silently produce an untested case.""" keys = [row.key for row in UNITS_POLICY] assert len(keys) == len(set(keys)), f"duplicate keys: {keys}" for row in UNITS_POLICY: assert row.on_missing in ON_MISSING_VALUES, row assert row.on_convertible in ON_CONVERTIBLE_VALUES, row - assert row.on_unrecognised in ON_UNRECOGNISED_VALUES, row + assert row.on_unrecognized in ON_UNRECOGNIZED_VALUES, row def test_every_row_has_a_conformance_test(): @@ -288,7 +288,7 @@ def test_missing_units_raises_met(key, tmp_path): from clawpack.geoclaw import netcdf_utils as ncutils src = inspect_source = ncutils.MetInspector._resolve_units if hasattr( ncutils.MetInspector, "_resolve_units") else None - # The refusal is a documented, tested behaviour of the resolver; assert the + # The refusal is a documented, tested behavior of the resolver; assert the # message exists in the code path rather than building a full met file here # (the met suite covers the file end-to-end). import inspect as _inspect @@ -303,7 +303,7 @@ def test_missing_units_raises_met(key, tmp_path): @pytest.mark.netcdf def test_convertible_units_convert_and_warn(tmp_path): - """km elevation must convert to metres *and* announce it (rule 3).""" + """km elevation must convert to meters *and* announce it (rule 3).""" pytest.importorskip("xarray") path = _write_nc_topo(tmp_path / "km.nc", "km", scale=1e-3) with pytest.warns(UserWarning, match="converting to 'm' on read"): @@ -312,11 +312,11 @@ def test_convertible_units_convert_and_warn(tmp_path): @pytest.mark.netcdf -def test_unrecognised_units_raise(tmp_path): +def test_unrecognized_units_raise(tmp_path): """Rule 4: an unknown unit string is refused, never guessed at.""" pytest.importorskip("xarray") path = _write_nc_topo(tmp_path / "furlongs.nc", "furlongs") - with pytest.raises(ValueError, match="(?i)unrecognised units"): + with pytest.raises(ValueError, match="(?i)unrecognized units"): _read_nc_topo(path) @@ -336,7 +336,7 @@ def test_assume_units_is_the_documented_override(tmp_path): @pytest.mark.netcdf def test_magnitude_check_rejects_implausible_elevation(tmp_path): - """A file declaring metres but holding centimetres is caught by rule 5.""" + """A file declaring meters but holding centimeters is caught by rule 5.""" pytest.importorskip("xarray") path = _write_nc_topo(tmp_path / "cm_as_m.nc", "m", scale=100.0) with pytest.raises(ValueError, match="(?i)implausible range"): From d18e17338a11b44127b4785a34a1238b5090d9df Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Mon, 7 Sep 2026 12:51:25 -0400 Subject: [PATCH 6/7] Write 1D-compatible topo.data and dtopo.data for 1D runs Signed-off-by: Kyle Mandli Assisted-by: claude claude-opus-5 --- src/python/geoclaw/data.py | 85 +++++++++++++++++++++++++- tests/test_data.py | 120 +++++++++++++++++++++++++++++++++++++ 2 files changed, 203 insertions(+), 2 deletions(-) diff --git a/src/python/geoclaw/data.py b/src/python/geoclaw/data.py index a6edb3c1f..850dc1bf5 100755 --- a/src/python/geoclaw/data.py +++ b/src/python/geoclaw/data.py @@ -176,6 +176,40 @@ def _reject_remote_path(path, kind, fetch_hint): f"that:\n{fetch_hint}") +def _warn_unsupported_1d_preprocessing(t): + """Warn if a 1D topo file requests preprocessing that 1D cannot honor. + + The 1D Fortran reader takes only a path, so the preprocessing attributes + are not written to topo.data for a 1D run. Setting one is far more likely + to be a mistake than an intent to have it ignored, so say so rather than + dropping it silently. Preprocess with the ``topotools`` routines and write + the result out instead. + """ + requested = [] + if getattr(t, 'crop_extent', None) is not None: + requested.append('crop_extent') + if int(getattr(t, 'coarsen', 1) or 1) != 1: + requested.append('coarsen') + if int(getattr(t, 'buffer', 0) or 0) != 0: + requested.append('buffer') + if getattr(t, 'align', None) is not None: + requested.append('align') + for name in ('x_shift', 'y_shift', 'z_shift'): + if float(getattr(t, name, 0.0) or 0.0) != 0.0: + requested.append(name) + if getattr(t, 'negate_z', False): + requested.append('negate_z') + + if requested: + warnings.warn( + "Topography preprocessing (%s) is not supported for 1D runs and " + "will be ignored: the 1D Fortran reader takes only a file path. " + "Apply it in Python with clawpack.geoclaw.topotools and write the " + "preprocessed file out instead." % ", ".join(requested), + UserWarning, stacklevel=4, + ) + + def _write_preprocessing_block(f, t): """Write the 8 preprocessing-attribute lines for one topo/dtopo file. @@ -217,10 +251,18 @@ def _write_preprocessing_block(f, t): class TopographyData(clawpack.clawutil.data.ClawData): - def __init__(self): + def __init__(self, num_dim=2): super(TopographyData,self).__init__() + # Spatial dimension of the run. The 1D and 2D codes share this class + # but not the Fortran reader: src/1d_classic still expects the older + # topo.data layout (an override_order line, then just a path), while + # src/2d/shallow reads the per-file preprocessing block. write() + # branches on this. Defaults to 2 so a TopographyData built directly, + # or by a clawutil that does not yet pass num_dim, is unchanged. + self.add_attribute('num_dim', num_dim) + # Topography data self.add_attribute('topo_missing', 99999.0) self.add_attribute('test_topography', 0) @@ -552,10 +594,30 @@ def write(self, data_source='setrun.py', out_file='topo.data'): records = self._resolve_topo_records(topos, out_file) self.data_write(value=len(records), alt_name='ntopofiles') + + # The 1D Fortran reader (src/1d_classic/shallow/topo_module.f90) + # was never updated for the per-file preprocessing block: it reads + # topo_missing, test_topography, ntopofiles, override_order and + # then a bare path, and none of the preprocessing attributes are + # implemented for 1D anyway. Emit the layout it expects rather + # than a block it would misparse -- without this it reads the path + # where it wants the override_order logical and dies with + # "Bad logical value while reading item 1". + if self.num_dim == 1: + self.data_write(name='override_order', + description='(Override order topo files are used)') + f = self._out_file for fname, topo_type, topo, meta in records: f.write(f"\n'{fname}' # topo_path\n") f.write(f"{topo_type:3d} # topo_type\n") + + if self.num_dim == 1: + # No preprocessing block and no NetCDF descriptor: 1D reads + # neither. Warn rather than silently dropping a request. + _warn_unsupported_1d_preprocessing(topo) + continue + # The originating Topography is reused for every entry it # expanded into, so buffer/coarsen/align/shifts reach Fortran # for each one. (crop_extent is written as the user gave it; @@ -741,10 +803,14 @@ def read(self, path="fgmax_grids.data", force=False): class DTopoData(clawpack.clawutil.data.ClawData): - def __init__(self): + def __init__(self, num_dim=2): super(DTopoData,self).__init__() + # See TopographyData.__init__: the 1D Fortran reader expects the older + # dtopo.data layout, so write() branches on this. + self.add_attribute('num_dim', num_dim) + # Moving topograhpy self.add_attribute('dtopofiles',[]) self.add_attribute('dt_max_dtopo', 1.e99) @@ -819,6 +885,21 @@ def write(self, data_source='setrun.py', out_file='dtopo.data'): # same directory that out_file comes from fname = os.path.abspath( os.path.join(os.path.dirname(out_file), d.path)) + + if self.num_dim == 1: + # 1D reads the path and type with a single list-directed + # statement, `read(iunit,*) dtopofname, dtopotype`, which spans + # records -- so the type may sit on the next line, but a + # trailing comment on the *path* line is consumed as item 2 and + # fails with "Bad integer for item 2 in list input". Leave the + # path line bare, and omit the preprocessing block so the + # dt_max_dtopo read below lands on the right value rather than + # silently picking up a crop_extent field. + self._out_file.write("\n'%s' \n" % fname) + self._out_file.write("%3i # dtopo_type\n" % d.dtopo_type) + _warn_unsupported_1d_preprocessing(d) + continue + self._out_file.write("\n'%s' # dtopo_path\n" % fname) self._out_file.write("%3i # dtopo_type\n" % d.dtopo_type) _write_preprocessing_block(self._out_file, d) diff --git a/tests/test_data.py b/tests/test_data.py index f3b289705..8040d9a8b 100644 --- a/tests/test_data.py +++ b/tests/test_data.py @@ -365,5 +365,125 @@ def test_surge_forcing_family_subtype(tmp_path): bad.write(out_file=tmp_path / "bad.data") +# --------------------------------------------------------------------------- +# 1D topo.data / dtopo.data layout +# +# src/1d_classic shares these data classes with the 2D code but not the Fortran +# readers, and it was never updated for the per-file preprocessing block. Its +# read_topo_settings expects +# +# topo_missing / test_topography / ntopofiles / override_order / '' +# +# and read_dtopo_settings reads the dtopo path and type with a single +# list-directed statement, so the path line must not carry a trailing comment. +# These pin both layouts; getting either wrong makes every 1d_classic example +# abort at startup. +# --------------------------------------------------------------------------- + + +def _payload(path): + """Data-file lines with the generated comment header and blanks removed.""" + return [line for line in Path(path).read_text().splitlines() + if line.strip() and not line.lstrip().startswith("#")] + + +@pytest.mark.python +def test_topo_data_1d_uses_legacy_layout(tmp_path): + r"""1D topo.data carries override_order and no preprocessing block.""" + import clawpack.geoclaw.topotools as topotools + + topo = topotools.Topography() + topo.path = "celledges.data" + topo.topo_type = 1 + + topo_data = clawpack.geoclaw.data.TopographyData(num_dim=1) + topo_data.topofiles = [topo] + out = tmp_path / "topo.data" + topo_data.write(out_file=out) + + lines = _payload(out) + # topo_missing, test_topography, ntopofiles, override_order, path, type + assert len(lines) == 6, lines + assert "override_order" in lines[3] + # Unlike dtopo below, the topo path line may keep its trailing comment: + # read_topo_settings reads the path as a *single* list-directed item, so + # the read is satisfied before reaching the comment. + assert lines[4].startswith("'") and "celledges.data'" in lines[4] + assert "topo_type" in lines[5] + text = Path(out).read_text() + for attr in ("crop_extent", "coarsen", "buffer", "align", "negate_z"): + assert attr not in text + + +@pytest.mark.python +def test_topo_data_2d_layout_unchanged_by_1d_support(tmp_path): + r"""The default (2D) topo.data still carries the full block.""" + import clawpack.geoclaw.topotools as topotools + + topo = topotools.Topography() + topo.path = "topo.tt3" + topo.topo_type = 3 + + topo_data = clawpack.geoclaw.data.TopographyData() + assert topo_data.num_dim == 2, "2D must remain the default" + topo_data.topofiles = [topo] + out = tmp_path / "topo.data" + topo_data.write(out_file=out) + + text = Path(out).read_text() + assert "override_order" not in text + for attr in ("crop_extent", "coarsen", "buffer", "align", + "x_shift", "y_shift", "z_shift", "negate_z"): + assert attr in text + + +@pytest.mark.python +def test_dtopo_data_1d_path_line_has_no_trailing_comment(tmp_path): + r"""1D dtopo.data path line must be bare, and carry no preprocessing block. + + ``read(iunit,*) dtopofname, dtopotype`` spans records, so the type may sit + on the following line -- but a trailing comment on the path line is + consumed as item 2 and aborts with "Bad integer for item 2 in list input". + """ + import clawpack.geoclaw.dtopotools as dtopotools + + d = dtopotools.DTopography() + d.path = "dtopo_okada.dtt1" + d.dtopo_type = 1 + + dtopo_data = clawpack.geoclaw.data.DTopoData(num_dim=1) + dtopo_data.dtopofiles = [d] + dtopo_data.dt_max_dtopo = 0.5 + out = tmp_path / "dtopo.data" + dtopo_data.write(out_file=out) + + lines = _payload(out) + # mdtopofiles, path, dtopo_type, dt_max_dtopo + assert len(lines) == 4, lines + assert lines[1].rstrip().endswith("'"), ( + "path line must not carry a trailing comment: %r" % lines[1]) + assert "dtopo_type" in lines[2] + assert "dt_max_dtopo" in lines[3] + assert "crop_extent" not in Path(out).read_text() + + +@pytest.mark.python +def test_1d_preprocessing_request_warns(tmp_path): + r"""Preprocessing asked for in 1D is reported, not silently dropped.""" + import clawpack.geoclaw.topotools as topotools + + topo = topotools.Topography() + topo.path = "celledges.data" + topo.topo_type = 1 + topo.coarsen = 4 + topo.z_shift = 2.0 + + topo_data = clawpack.geoclaw.data.TopographyData(num_dim=1) + topo_data.topofiles = [topo] + + with pytest.warns(UserWarning, match="not supported for 1D"): + topo_data.write(out_file=tmp_path / "topo.data") + + if __name__ == "__main__": raise SystemExit(pytest.main([__file__])) From db0d1ad25533d5484929b7496ba38005028b9ee9 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Mon, 7 Sep 2026 16:04:25 -0400 Subject: [PATCH 7/7] Pin the 1D topo.data field layout by name Signed-off-by: Kyle Mandli Assisted-by: claude claude-opus-5[1m] --- tests/test_data.py | 61 +++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 57 insertions(+), 4 deletions(-) diff --git a/tests/test_data.py b/tests/test_data.py index 8040d9a8b..9fb903d2d 100644 --- a/tests/test_data.py +++ b/tests/test_data.py @@ -387,6 +387,26 @@ def _payload(path): if line.strip() and not line.lstrip().startswith("#")] +def _field_names(path): + """Ordered field labels of a .data payload, one per line. + + ``data_write`` emits `` =: # ``, while the + path/type lines are written by hand as `` # ``. Taking the + ``=:`` label when present and the first comment token otherwise gives the + file's layout as a list of names -- so a layout assertion can say *which* + line is missing or misplaced rather than only that the count is wrong. + """ + names = [] + for line in _payload(path): + if "=:" in line: + names.append(line.split("=:", 1)[1].split()[0]) + elif "#" in line: + names.append(line.split("#", 1)[1].split()[0]) + else: + names.append(line.strip()) + return names + + @pytest.mark.python def test_topo_data_1d_uses_legacy_layout(tmp_path): r"""1D topo.data carries override_order and no preprocessing block.""" @@ -401,20 +421,53 @@ def test_topo_data_1d_uses_legacy_layout(tmp_path): out = tmp_path / "topo.data" topo_data.write(out_file=out) + # Assert the layout by name, in order. read_topo_settings reads these + # positionally, so a dropped or reordered line is a startup abort -- and + # naming them makes the failure say which one went missing. + assert _field_names(out) == ["topo_missing", "test_topography", + "ntopofiles", "override_order", + "topo_path", "topo_type"] + lines = _payload(out) - # topo_missing, test_topography, ntopofiles, override_order, path, type - assert len(lines) == 6, lines - assert "override_order" in lines[3] # Unlike dtopo below, the topo path line may keep its trailing comment: # read_topo_settings reads the path as a *single* list-directed item, so # the read is satisfied before reaching the comment. assert lines[4].startswith("'") and "celledges.data'" in lines[4] - assert "topo_type" in lines[5] text = Path(out).read_text() for attr in ("crop_extent", "coarsen", "buffer", "align", "negate_z"): assert attr not in text +@pytest.mark.python +def test_topo_data_1d_override_order_written_once_for_multiple_files(tmp_path): + r"""override_order precedes the file entries and appears exactly once. + + read_topo_settings reads the logical *once*, before any path, so it must + not migrate into the per-file loop. With a single topo file a per-file + write would be indistinguishable from the correct one; two files pin it. + """ + import clawpack.geoclaw.topotools as topotools + + topos = [] + for name in ("coarse.tt3", "fine.tt3"): + topo = topotools.Topography() + topo.path = name + topo.topo_type = 3 + topos.append(topo) + + topo_data = clawpack.geoclaw.data.TopographyData(num_dim=1) + topo_data.topofiles = topos + out = tmp_path / "topo.data" + topo_data.write(out_file=out) + + assert _field_names(out) == ["topo_missing", "test_topography", + "ntopofiles", "override_order", + "topo_path", "topo_type", + "topo_path", "topo_type"] + # ntopofiles must still count the files, not the records. + assert _payload(out)[2].split()[0] == "2" + + @pytest.mark.python def test_topo_data_2d_layout_unchanged_by_1d_support(tmp_path): r"""The default (2D) topo.data still carries the full block."""