Skip to content

Topo crop_extent / antimeridian: stop failing silently, and make the working paths reachable - #739

Open
mandli wants to merge 2 commits into
clawpack:masterfrom
mandli:topo-crop-silent-failures
Open

Topo crop_extent / antimeridian: stop failing silently, and make the working paths reachable#739
mandli wants to merge 2 commits into
clawpack:masterfrom
mandli:topo-crop-silent-failures

Conversation

@mandli

@mandli mandli commented Sep 5, 2026

Copy link
Copy Markdown
Member
  1. Make the silent failures loud — 13 cases that produced a wrong answer with no message.
  2. Make the remote and cross-seam paths reachable from setrun — two failures reported from the field, where the capability existed, was tested, and could not be reached from ordinary user code.

The second commit updates one error message introduced by the first, so they read best in order.


Claude Description of PR

Part 1 — every one of these produced a wrong answer with no message

None of them raised, warned, or logged. Several returned an empty or full-file grid where a crop was asked for; one made GeoClaw abort on correctly specified input.

# Before Now
1 crop_extent=[170, -170] produced an empty Topography (Z.shape=(11, 0)); .extent then raised an opaque numpy "zero-size array to reduction" far from the cause ValueError naming both spellings
2 A crop overlapping the extent but narrower than one cell returned the full file ValueError reporting dx/dy
3 crop_extent=[-211, -99] on a [-180, 180] file silently clipped 112° → 81° Still clipped (legitimate), but warns with the amount
4 Non-overlapping crop: Python silently kept the full file, Fortran does stop 1 Fallback kept (a test pins it) but warns, on both the ASCII and NetCDF paths
5 crop() reported via print() — uncatchable, unfilterable warnings.warn
6 read(unstructured=True, crop_extent=...) died with TypeError: list indices must be integers or slices, not tuple NotImplementedError, matching what crop() already did
7 topo_type=5 (GeoTIFF) written to topo.data; Fortran has no case(5) and aborts Warns at write time
8 topo_type=None hit the :3d format and raised a TypeError naming neither file nor attribute — reachable via the recommended Topography-object path Falls back to determine_topo_type, else raises naming the attribute
9 A cross-seam crop_extent on a type-4 file wrote crop_bounds = 170.0 -170.0; Fortran resolves that to mx=0, my=0 — an empty topography, no error Raises (and see Part 2: the continuous spelling now just works)
10 read_header() then a cropped type-4 read() left .extent reporting the full file — and .extent is what _compute_priority_order and the plotting routines use _extent/_delta invalidated alongside _X/_Y
11 (Fortran) buffer silently ignored for every file carrying descriptor crop_bounds — i.e. every file topo_entries() writes nbuf4 set on that branch too

Item 10 was suspected by reading, not running; it is confirmed here with a NetCDF fixture (.extent returned [-100, -80, 20, 30] for data covering [-95, -90, 22, 25]).

The Fortran hunk (item 11), measured

One line in topo_module.f90 — file I/O and index arithmetic, no solver numerics. It was asked for as a measurement rather than an assurance, and it is not small. On a case whose domain sits in the buffer ring:

before:  **** topo arrays do not cover domain
         **** area of overlap =  14.6025
         **** area of domain  =  17.4900     -> GeoClaw stops, 21 lines of output
after:   runs to completion, 80 lines

83.5% coverage; the missing 16.5% is exactly the buffer ring that was discarded. The failure mode was a hard abort on correctly specified input, not a small numerical shift.


Part 2 — the capability existed; it just wasn't reachable

Reported from the field, writing the most natural possible setrun code:

topo.path = 'https://.../ETOPO_2022_v1_30s_N90W180_bed.nc'   # THREDDS
topo.topo_type = 4
topo.crop_extent = [-160, -120, -60, 0]
topo.buffer = 1
topo.coarsen = 20
rundata.topo_data.topofiles.append(topo)
  1. make dataFileNotFoundError: /run/dir/https:/www.ngdc.noaa.gov/... — a path the user never typed.
  2. With a local global file and crop_extent = [-190, -120, -60, 0] to span the date line → ValueError: crop_bounds lon [-190.0, -120.0] exceed file extent [-180.0, 180.0].

Part 1 fixes neither. The second crop is ascending, so the new descending-crop check never fires.

  • NetCDFInspector.__init__ has guarded against the https://https:/ collapse since Add fetch_remote_topo; deprecate read_netcdf; refactor etopotools #726 — but TopographyData.write() ran the path through os.path.abspath first, corrupting it before the guard saw it. Same bug in DTopoData.write().
  • _compute_lon_entries already covers a cross-seam crop. For (-180, 180, -190, -120) it returns (-180, -120, offset 0) and (170, 180, offset -360). topo_entries() wraps it and _normalize_topofiles already accepts its output as a first-class entry form. Nothing called it from write().

Remote URLs rejected with the working recipe

The scheme-anchored regex from NetCDFInspector.__init__ is now netcdf_utils.is_remote_url(), shared rather than repeated. Both .data writers check before abspath and raise with fetch_remote_topo's own documented pattern (fetch the hyperslab → write a local file → append that).

Pass-through was considered and rejected: it needs netcdf-fortran built with DAP, which is not guaranteed, and it makes every run depend on the server.

A cross-seam crop is split across the seam automatically

TopographyData.write() becomes resolve-then-writentopofiles is written before the blocks, and the entry count is no longer known until every file is resolved. A type-4 crop running off the file's longitude extent now produces one entry per side:

2                    =: ntopofiles
  ... lon_wrap_offset = 0.0      crop_bounds = -180.0 -120.0 -60.0 0.0
  ... lon_wrap_offset = -360.0   crop_bounds =  170.0  180.0 -60.0 0.0

Three things this had to get right:

  • buffer and coarsen reach every entry. topo_entries() builds Topography objects carrying only _netcdf_meta, so routing through it naively writes buffer = 0 — which would have silently undone item 11 above, the change that made buffer work for descriptor crops at all.
  • No full-file scan. topo_entries() inspects with crop_bounds unset (it must — a wrapping crop lies outside the file extent by construction), so its fill and magnitude checks read the entire variable and reject NaN anywhere in it. On a global DEM over OPeNDAP that turns make data into a full download. inspect_topo/topo_entries gain fill_scan (default True, so no existing caller changes) and write() passes False, preserving exactly its previous cost and behaviour.
  • Latitude is still validated. Longitude wraps; latitude has no seam, so a crop running off the file in latitude is now an explicit error.

lon_wrap_offset comes from the computed entry rather than the hard-coded 0.0 — wrong even for a single entry, since [185, 195] on a [-180, 180] file needs offset +360.

Non-wrapping cases take the identical code path as before; the resolver only diverges once it has established that the crop runs off the file.


Verification

A byte-exact topo.data golden suite, written before the refactor. write() runs on every make data, and the existing tests check substrings and line counts only — the write paths could diverge in field width, whitespace or float formatting and still pass. Five cases are pinned byte-for-byte.

I generated them from the pre-refactor code, regenerated after, and diffed the two sets: identical. So the restructure is a demonstrated pure refactor, not an asserted one. (I first regenerated them post-refactor, which would have defeated the purpose; the diff is what caught it.)

25 new Python tests, plus test_crop_no_overlap_keeps_full_grid updated to assert the warning it now emits. Every one was checked for the ability to fail by reverting the source: 14 of 16 in Part 1 and all 9 in Part 2. The two that pass either way are the two that should — a guard against the new clipping warning becoming noise, and one pinning the existing, correct buffer-before-coarsen order.

Two Part 2 tests reproduce the reported errors verbatim. One is worth calling out: test_wrapping_write_does_not_scan_the_whole_file spies on _check_fill_in_crop and asserts it is never called — that regression is invisible on a small local fixture and ruinous on a remote DEM.

A new Fortran regression suite, tests/regression/topo_crop/ — the first tests that compile GeoClaw and feed it a NetCDF topo file through a descriptor from topo_entries(). Five cases, including two that feed GeoClaw two entries for the same file with different lon_wrap_offset: the Python side was tested, the Fortran side was assumed. One asserts a date-line-straddling domain is fully covered; the other asserts both grids are reported in domain coordinates in fort.geo, so the wrapped entry appears at -190 rather than at its file position of +170. Verified by sabotage — forcing the offset to 0.0 fails both (assert -180.0 == -190.0).

Three fixture defects were found and fixed by that exercise, each of which would have made a test that pins nothing:

  • The unstructured test used a headed .tt3, so without the fix it died in numpy.loadtxt on the header — an unrelated parse error, not the TypeError it claims to pin. Now a real .xyz, and the pre-fix failure is verified.
  • A NetCDF fixture used arange(..., 0.05), whose spacing comes out as 4.9999999999997e-2; the inclusive crop_bounds comparison then landed one index differently on each edge. DELTA is now 0.125, exact in binary.
  • A test domain whose bounds coincided exactly with the topo bounds passed GeoClaw's float coverage test and then crashed intermittently.

Results: 546 passed / 4 skipped / 1 xfailed; 12 regression tests (7 met_forcing + 5 topo_crop). All met_forcing goldens byte-identical.

buffer before coarsen — correct as implemented, now tested and documented

Confirmed by reading both sides: _crop_indices and Fortran apply_align_buffer_coarsen agree, and buffer counts coarsened output pointsbuffer=2, coarsen=4 adds 2 points per edge, not 8. Nothing pinned that (the existing combined test is a relative netCDF-vs-ASCII equality, which passes for either meaning). Now asserted directly and documented in read(), crop() and the class docstring rather than only in a private helper.

Documentation

The Topography docstring gains two sections aimed at the question of what a cropped Topography represents across the seam. The short version: a Topography never wraps — wrapping is applied by the Fortran reader from a lon_wrap_offset that only exists in a NetCDF descriptor. The three cases are spelled out, including that reading and writing deliberately differ (reading clips and warns; writing splits across the seam), plus the crop → align → buffer → coarsen order.

Behaviour changes worth a release note

  • A descending crop_extent now raises (previously an empty grid).
  • A sub-cell crop_extent now raises (previously the full file).
  • unstructured=True with preprocessing raises NotImplementedError (previously a TypeError, so no working code is affected).
  • A URL in topofiles/dtopofiles now raises with instructions (previously a FileNotFoundError naming a mangled path).
  • A cross-seam crop_extent in continuous coordinates now works, writing two entries where it previously raised.
  • The Fortran nbuf4 fix changes what data a descriptor-cropped type-4 read returns, by adding the buffer that was requested.

Notes for reviewers

  • The goldens are .txt, not .data: .gitignore excludes *.data, so goldens named for what they are would have been silently left out of the commit and the test would have failed for everyone else. Same convention as the met_forcing goldens.
  • While building the Fortran harness I lost time to an intermittent segfault that looked like a topo bug. It was not: num_waves defaults to 1 while the GeoClaw Riemann solver always returns 3 and uses f-waves, so omitting num_waves = 3 / limiter = ['mc']*3 / use_fwaves = True lets the solver write past the end of the wave arrays. With them set: 0 crashes in 20 runs, where before it was 5–12 in 20. Recorded as a comment in the new setrun.py rather than left as an unexplained incantation.

Follow-ups, recorded not fixed

  • Scope the fill/magnitude scan to each entry's own crop instead of skipping it. fill_scan=False preserves today's behaviour but is a safety downgrade relative to scanning the crop properly. Belongs with topo-input-unify.
  • FrictionData.write is broken: it writes "'%s' %s\n " % fname — two format specifiers, one argument, so it raises TypeError whenever friction_files is non-empty, and that path has zero test coverage. Left alone because the feature is a known-dormant stub (friction_module.f90 prints "File based friction specification unimplemented") with its own tracked design doc.

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 <kyle.mandli@gmail.com>
Assisted-by: claude claude-opus-5[1m]
@mandli
mandli force-pushed the topo-crop-silent-failures branch from 63e8cf6 to c893f23 Compare September 5, 2026 01:04
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 <kyle.mandli@gmail.com>
Assisted-by: claude claude-opus-5[1m]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant