From aaa11607a11485fafe15a355aaf6a5e8cb249e92 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:13:25 -0400 Subject: [PATCH] Stop surface_distance outputs adopting the dask graph token as .name (#3708) surface_distance, surface_allocation and surface_direction built their result with xr.DataArray(..., coords=, dims=, attrs=) and left name unset. DataArray.__init__ then falls back to getattr(data, "name"), which on a dask array is the graph key, so the dask backends returned '_trim-' (bounded map_overlap route), 'xrspatial.surface_*-' (unbounded iterative route) or 'asarray-' (dask+cupy unbounded) while numpy and cupy returned None. The divergence is user-visible through .to_dataset(): the numpy result raises "unable to convert unnamed DataArray", the dask result silently creates a variable named after the hash, and the hash moves when chunking or max_distance changes. Route the three functions through a shared _wrap_result() helper that resets .name to None after construction, matching the proximity/allocation/direction trio this module mirrors. Same bug class as cost_distance #3344 and pathfinding #3652. Tests cover .name parity over 4 backends x 3 functions x bounded/unbounded max_distance, plus an attrs/coords/dims/dtype preservation guard. The 12 dask cases fail without the fix. Also records the metadata sweep result for this module in .claude/sweep-metadata-state.csv. --- .claude/sweep-metadata-state.csv | 1 + xrspatial/surface_distance.py | 40 ++++++++------- xrspatial/tests/test_surface_distance.py | 65 ++++++++++++++++++++++++ 3 files changed, 88 insertions(+), 18 deletions(-) diff --git a/.claude/sweep-metadata-state.csv b/.claude/sweep-metadata-state.csv index 1872a9a6e..5c5455cdc 100644 --- a/.claude/sweep-metadata-state.csv +++ b/.claude/sweep-metadata-state.csv @@ -17,6 +17,7 @@ proximity,2026-05-29,2723,MEDIUM,4;5,"Audited 2026-05-29 (agent-a61dbadc2452a200 rasterize,2026-06-09,3087,MEDIUM,1,GeoDataFrame .crs dropped on no-like path (Cat 1); fixed via #3087 emitting attrs crs/crs_wkt when output has no CRS. like-path attrs/coords/dims/nodata verified live on all 4 backends (CUDA available); Cats 2-5 clean. reproject,2026-06-12,3262,MEDIUM,4,"Re-audited 2026-06-12 (agent-ae420c90e50a23c5c worktree, branch deep-sweep-metadata-reproject-2026-06-12). CUDA available; all 4 backends (numpy/cupy/dask+numpy/dask+cupy) run live end-to-end for reproject() and merge(). Cat 1 attrs (crs/nodata/res/transform/_FillValue/nodatavals refreshed or carried, crs_wkt dropped), Cat 2 coords (pixel-center verified numerically, scalar time + band coord carry, float64), Cat 3 dims (lat/lon names, band-first (band,y,x) round-trip), Cat 4 int16 sentinel parity, and Cat 5 cross-backend attr parity all identical across the 4 backends for reproject(); vertical_crs=4979/vertical_datum verified on numpy + dask; geoid_height_raster carries input attrs per its documented contract. NEW MEDIUM finding #3262 (Cat 4): merge() hardcoded float64 output on every path (_merge_inmemory, _merge_dask template/meta, empty-chunk fills) while reproject() round-trips integer dtypes on all 5 paths (#2185/#2505/#3093/#3096); undocumented and unpinned by tests, so an int16/uint8 mosaic silently promoted (8x memory for uint8, GeoTIFF round-trip changes file dtype). Fix on this branch: shared-integer-dtype inputs now cast back via the reproject round/clip/cast convention (_cast_merged_dtype), output nodata resolved with _detect_nodata dtype hint (NaN->sentinel swap per #2185, explicit out-of-range raises per #2572), dask template/meta + empty-chunk fills use the output dtype (#3096 trap), docstring documents the rule; mixed/float inputs keep float64. 13 new tests in TestMergeIntegerDtype (eager/dask/dask-empty-chunk/cupy, sentinel defaults, mean rounding, out-of-range raise); full reproject suite 514 passed. LOW (documented, not fixed): reproject() docstring says dask inputs are fully lazy but the dask+cupy VRAM-fitting fast path returns an eager cupy array (codified in tests; doc nit). Prior LOW from 2026-06-09 (geoid_height ndarray return for DataArray input) unchanged." resample,2026-05-27,2542,MEDIUM,2;4;5,"Audited 2026-05-27 (agent-a8135a6a246ecb93c worktree, branch deep-sweep-metadata-resample-2026-05-27). Cat 2 MEDIUM + Cat 4 MEDIUM + Cat 5 MEDIUM all rolled into issue #2542. (a) 2D non-identity path dropped scalar non-dim coords like rioxarrays spatial_ref and squeezed time/band selectors; identity path (scale==1.0, agg.copy()) and 3D path (per-band xr.concat) preserved them, so the bug was path-inconsistent (Cat 5). (b) _resolve_nodata reads attrs[nodata] as a fallback sentinel but the output post-processing only refreshed _FillValue and nodatavals, leaving attrs[nodata]=-9999 alongside data that was now NaN. Fix in resample(): refresh attrs[nodata] to NaN whenever the input had it, and carry across zero-dim non-dim coords on the 2D non-identity path. 7 new tests in TestMetadataPropagation cover nodata-attr refresh, spatial_ref/scalar coord carry, identity-vs-downsample coord parity, and the explicit choice to drop spatially-shaped extra coords. 4-backend (numpy/cupy/dask+numpy/dask+cupy) parity verified for spatial_ref carry; nodata-attr refresh verified on numpy/cupy/dask+numpy (dask+cupy non-NaN nodata masking hits a pre-existing xarray xr.where + cupy.astype quirk unrelated to this audit). Full resample test suite (175 passed) clean." +surface_distance,2026-08-16,3708,MEDIUM,5,"Audited 2026-08-16 (agent-ab3afc98fdd524ae1 worktree, branch deep-sweep-metadata-surface_distance-2026-08-16). CUDA available; all 4 backends (numpy/cupy/dask+numpy/dask+cupy) run live end-to-end for surface_distance/surface_allocation/surface_direction, across both dask routes (bounded map_overlap and unbounded iterative tile Dijkstra). Cat 1 attrs, Cat 2 coords, Cat 3 dims all clean: the three public functions re-emit coords=raster.coords, dims=raster.dims, attrs=raster.attrs, so res/crs/transform/nodatavals/_FillValue, extra scalar coords (spatial_ref), coord values and dtypes, and custom dim names via x=/y= (lat/lon) all survive identically on every backend. Output dtype is float32 on all four backends, matching the docstring. NEW MEDIUM finding #3708 (Cat 5): the three functions left name unset at the xr.DataArray constructor, so DataArray.__init__ fell back to getattr(data, 'name') and the dask backends adopted the graph key as .name -- '_trim-' from the bounded map_overlap route, 'xrspatial.surface_*-' from the unbounded iterative route, and 'asarray-' from the dask+cupy unbounded route -- while numpy and cupy returned None. User-visible via .to_dataset(): numpy raises 'unable to convert unnamed DataArray', dask silently creates a variable named '_trim-fe8156...' whose hash moves with chunking and max_distance. Recurring bug class: slope #2837, aspect #2841, focal #2733, viewshed #2743, zonal #2611, cost_distance #3344, pathfinding #3652. Fix in PR #3709: a shared _wrap_result() helper that sets result.name = None after construction, matching the proximity/allocation/direction trio this module mirrors. 36 new tests (name parity across 4 backends x 3 functions x bounded/unbounded, plus an attrs/coords/dims/dtype preservation guard); the 12 dask cases fail without the fix. Full suite 73 passed. LOW, documented not fixed: attrs are copied verbatim so a user-supplied nodatavals=(-9999,) / _FillValue=-9999 stays on an output that actually uses NaN as its nodata sentinel, and a 'units' attr carries over onto surface_direction output measured in degrees -- both are the library-wide attrs=raster.attrs convention (proximity and cost_distance behave identically), not surface_distance-specific. Also noted out of scope: the module's _dask_task_name_kwargs task names are overwritten by map_overlap's '_trim-' prefix on the bounded route, so the #3256 task-naming work does not reach that path. No CRITICAL or HIGH findings." viewshed,2026-05-29,2743,MEDIUM,4;5,output .name differed across backends (None/viewshed/dask-token) and dtype float32 on GPU vs float64 on CPU; added name= param and forced float64 on all backends; attrs/coords/dims already preserved visibility,2026-06-10,3193,HIGH,5,"cupy backend crash in cumulative_viewshed/visibility_frequency (count np vs cupy add) -> no result/metadata emitted; fixed by cupy count branch + cupy tests. numpy/dask preserve coords/dims/attrs incl crs; visibility_frequency keeps attrs through astype/divide. line_of_sight Dataset drops crs/transform (LOW, transect not raster, documented only)." zonal,2026-05-29,2611,MEDIUM,5,"Audited 2026-05-29 (agent-ae8d8b65cc3a5c40a worktree, branch deep-sweep-metadata-zonal-2026-05-29). CUDA available; all 4 backends (numpy/cupy/dask+numpy/dask+cupy) run live. 5 DataArray-returning functions checked end-to-end: apply, regions, hypsometric_integral, trim, crop. attrs (res/crs/transform/nodatavals), dims, and coords preserved correctly on all 4 backends for every function; trim/crop slice coords with no half-pixel drift. stats() and crosstab() return DataFrames by design so Cat 1-3 DataArray checks N/A. NEW MEDIUM finding #2611 (Cat 5): apply() never set output .name, so numpy/cupy returned None while dask+numpy/dask+cupy inherited a non-deterministic internal dask task name (e.g. _chunk_fn-). regions/hypsometric_integral/trim/crop all set deterministic names; apply was the outlier. Fix in PR #2611/#2622: add name param (default None) and assign result.name after DataArray construction (setting name= at construction does not override the dask graph name). New parametrized test test_apply_name_consistent_across_backends covers default-None and explicit-name on all 4 backends. Full zonal suite 213 passed. No other CRITICAL/HIGH/MEDIUM findings; no LOW findings to document." diff --git a/xrspatial/surface_distance.py b/xrspatial/surface_distance.py index c91aeb13e..f037caebc 100644 --- a/xrspatial/surface_distance.py +++ b/xrspatial/surface_distance.py @@ -1415,6 +1415,25 @@ def _compute(raster, elevation, x, y, target_values, max_distance, # --------------------------------------------------------------------------- +def _wrap_result(result_data, raster): + """Wrap raw output in a DataArray carrying the input's spatial metadata. + + The name is reset after construction: on dask backends + ``xr.DataArray(..., name=None)`` adopts the dask array's graph key + (``_trim-`` from map_overlap, ``xrspatial.surface_*-`` from + the iterative path) as ``.name``, while numpy and cupy return None. + Issue #3708; same fix as cost_distance #3344 and pathfinding #3652. + """ + result = xr.DataArray( + result_data, + coords=raster.coords, + dims=raster.dims, + attrs=raster.attrs, + ) + result.name = None + return result + + @supports_dataset def surface_distance( raster: xr.DataArray, @@ -1468,12 +1487,7 @@ def surface_distance( raster, elevation, x, y, target_values, max_distance, connectivity, method, DISTANCE, ) - return xr.DataArray( - result_data, - coords=raster.coords, - dims=raster.dims, - attrs=raster.attrs, - ) + return _wrap_result(result_data, raster) @supports_dataset @@ -1510,12 +1524,7 @@ def surface_allocation( raster, elevation, x, y, target_values, max_distance, connectivity, method, ALLOCATION, ) - return xr.DataArray( - result_data, - coords=raster.coords, - dims=raster.dims, - attrs=raster.attrs, - ) + return _wrap_result(result_data, raster) @supports_dataset @@ -1553,9 +1562,4 @@ def surface_direction( raster, elevation, x, y, target_values, max_distance, connectivity, method, DIRECTION, ) - return xr.DataArray( - result_data, - coords=raster.coords, - dims=raster.dims, - attrs=raster.attrs, - ) + return _wrap_result(result_data, raster) diff --git a/xrspatial/tests/test_surface_distance.py b/xrspatial/tests/test_surface_distance.py index 74ae62fbd..ad1c24461 100644 --- a/xrspatial/tests/test_surface_distance.py +++ b/xrspatial/tests/test_surface_distance.py @@ -675,6 +675,71 @@ def test_geodesic_basic(): assert 100000 < sd[pos] < 130000 # roughly 100-130 km +# --------------------------------------------------------------------------- +# Metadata propagation (issue #3708) +# --------------------------------------------------------------------------- + + +def _metadata_backends(): + backends = ['numpy'] + if da is not None: + backends.append('dask+numpy') + if has_cuda_and_cupy(): + backends.append('cupy') + if da is not None: + backends.append('dask+cupy') + return backends + + +@pytest.mark.parametrize("func", [surface_distance, surface_allocation, + surface_direction]) +@pytest.mark.parametrize("max_distance", [3.0, np.inf], + ids=['bounded', 'unbounded']) +@pytest.mark.parametrize("backend", _metadata_backends()) +def test_output_name_consistent_across_backends(backend, max_distance, func): + """Outputs must not adopt the dask graph token as .name. + + Without the post-construction reset the dask backends returned + '_trim-' (bounded map_overlap route), + 'xrspatial.surface_*-' (unbounded iterative route) or + 'asarray-' (dask+cupy unbounded), while numpy and cupy returned + None. Same bug class as cost_distance #3344 and pathfinding #3652. + """ + source = np.zeros((6, 6), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.arange(36, dtype=np.float64).reshape(6, 6) * 0.1 + raster = _make_raster(source, backend=backend) + elevation = _make_raster(elev, backend=backend) + + result = func(raster, elevation, max_distance=max_distance) + assert result.name is None + + +@pytest.mark.parametrize("func", [surface_distance, surface_allocation, + surface_direction]) +@pytest.mark.parametrize("backend", _metadata_backends()) +def test_output_preserves_attrs_coords_dims(backend, func): + """attrs, coords and dims come through unchanged on every backend.""" + source = np.zeros((6, 6), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.arange(36, dtype=np.float64).reshape(6, 6) * 0.1 + raster = _make_raster(source, backend=backend) + elevation = _make_raster(elev, backend=backend) + raster.attrs.update({'crs': 3857, 'nodatavals': (-9999.0,), + 'transform': (1.0, 0.0, 0.0, 0.0, -1.0, 0.0)}) + raster = raster.assign_coords(spatial_ref=0) + + result = func(raster, elevation, max_distance=3.0) + + assert result.dims == raster.dims + assert result.attrs == raster.attrs + assert set(result.coords) == set(raster.coords) + for name in raster.coords: + np.testing.assert_array_equal(result.coords[name].values, + raster.coords[name].values) + assert result.dtype == np.float32 + + # --------------------------------------------------------------------------- # Memory guard # ---------------------------------------------------------------------------