From 80de4c993a8734ce8b8965ffe30e4bd7afb1dbf0 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:19:45 -0400 Subject: [PATCH 1/2] Close the surface_distance backend test gaps and fix dask direction (#3725, #3713) The test file covered three public functions across four backends with 37 tests, but never ran anything on dask+cupy and only ever called surface_direction() on eager numpy. That second hole was hiding a real bug: _run_tile() records global source indices while _finalize_direction() built its pixel grid from block-local ones, so every tile except the top-left measured its bearings from the wrong origin. On a flat 6x6 raster chunked (3, 3) the iterative path returned the top-left tile's bearing field tiled across the whole output. _finalize_direction() and _extract_output() now take row/col offsets, and _assemble_sd() passes the tile's own offsets so both sides of the subtraction live in the same coordinate space. The eager numpy and bounded map_overlap paths keep the default of 0 because they seed local indices. New tests cover the holes the audit found: dask+cupy on all three public functions, surface_direction parity across every backend on both the bounded and iterative branches, Inf and all-NaN elevation, 1x1 and Nx1 rasters, connectivity=4 off the eager path, target_values on cupy, Dataset input through @supports_dataset, the geodesic NotImplementedError per backend, the dims rejection, and attrs/coords/dim preservation. Tests 37 -> 106, all passing with CUDA present. Branch coverage of xrspatial.surface_distance 80% -> 87% under NUMBA_DISABLE_JIT=1. --- xrspatial/surface_distance.py | 22 +- xrspatial/tests/test_surface_distance.py | 414 +++++++++++++++++++++++ 2 files changed, 431 insertions(+), 5 deletions(-) diff --git a/xrspatial/surface_distance.py b/xrspatial/surface_distance.py index c91aeb13e..b2b3c8fcc 100644 --- a/xrspatial/surface_distance.py +++ b/xrspatial/surface_distance.py @@ -357,13 +357,20 @@ def _finalize_alloc(alloc, dist, max_distance): def _finalize_direction(src_row, src_col, dist, cellsize_x, cellsize_y, - max_distance): + max_distance, row_offset=0, col_offset=0): """Compute compass bearing from each pixel to its allocated source. Uses pixel index differences scaled by cell size. + + ``src_row``/``src_col`` may hold indices in a wider coordinate space + than the block being finalized -- the dask iterative path records + global indices while each block sees only its own shape. The + offsets shift the pixel grid into that same space so the difference + is taken between comparable indices. """ H, W = dist.shape - row_idx, col_idx = np.meshgrid(np.arange(H), np.arange(W), indexing='ij') + row_idx, col_idx = np.meshgrid(np.arange(H) + row_offset, + np.arange(W) + col_offset, indexing='ij') # Coordinate differences (source - pixel) dx = (src_col.astype(np.float64) - col_idx) * cellsize_x @@ -381,7 +388,8 @@ def _finalize_direction(src_row, src_col, dist, cellsize_x, cellsize_y, def _extract_output(dist, alloc, src_row, src_col, - cellsize_x, cellsize_y, max_distance, mode): + cellsize_x, cellsize_y, max_distance, mode, + row_offset=0, col_offset=0): """Select and finalize the requested output from raw Dijkstra arrays.""" if mode == DISTANCE: return _finalize_dist(dist, max_distance) @@ -389,7 +397,8 @@ def _extract_output(dist, alloc, src_row, src_col, return _finalize_alloc(alloc, dist, max_distance) else: return _finalize_direction(src_row, src_col, dist, - cellsize_x, cellsize_y, max_distance) + cellsize_x, cellsize_y, max_distance, + row_offset, col_offset) # --------------------------------------------------------------------------- @@ -1174,8 +1183,11 @@ def _tile_fn(source_block, elev_block, block_info=None): max_distance, dy, dx, dd, row_offset, col_offset, ) + # _run_tile records global src_row/src_col, so the pixel grid has + # to be built in global coordinates too. return _extract_output(dist, alloc_arr, srow, scol, - cellsize_x, cellsize_y, max_distance, mode) + cellsize_x, cellsize_y, max_distance, mode, + row_offset, col_offset) return da.map_blocks( _tile_fn, diff --git a/xrspatial/tests/test_surface_distance.py b/xrspatial/tests/test_surface_distance.py index 74ae62fbd..c3b2db616 100644 --- a/xrspatial/tests/test_surface_distance.py +++ b/xrspatial/tests/test_surface_distance.py @@ -1,5 +1,7 @@ """Tests for xrspatial.surface_distance.""" +import warnings + import numpy as np import pytest import xarray as xr @@ -792,3 +794,415 @@ def test_error_message_mentions_grid_size(self): surface_distance(raster, elevation) with pytest.raises(MemoryError, match="dask"): surface_distance(raster, elevation) + + +# --------------------------------------------------------------------------- +# Backend coverage — every public function on every backend +# --------------------------------------------------------------------------- + +ALL_BACKENDS = ['numpy', 'dask+numpy', 'cupy', 'dask+cupy'] + +# Bounded vs iterative: with 6x6 chunks, max_distance=4 gives an overlap +# depth of 5, which fits inside a chunk, so the dask backends take the +# map_overlap branch. max_distance=inf forces the iterative tile sweep. +_BOUNDED = 4.0 + + +def _needs(backend): + if backend in ('dask+numpy', 'dask+cupy') and da is None: + pytest.skip("dask not installed") + if backend in ('cupy', 'dask+cupy') and not has_cuda_and_cupy(): + pytest.skip("cupy/cuda not available") + + +def _parity_scene(): + """A scene with relief and two distinguishable sources.""" + source = np.zeros((12, 12), dtype=np.float64) + source[2, 3] = 1.0 + source[9, 8] = 2.0 + elev = np.random.default_rng(42).uniform(0, 60, (12, 12)) + return source, elev + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +@pytest.mark.parametrize("backend", ALL_BACKENDS[1:]) +def test_bounded_matches_numpy(func, backend): + """Bounded (map_overlap) path matches the numpy baseline everywhere. + + surface_direction had no backend coverage at all before this test. + """ + _needs(backend) + source, elev = _parity_scene() + + expected = _compute(func(_make_raster(source), _make_raster(elev), + max_distance=_BOUNDED)) + actual = _compute(func(_make_raster(source, backend, chunks=(6, 6)), + _make_raster(elev, backend, chunks=(6, 6)), + max_distance=_BOUNDED)) + + np.testing.assert_allclose(actual, expected, rtol=1e-5, equal_nan=True) + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +@pytest.mark.parametrize("backend", ['dask+numpy', 'dask+cupy']) +def test_iterative_matches_numpy_all_modes(func, backend): + """Unbounded (iterative tile Dijkstra) path matches numpy. + + Regression test for the direction mode: ``_run_tile`` records global + source indices, so ``_finalize_direction`` has to build its pixel + grid in global coordinates too. Before the fix every tile but the + top-left one reported a bearing measured from its own origin. + """ + _needs(backend) + source, elev = _parity_scene() + + expected = _compute(func(_make_raster(source), _make_raster(elev))) + with pytest.warns(UserWarning, match="iterative"): + actual = _compute(func(_make_raster(source, backend, chunks=(6, 6)), + _make_raster(elev, backend, chunks=(6, 6)))) + + np.testing.assert_allclose(actual, expected, rtol=1e-4, equal_nan=True) + + +def test_iterative_direction_uses_global_pixel_indices(): + """Explicit check of the bearings the iterative path produces. + + A single source in the top-left chunk of a flat 6x6 raster gives an + unambiguous bearing for every pixel. The bug returned the top-left + chunk's bearing field tiled across the whole raster, so (3, 0) read + as 0 ("you are the source") instead of 360 (source due north). + """ + if da is None: + pytest.skip("dask not installed") + + source = np.zeros((6, 6), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.zeros((6, 6), dtype=np.float64) + + with pytest.warns(UserWarning, match="iterative"): + result = _compute(surface_direction( + _make_raster(source, 'dask+numpy', chunks=(3, 3)), + _make_raster(elev, 'dask+numpy', chunks=(3, 3)), + )) + + # Only the source itself reports 0. + assert result[0, 0] == 0.0 + assert np.count_nonzero(result == 0.0) == 1 + # Column 0 below the source: the source is due north. + np.testing.assert_allclose(result[1:, 0], 360.0, atol=1e-4) + # Row 0 right of the source: the source is due west. + np.testing.assert_allclose(result[0, 1:], 270.0, atol=1e-4) + # Perfect diagonal: north-west. + for i in range(1, 6): + assert result[i, i] == pytest.approx(315.0, abs=1e-4) + + +@pytest.mark.skipif(not has_cuda_and_cupy(), reason="cupy/cuda not available") +@pytest.mark.skipif(da is None, reason="dask not installed") +def test_dask_cupy_returns_dask_cupy_array(): + """dask+cupy input keeps both wrappers on the way out.""" + source = np.zeros((12, 12), dtype=np.float64) + source[5, 5] = 1.0 + elev = np.zeros((12, 12), dtype=np.float64) + + result = surface_distance( + _make_raster(source, 'dask+cupy', chunks=(6, 6)), + _make_raster(elev, 'dask+cupy', chunks=(6, 6)), + max_distance=_BOUNDED, + ) + assert isinstance(result.data, da.Array) + assert is_cupy_array(result.data.compute()) + + +@pytest.mark.skipif(not has_cuda_and_cupy(), reason="cupy/cuda not available") +def test_cupy_direction_returns_cupy_array(): + """The cupy direction branch converts its numpy result back to cupy.""" + source = np.zeros((5, 5), dtype=np.float64) + source[2, 2] = 1.0 + elev = np.zeros((5, 5), dtype=np.float64) + + result = surface_direction(_make_raster(source, 'cupy'), + _make_raster(elev, 'cupy')) + assert is_cupy_array(result.data) + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_target_values_all_backends(backend): + """target_values filtering behaves the same on every backend.""" + _needs(backend) + source = np.zeros((6, 6), dtype=np.float64) + source[0, 0] = 1.0 + source[5, 5] = 2.0 + elev = np.zeros((6, 6), dtype=np.float64) + + expected = _compute(surface_distance(_make_raster(source), + _make_raster(elev), + target_values=[2])) + raster = _make_raster(source, backend, chunks=(3, 3)) + elevation = _make_raster(elev, backend, chunks=(3, 3)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + actual = _compute(surface_distance(raster, elevation, + target_values=[2])) + + # Only the pixel holding value 2 is a source. + assert actual[5, 5] == 0.0 + assert actual[0, 0] == pytest.approx(5 * np.sqrt(2), abs=1e-4) + np.testing.assert_allclose(actual, expected, rtol=1e-5, equal_nan=True) + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_no_sources_all_nan_all_backends(backend): + """An all-zero source raster yields NaN everywhere on every backend.""" + _needs(backend) + zeros = np.zeros((4, 4), dtype=np.float64) + raster = _make_raster(zeros, backend, chunks=(2, 2)) + elevation = _make_raster(zeros, backend, chunks=(2, 2)) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + assert np.all(np.isnan(_compute(surface_distance(raster, elevation)))) + assert np.all(np.isnan(_compute(surface_allocation(raster, + elevation)))) + assert np.all(np.isnan(_compute(surface_direction(raster, + elevation)))) + + +# --------------------------------------------------------------------------- +# Edge cases — Inf and all-NaN elevation +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_inf_elevation_is_a_barrier(backend): + """Non-finite elevation blocks paths whether it is NaN or Inf.""" + _needs(backend) + source = np.zeros((1, 5), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.zeros((1, 5), dtype=np.float64) + elev[0, 2] = np.inf + + raster = _make_raster(source, backend, chunks=(1, 5)) + elevation = _make_raster(elev, backend, chunks=(1, 5)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + result = _compute(surface_distance(raster, elevation)) + + assert result[0, 0] == 0.0 + assert result[0, 1] == pytest.approx(1.0, abs=1e-5) + # The Inf cell and everything behind it are unreachable. + assert np.all(np.isnan(result[0, 2:])) + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_all_nan_elevation(backend): + """An entirely impassable elevation surface reaches nothing.""" + _needs(backend) + source = np.zeros((4, 4), dtype=np.float64) + source[1, 1] = 1.0 + elev = np.full((4, 4), np.nan, dtype=np.float64) + + raster = _make_raster(source, backend, chunks=(2, 2)) + elevation = _make_raster(elev, backend, chunks=(2, 2)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + result = _compute(surface_distance(raster, elevation)) + + assert np.all(np.isnan(result)) + + +# --------------------------------------------------------------------------- +# Geometric degeneracies +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_single_pixel_raster(backend): + """A 1x1 raster: the pixel is either its own source or unreachable.""" + _needs(backend) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + hit = _compute(surface_distance( + _make_raster(np.array([[1.0]]), backend, chunks=(1, 1)), + _make_raster(np.array([[0.0]]), backend, chunks=(1, 1)), + )) + miss = _compute(surface_distance( + _make_raster(np.array([[0.0]]), backend, chunks=(1, 1)), + _make_raster(np.array([[0.0]]), backend, chunks=(1, 1)), + )) + + assert hit.shape == (1, 1) + assert hit[0, 0] == 0.0 + assert np.isnan(miss[0, 0]) + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_column_strip(backend): + """An Nx1 column exercises the kernel with no horizontal neighbours.""" + _needs(backend) + source = np.zeros((5, 1), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.zeros((5, 1), dtype=np.float64) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + result = _compute(surface_distance( + _make_raster(source, backend, chunks=(2, 1)), + _make_raster(elev, backend, chunks=(2, 1)), + )) + + np.testing.assert_allclose(result.ravel(), [0.0, 1.0, 2.0, 3.0, 4.0], + rtol=1e-5) + + +# --------------------------------------------------------------------------- +# Parameter coverage +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_connectivity_4_all_backends(backend): + """4-connectivity is honoured by every backend, not just eager numpy.""" + _needs(backend) + source = np.zeros((6, 6), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.zeros((6, 6), dtype=np.float64) + + raster = _make_raster(source, backend, chunks=(3, 3)) + elevation = _make_raster(elev, backend, chunks=(3, 3)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + conn4 = _compute(surface_distance(raster, elevation, connectivity=4)) + conn8 = _compute(surface_distance(raster, elevation, connectivity=8)) + + # Without diagonals the corner pixel costs two cardinal steps. + assert conn4[1, 1] == pytest.approx(2.0, abs=1e-5) + assert conn8[1, 1] == pytest.approx(np.sqrt(2), abs=1e-5) + # 4-connectivity can never beat 8-connectivity. + assert np.all(conn4 >= conn8 - 1e-5) + + +@pytest.mark.parametrize("backend", ['cupy', 'dask+numpy', 'dask+cupy']) +def test_geodesic_unsupported_backends_raise(backend): + """Geodesic mode is numpy-only and says so per backend.""" + _needs(backend) + raster = _make_raster(np.zeros((4, 4)), backend, chunks=(2, 2)) + elevation = _make_raster(np.zeros((4, 4)), backend, chunks=(2, 2)) + + with pytest.raises(NotImplementedError, match="geodesic"): + surface_distance(raster, elevation, method='geodesic') + + +def test_invalid_dims(): + """A raster whose dims do not match x/y is rejected by name.""" + source = _make_raster(np.zeros((3, 3))).rename({'y': 'x', 'x': 'y'}) + elevation = _make_raster(np.zeros((3, 3))) + with pytest.raises(ValueError, match=r"raster.dims should be"): + surface_distance(source, elevation) + + +# --------------------------------------------------------------------------- +# Dataset support +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +def test_dataset_input(func): + """@supports_dataset maps the function over every data variable.""" + source = np.zeros((1, 4), dtype=np.float64) + source[0, 0] = 1.0 + elevation = _make_raster(np.zeros((1, 4))) + + ds = xr.Dataset( + {'a': _make_raster(source), 'b': _make_raster(source * 3.0)}, + attrs={'source': 'test'}, + ) + result = func(ds, elevation) + + assert isinstance(result, xr.Dataset) + assert list(result.data_vars) == ['a', 'b'] + assert result.attrs == {'source': 'test'} + for name in ('a', 'b'): + expected = func(ds[name], elevation) + np.testing.assert_allclose(_compute(result[name]), + _compute(expected), equal_nan=True) + + +def test_dataset_allocation_keeps_per_variable_values(): + """Allocation over a Dataset reports each variable's own source values.""" + source = np.zeros((1, 4), dtype=np.float64) + source[0, 0] = 1.0 + elevation = _make_raster(np.zeros((1, 4))) + + ds = xr.Dataset({'a': _make_raster(source), + 'b': _make_raster(source * 3.0)}) + result = surface_allocation(ds, elevation) + + np.testing.assert_allclose(_compute(result['a']), np.full((1, 4), 1.0)) + np.testing.assert_allclose(_compute(result['b']), np.full((1, 4), 3.0)) + + +# --------------------------------------------------------------------------- +# Metadata preservation +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_attrs_and_coords_preserved(func, backend): + """Input attrs, coords and dims survive the round trip. + + The module reads ``res`` off the input attrs to build its edge costs, + so losing them downstream would break any chained call. + """ + _needs(backend) + source = np.zeros((6, 6), dtype=np.float64) + source[1, 1] = 1.0 + elev = np.zeros((6, 6), dtype=np.float64) + + raster = _make_raster(source, backend, chunks=(3, 3)) + raster.attrs['crs'] = 4326 + raster.attrs['units'] = 'meters' + elevation = _make_raster(elev, backend, chunks=(3, 3)) + + result = func(raster, elevation, max_distance=_BOUNDED) + + assert result.attrs == raster.attrs + assert result.dims == raster.dims + assert result.shape == raster.shape + np.testing.assert_array_equal(result.y.values, raster.y.values) + np.testing.assert_array_equal(result.x.values, raster.x.values) + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +def test_custom_dim_names_preserved(func): + """lat/lon dims are not silently renamed to y/x.""" + source = np.zeros((4, 4), dtype=np.float64) + source[0, 0] = 1.0 + raster = xr.DataArray( + source, dims=['lat', 'lon'], + coords={'lat': np.arange(4, dtype=np.float64), + 'lon': np.arange(4, dtype=np.float64)}, + attrs={'res': (1.0, 1.0)}, + ) + elevation = raster.copy(data=np.zeros((4, 4))) + + result = func(raster, elevation, x='lon', y='lat') + + assert result.dims == ('lat', 'lon') + np.testing.assert_array_equal(result.lat.values, raster.lat.values) From e4b2da3bbdcb64d74688a020a1f482d04a0ca307 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:20:29 -0400 Subject: [PATCH 2/2] Record the surface_distance test-coverage sweep (#3725) --- .claude/sweep-test-coverage-state.csv | 1 + 1 file changed, 1 insertion(+) diff --git a/.claude/sweep-test-coverage-state.csv b/.claude/sweep-test-coverage-state.csv index c57ee8e5c..cef6da301 100644 --- a/.claude/sweep-test-coverage-state.csv +++ b/.claude/sweep-test-coverage-state.csv @@ -29,6 +29,7 @@ rasterize,2026-06-18,2614;3102;3105;3296;3383,HIGH,4,,"Pass 7 (2026-06-18, deep- reproject,2026-06-09,2618;3050;3100;3101;3141,MEDIUM,1,,"CI follow-up same day: first CI run of the threaded streaming branch hard-crashed macos-arm64 py3.14 (SIGABRT in numba call_cfunc, two ThreadPoolExecutor threads concurrently inside try_numba_transform/tmerc_inverse) -- the projection kernels are @njit(parallel=True) and numba's workqueue threading layer aborts on concurrent entry; filed source bug #3141. Test fix: threaded parity test now uses transform_precision=0 (per-thread pyproj Transformer, no numba), the NaN multi-tile test and 3-D xfail forced serial (max_memory=1) so the numba fast path stays covered without concurrent entry. windows-3.14 failure was fail-fast collateral (its suite fully passed). Pass 2026-06-09 (deep-sweep test-coverage): delta re-sweep one day after the 2026-06-08 pass; module modified today by #3077 (datum-probe warning silencing) and #3081 (merge output-size guard backend-aware) -- both landed WITH their own tests (TestDatumProbeNoProjWarning; TestSecurityGuards merge-guard trio incl. the monkeypatched in-memory raise), so the delta added no gap; the guard branching is is_dask-only, so cupy eager shares the tested numpy branch (no per-backend guard test needed). Found one MEDIUM Cat 1 gap every prior pass missed: the 5th dispatch branch of reproject() -- the streaming fallback (_reproject_streaming / _process_tile_batch / _parse_max_memory, taken when source >512MB and dask is not importable) -- had zero coverage anywhere; _parse_max_memory only runs on that branch so the existing max_memory kwarg tests never reached it. Filed #3101, added test_reproject_streaming_3101.py (15 tests: parity vs in-memory numpy for threaded / serial(max_memory=1) / single-tile / nearest+NaN, plus 10 _parse_max_memory unit cases). Probe surfaced source bug #3100: streaming assembly allocates a 2-D output buffer but 3-D sources yield (h,w,b) tiles -> ValueError broadcast in both assembly loops; pinned with strict xfail, source fix left to #3100 (test-only PR, source untouched). CPU-only path so no GPU tests needed (CUDA host; file ran 14 passed + 1 xfailed). LOW carried (documented, not fixed): reproject(name=) / merge(name=) override values untested (only merge name fallback covered); non-square-cellsize successful anisotropic run; dask.bag distributed branch of _reproject_streaming still unexercised (needs a live distributed client). || PREVIOUS: Pass 2026-06-08 (deep-sweep test-coverage): #3050 closes the one live gap found this pass. reproject()'s dask+cupy backend was parity-tested only with resampling='cubic' (TestCupyPyprojFallbackParity::test_projected_to_projected_dask_cupy_match); nearest/bilinear were covered on numpy (end-to-end) and eager cupy (parametrized test_projected_to_projected_numpy_cupy_match) but never on the dask+cupy chunk-assembly path. Parametrized that test over ['nearest','bilinear','cubic']; all 3 RUN+PASS on a CUDA host. Cat 4 MEDIUM (resampling-mode parameter coverage on the dask+cupy backend). Test-only, source untouched. Re-confirmed _merge.merge() has NO genuine cupy/dask+cupy backend (_merge_inmemory/_merge_dask use _merge_arrays_numpy + raster.values; _merge_arrays_cupy is imported but never dispatched = dead code, not a test gap) matching the prior pass's observation. reproject() otherwise saturated across all 4 backends, NaN/Inf/all-NaN, degenerate shapes, metadata, vertical, bounds_policy, integer nodata. LOW (documented, not filed): dask+cupy resampling-mode parity is the only per-mode-per-backend cell that had been missing. || PREVIOUS: Pass 2026-05-29: reproject already has a deep suite (369 tests in test_reproject.py + coverage/gate files) covering all 4 backends, NaN/Inf/all-NaN/all-Inf, 1x1/2x2, metadata, vertical shift, bounds_policy x backends, integer nodata x backends. Gaps found: Cat 3 HIGH single-row (1xN) and single-col (Nx1) strip rasters never tested (hit size<2 branch of _validate_regular_axis + degenerate resampling axis); Cat 3 MEDIUM constant-value/zero-gradient raster never reprojected. Added TestDegenerateShapeReproject (12 tests): 1xN+Nx1 strips x numpy/dask/cupy/dask+cupy, constant raster numpy value-preservation + cross-backend parity. All 12 executed and passed on a CUDA host. Test-only, no source change (#2618). LOW (documented only): _merge._merge_arrays_cupy imported but never called by merge() (host-bounces via _merge_arrays_numpy) - dead-code source observation not a test gap; non-square cellsize reproject only covered via resolution-tuple validation errors not a successful anisotropic run." resample,2026-05-29,2547;2615,HIGH,1;2;3;5,,"Pass 2 (2026-05-29): added test_resample_cupy_agg_fallback_2615.py (6 tests, all passing on CUDA host). Closes Cat 1 MEDIUM backend-coverage gap: the cupy eager aggregate CPU fallback for average/min/max at a NON-integer downsample factor (_run_cupy fy==int(fy) branch in resample.py ~L957-973) was never exercised; existing TestCuPyParity used 12x12 scale 0.5 (integer factor 2 -> GPU reshape path) and only median/mode hit the host fallback. New tests use 10x10 scale 0.3 (factor 3.33) for average/min/max parity vs numpy plus a NaN-masked variant. Issue #2615. Module is otherwise very thoroughly covered (test_resample.py + 3 supplementary files); no remaining HIGH gaps found. Pass 1 (2026-05-27): added test_resample_coverage_2026_05_27.py with 70 tests (68 passing, 2 skipped). Closes Cat 3 HIGH Nx1 single-column gap across numpy/cupy/dask+numpy/dask+cupy x 8 methods (nearest/bilinear/cubic/average/min/max/median/mode) plus Nx1 upsample-nearest parity and Nx1 cross-backend aggregate parity. Closes Cat 2 MEDIUM NaN-parity gap on cupy and dask+cupy (existing TestCuPyParity/TestDaskCuPyParity used random data without NaN; the weight-mask gate and spline-prepad had no GPU NaN coverage). Closes Cat 3 MEDIUM all-equal-value raster across 8 methods (downsample) and 3 interp methods (upsample) plus a constant-with-NaN aggregate variant. Closes Cat 5 MEDIUM non-default dim-name propagation: lat/lon, latitude/longitude, and (channel, lat, lon) 3D round-trip without being renamed to y/x; per-dim attrs (units) preserved. Closes Cat 3 MEDIUM empty-raster behaviour pin: 0-row and 0-col rasters raise (currently IndexError) -- contract covered. Filed source-bug issue #2547: cubic on dask backends fails for Nx1 / arrays smaller than depth=16; the 2 skipped tests in this file gate on that fix landing. Source untouched." slope,2026-05-29,2697,MEDIUM,3,,"PR #2703: added degenerate-shape tests (1x1/1xN/Nx1) for all 4 planar backends + geodesic; no live bug, pins all-NaN+shape contract. CUDA host: cupy/dask+cupy ran. Backend/NaN/param/metadata coverage already complete." +surface_distance,2026-08-16,3725,HIGH,1;2;3;4;5,87,"Deep-sweep 2026-08-16 test-coverage on a CUDA host (cupy + dask+cupy RAN, not skipped). Before: 37 tests, branch_cov 80 (570/695). After: 106 tests, branch_cov 87. Cat 1 HIGH: _surface_distance_dask_cupy (:1232-1289) had ZERO coverage (neither the bounded map_overlap branch nor the unbounded convert-to-dask+numpy branch), and surface_direction was only ever called on eager numpy -- no cupy, no dask, no dask+cupy. That second hole was hiding a real SOURCE bug (#3713, filed by the sibling documentation sweep, no PR then): _run_tile records GLOBAL src_row/src_col while _finalize_direction built its pixel grid from BLOCK-LOCAL np.arange(H)/np.arange(W), so the dask iterative path returned the top-left tile's bearing field tiled across the whole raster (flat 6x6 chunked (3,3): 360 deg max error, four pixels reporting 0 'you are the source'). Fixed in the same PR by threading row_offset/col_offset through _extract_output -> _finalize_direction from _assemble_sd; eager numpy and bounded map_overlap keep offset 0 because they seed local indices. Fail-before/pass-after verified. Also covered: cupy target_values seeding (:560-561) and cupy no-source early return (:571-573); dask iterative + target_values (:698-701, :956-959); connectivity=4 off the eager path (:919-932). Cat 2 MEDIUM: Inf elevation and all-NaN elevation untested -- both behave correctly (non-finite = barrier) on all 4 backends, coverage gap not a bug. Cat 3 MEDIUM: no 1x1 and no Nx1 column-strip tests (only 1xN rows existed); both correct on all 4 backends. Cat 4 MEDIUM: all three public funcs carry @supports_dataset and none was ever called with a Dataset (same gap as the proximity trio) -- works correctly; error paths untested (raster.dims rejection :1306, geodesic NotImplementedError per backend :1371/:1379/:1401). Cat 5 HIGH: nothing asserted attrs/coords/dims preservation even though _compute reads res off attrs for its edge costs -- all 4 backends preserve them correctly. LOW (documented, not fixed): non-square cells (res_x != res_y) untested but correct (verified N=2.0/E=1.0 on a 2x1 grid); _available_memory_bytes psutil + 2GB-default fallbacks (:107-115), _available_gpu_memory_bytes failure path (:129-130), _check_gpu_memory GPU-branch raise (:156,:159), _precompute_dd_grid memory guard (:423-426) all unexercised; _dijkstra_geodesic visited/over-budget/non-finite branches (:295,:299, :313) unreached. Pre-existing flake8 F841 at test file line 236 (sd assigned but unused in test_allocation_consistency) left alone per surgical-change rule. Post-fix branch_cov of 87 is still measured with NUMBA_DISABLE_JIT=1, under which the cupy tests fail, so true with-GPU coverage is higher; full file is 106/106 GREEN with JIT on." templates,2026-06-30,3580,MEDIUM,1;4,,"Deep-sweep test-coverage re-run on a CUDA host (cuda available). Module is already heavily tested (281 tests): 4-backend matrix (numpy/dask+numpy/cupy/dask+cupy) + dask alias + bad-backend all green; preserve area/shape across backends; single-pixel + Nx1/1xN strips; cell-cap and chunk-count guards; padding/tiling helpers; country/region/city resolution + aliases; CF metadata + no-pyproj fallback. Cat 2 N/A (procedural generator, no raster input). Found two MEDIUM parameter-coverage gaps, no source bug. Cat 4: chunks=tuple only exercised via internal _estimate_n_chunks, never end-to-end through from_template (int/'auto' were); the dask chunk-count guard message on the explicit height/width path was untested (only the resolution-path message and the eager cell-cap message named the knob). Added test-only test_chunks_tuple_through_public_api and test_explicit_shape_chunk_count_message_names_height_width; both RAN and PASSED on the CUDA host (281 passed). LOW (documented, no test): non-NaN fill is only asserted on eager numpy (fill=0); probed live and works on all 4 backends but cross-backend fill value parity is not asserted. PR #3580 opened with the two tests." viewshed,2026-05-29,2693,HIGH,1;2;5,,"Pass 1 (2026-05-29): added 4 new test groups to test_viewshed.py (13 new tests + 1 xfail, all passing/xfailing on a CUDA+RTX host). Closes Cat 1 HIGH backend-coverage gap: the dask+cupy dispatch path in _viewshed_dask (Tier B) and _viewshed_windowed (max_distance) was registered but never invoked by any test -- added test_viewshed_dask_cupy_flat (analytical-angle parity, atol 0.03) and test_viewshed_dask_cupy_max_distance (windowed GPU run; observer cell 180, corners INVISIBLE). Both use non-zero flat terrain (1.3) because the RTX mesh builder rejects an all-zero raster (#1378). Closes Cat 5 HIGH metadata-preservation gap: only the numpy test_viewshed called general_output_checks; the cupy/dask/dask+cupy and max_distance paths never asserted attrs/coords/dims/array-type preservation. Added parametrised test_viewshed_metadata_preserved over {numpy,cupy,dask+numpy,dask+cupy} x {full, max_distance=2.0}: asserts attrs==, dims==, shape==, x/y coords allclose; runs general_output_checks (full type parity) for all backends except dask+cupy. Closes Cat 2 HIGH NaN-input gap and surfaced source bug #2693: viewshed on a numpy raster crashes with ValueError 'node not found' from _delete_from_tree when a NaN cell sits at certain positions (e.g. (2,4) in a 5x5 with observer at (2,2)), while NaN at (1,1)/(0,0)/(4,4) runs fine. Added test_viewshed_nan_input_supported_positions (parametrised working positions, asserts observer=180 and NaN cell is INVISIBLE/NaN) plus test_viewshed_nan_input_crashing_position (xfail strict, raises, links #2693). Noted but NOT fixed (source change out of scope for test sweep): the dask+cupy backend does not preserve the cupy backing -- _viewshed_dask computes then rewraps via da.from_array(result_np), so the output computes to numpy not cupy; general_output_checks is skipped for dask+cupy for that reason (candidate for the metadata/backend-parity sweep). LOW (documented only): non-square cell sizes; 1x1 and 1xN geometry covered behaviourally by probing (run without error). Test-only PR; viewshed.py untouched." visibility,2026-06-10,3192,HIGH,1;2;4,,"cupy cumulative_viewshed/visibility_frequency broken (numpy count + cupy viewshed) -> issue #3192 (dup #3193), fix in flight in #3205 with its own cupy parity tests, xfail pins dropped to avoid an XPASS race; added cupy _extract_transect+line_of_sight parity, NaN LOS, Fresnel-blocked branch; dask+metadata already covered"