Stop surface_distance outputs adopting the dask graph token as .name (#3708)

.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-<hash>' from the bounded map_overlap route, 'xrspatial.surface_*-<hash>' from the unbounded iterative route, and 'asarray-<hash>' 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-<hash>). 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."

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-<hash>`` from map_overlap, ``xrspatial.surface_*-<hash>`` 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)

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-<hash>' (bounded map_overlap route),
+    'xrspatial.surface_*-<hash>' (unbounded iterative route) or
+    'asarray-<hash>' (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
 # ---------------------------------------------------------------------------