Fix surface_direction bearings on north-up rasters and on chunked dask input

.claude/sweep-accuracy-state.csv

@@ -35,6 +35,7 @@ reproject,2026-07-27,3697,HIGH,4;6,"2026-07-27 sweep: filed #3697 (HIGH, cat 4;6
 resample,2026-05-29,2610,HIGH,3;5,"dask interp (nearest/bilinear) overlap depth=1 too small on downsample; block-centered source coord landed past chunk, map_coordinates clamped to edge -> wrong seam rows. Fixed PR #2627 via per-axis _downsample_radius. cupy+dask+cupy verified."
 sieve,2026-04-13T12:00:00Z,,,,Union-find CCL correct. NaN excluded from labeling. All backends funnel through _sieve_numpy.
 sky_view_factor,2026-07-03,3626,MEDIUM,4,"Re-audit post-#1407 fix (ground distance OK). MEDIUM #3626: ray azimuths uniform in cell-INDEX angle not ground azimuth, so anisotropic cellsizes bias the azimuthal quadrature (45-deg ramp, true SVF 0.75: 0.728/0.773/0.693 at 1:1/1:2/2:1 aspect; +0.105 pure weighting error at 4:1); fixed by aiming rays in ground space (cos(phi)/cx, sin(phi)/cy normalized), square-cell results bit-identical. Cats 1-3,5 clean: float64 throughout, NaN center->NaN out, flat=1.0 exactly all 4 backends, flip/rot90/translate invariants at machine eps, numpy==dask==cupy==dask+cupy (max 2.2e-16) incl NaN holes straddling chunk boundaries, map_overlap depth==max_radius==stencil radius, depth>chunksize works (dask 2025.7). Cat6 skipped: richdem-unavailable rvt-unavailable grass-unavailable. LOWs (not fixed): NaN neighbor silently truncates ray (horizon beyond nodata ignored, undocumented but keeps numpy/dask edge semantics consistent); nearest-cell ray rounding gives ~-2.9% jitter bias on 45-deg slopes even square cells (inherent to integer-cell marching, reference tools interpolate along ray). For test-coverage sweep: cross-backend tests had no NaN-input and no anisotropic-cellsize cases (added aniso numpy/dask/cupy parity in #3626 PR); for style sweep: unused import _boundary_to_dask (F401)."
+surface_distance,2026-08-16,3719;3721;3733,HIGH,2;4;5,"3 HIGH found, all reproduced on this host; CUDA available so cupy and dask+cupy ran for real. #3719 (PR): surface_direction built bearings from row/col indices scaled by abs(cellsize), so a north-up (descending-y) raster came back mirrored north/south and disagreed with proximity.direction on identical input - same bug class as #2896; fix passes signed cell sizes (every edge-cost site squares or abs()es them). Same PR: the dask iterative path fed _finalize_direction global src indices against a block-local np.arange grid, so every chunk but (0,0) reported standalone-raster bearings (58/80 pixels wrong, max diff 360) - default path since max_distance defaults to inf; rebasing to block-local indices then exposed a third defect, the no-source mask src_row<0 firing on legitimately negative rebased indices, now keyed on dist alone. #3721 (PR): cupy _surface_distance_cupy caps parallel relaxation at H+W passes, but the bound is the longest shortest-path hop count (worst case H*W); on a serpentine NaN maze the GPU reported 99 of 136 reachable pixels as NaN at 16x16 and 247 of 300 at 24x24, silently, with no warning. Cat1 clean (float64 Dijkstra, float32 output documented). Cat2 clean - checked the int(inf) rejection-guard class in _run_tile's int(ssr[...]) and it is unreachable (finite dist always implies an assigned src). Cat3 clean - map_overlap depth max_distance/min_cellsize is provably sufficient since each pixel-radius step costs at least min_cellsize horizontally; corner index picks in _can_skip_sd/_compute_seeds_sd verified. Cat4: geodesic haversine radius 6378137 matches proximity.py's convention (clean), but geodesic surface_direction scales index offsets by degrees with no cos(lat) correction, so bearings tilt away from the equator - MEDIUM, deferred to #3733 because the fix needs an azimuth convention decision. Cat6 reference-unavailable (no terrain surface-distance oracle installed). Explicitly NOT filed: the dask iterative sweep cap max(n_tile_y,n_tile_x)+10 survived serpentine mazes at 4 chunk sizes; cupy-vs-numpy allocation differences on flat terrain are equidistant tie-breaks, not errors. Test-coverage gaps for the test-coverage sweep: surface_direction had no dask, no cupy, and no descending-y coverage at all before this PR, and no test uses a north-up raster; connectivity=4 and geodesic are numpy-only. Doc gap: surface_allocation does not document a tie-break rule the way proximity.allocation does, and the geodesic NotImplementedError on cupy/dask is absent from the docstrings. Pre-existing flake8 F841 at test_surface_distance.py:234 left alone."
 terrain,2026-04-10T12:00:00Z,,,,Perlin/Worley/ridged noise correct. Dask chunk boundaries produce bit-identical results. No precision issues.
 terrain_metrics,2026-04-30,,LOW,2;5,"LOW: Inf input not rejected, propagates as Inf (consistent across backends but undocumented). LOW: dask+cupy non-nan boundary path double-pads (wasted compute, central output values still correct). No CRIT/HIGH; tests cover NaN propagation, all 4 backends, all 4 boundary modes, dtype acceptance."
 viewshed,2026-05-29,2691,HIGH,3;5,max_distance window sized from coarser axis clipped cells on anisotropic rasters (PR #2702). LOW unfixed: distance_sweep ring radius same max(res) pattern but max_distance arg always None; _calculate_event_row_col line 880 abs(x>1) precedence bug is a broken guard only. cuda+rtx paths validated.

xrspatial/surface_distance.py

@@ -331,6 +331,27 @@ def _dijkstra_geodesic(elev_data, height, width, max_distance,
 # ---------------------------------------------------------------------------
 
 
+def _coord_step_sign(raster, dim):
+    """Sign of the step along *dim* as the array index increases.
+
+    Returns -1.0 for a descending coordinate (the north-up convention on
+    the y axis), +1.0 otherwise.  Rasters without a usable coordinate on
+    *dim* fall back to +1.0, which is index order.
+
+    The axis is assumed monotonic, the same assumption ``calc_res``
+    already makes when it divides the coordinate span by ``n - 1``.
+    """
+    if dim not in raster.coords:
+        return 1.0
+    values = np.asarray(raster.coords[dim].values)
+    if values.ndim != 1 or values.size < 2:
+        return 1.0
+    if not np.issubdtype(values.dtype, np.number):
+        return 1.0
+    step = float(values[-1]) - float(values[0])
+    return -1.0 if step < 0 else 1.0
+
+
 def _init_arrays(H, W):
     """Create and initialize output arrays for the Dijkstra kernel."""
     dist = np.full((H, W), np.inf, dtype=np.float64)
@@ -360,7 +381,11 @@ def _finalize_direction(src_row, src_col, dist, cellsize_x, cellsize_y,
                         max_distance):
     """Compute compass bearing from each pixel to its allocated source.
 
-    Uses pixel index differences scaled by cell size.
+    Uses pixel index differences scaled by cell size.  *cellsize_x* and
+    *cellsize_y* are signed: their sign is the direction the matching
+    coordinate runs as the index rises, so a north-up raster (descending
+    y) gets a negative *cellsize_y* and the bearings come out in true
+    compass space rather than mirrored in index space.
     """
     H, W = dist.shape
     row_idx, col_idx = np.meshgrid(np.arange(H), np.arange(W), indexing='ij')
@@ -374,8 +399,14 @@ def _finalize_direction(src_row, src_col, dist, cellsize_x, cellsize_y,
         np.zeros((H, W), dtype=np.float64), dy,
     )
 
-    # Mask unreachable and no-source pixels
-    mask = np.isinf(dist) | (dist > max_distance) | (src_row < 0)
+    # Mask unreachable pixels.  Every pixel with a finite in-budget
+    # distance was assigned a source when that distance was written, so
+    # `dist` alone identifies them, and the direction output ends up with
+    # the same NaN mask as the distance and allocation outputs.  Testing
+    # `src_row < 0` instead would be wrong on the tiled dask path, where
+    # a source in a chunk above or to the left rebases to a legitimately
+    # negative block-local index.
+    mask = np.isinf(dist) | (dist > max_distance)
     result[mask] = np.nan
     return result
 
@@ -610,8 +641,10 @@ def _surface_distance_cupy(source_data, elev_data, cellsize_x, cellsize_y,
         dy_np = cp.asnumpy(dy_coord)
         zeros = np.zeros((H, W), dtype=np.float64)
         result = _vectorized_calc_direction(zeros, dx_np, zeros, dy_np)
-        mask_np = cp.asnumpy(cp.isinf(dist) | (dist > max_distance)
-                             | (srow < 0))
+        # Same rule as _finalize_direction: a finite in-budget distance
+        # always came with an assigned source, so `dist` alone marks the
+        # unreachable pixels and all three outputs share one NaN mask.
+        mask_np = cp.asnumpy(cp.isinf(dist) | (dist > max_distance))
         result[mask_np] = np.nan
         out = cp.asarray(result)
 
@@ -1174,6 +1207,15 @@ def _tile_fn(source_block, elev_block, block_info=None):
             max_distance, dy, dx, dd, row_offset, col_offset,
         )
 
+        # _run_tile tracks sources in global raster indices so seeds stay
+        # comparable across tiles, but _finalize_direction measures the
+        # offset against this block's own np.arange grid.  Rebase to
+        # block-local indices so the two sides agree; a source outside
+        # this block lands outside [0, h) x [0, w), which is exactly the
+        # relative offset the bearing needs.
+        srow = srow - row_offset
+        scol = scol - col_offset
+
         return _extract_output(dist, alloc_arr, srow, scol,
                                cellsize_x, cellsize_y, max_distance, mode)
 
@@ -1315,6 +1357,13 @@ def _compute(raster, elevation, x, y, target_values, max_distance,
     cellsize_x = abs(float(cellsize_x))
     cellsize_y = abs(float(cellsize_y))
 
+    # The DIRECTION output is a compass bearing, so it needs to know which
+    # way each axis runs; the edge costs below only ever use squared or
+    # already-absolute cell sizes, so the signed values are what the
+    # backends carry (see _finalize_direction).
+    signed_cellsize_x = cellsize_x * _coord_step_sign(raster, x)
+    signed_cellsize_y = cellsize_y * _coord_step_sign(raster, y)
+
     target_values = np.asarray(target_values, dtype=np.float64)
     max_distance_f = float(max_distance)
 
@@ -1371,28 +1420,28 @@ def _compute(raster, elevation, x, y, target_values, max_distance,
             raise NotImplementedError(
                 "geodesic mode is not yet supported for CuPy arrays")
         result_data = _surface_distance_cupy(
-            source_data, elev_data, cellsize_x, cellsize_y,
+            source_data, elev_data, signed_cellsize_x, signed_cellsize_y,
             max_distance_f, target_values, dy_arr, dx_arr, dd_arr, mode,
         )
     elif _is_dask_cupy_flag:
         if use_geodesic:
             raise NotImplementedError(
                 "geodesic mode is not yet supported for Dask+CuPy arrays")
         result_data = _surface_distance_dask_cupy(
-            source_data, elev_data, cellsize_x, cellsize_y,
+            source_data, elev_data, signed_cellsize_x, signed_cellsize_y,
             max_distance_f, target_values, dy_arr, dx_arr, dd_arr, mode,
         )
     elif isinstance(source_data, np.ndarray):
         if isinstance(elev_data, np.ndarray):
             result_data = _surface_distance_numpy(
-                source_data, elev_data, cellsize_x, cellsize_y,
+                source_data, elev_data, signed_cellsize_x, signed_cellsize_y,
                 max_distance_f, target_values, dy_arr, dx_arr, dd_arr,
                 dd_grid, use_geodesic, mode,
             )
         else:
             elev_np = np.asarray(elev_data)
             result_data = _surface_distance_numpy(
-                source_data, elev_np, cellsize_x, cellsize_y,
+                source_data, elev_np, signed_cellsize_x, signed_cellsize_y,
                 max_distance_f, target_values, dy_arr, dx_arr, dd_arr,
                 dd_grid, use_geodesic, mode,
             )
@@ -1401,7 +1450,7 @@ def _compute(raster, elevation, x, y, target_values, max_distance,
             raise NotImplementedError(
                 "geodesic mode is not yet supported for Dask arrays")
         result_data = _surface_distance_dask(
-            source_data, elev_data, cellsize_x, cellsize_y,
+            source_data, elev_data, signed_cellsize_x, signed_cellsize_y,
             max_distance_f, target_values, dy_arr, dx_arr, dd_arr, mode,
         )
     else:

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
@@ -299,6 +301,179 @@ def test_direction_cardinal_points():
     assert sd_dir[2, 1] == pytest.approx(360.0, abs=1.0)
 
 
+def test_direction_follows_y_axis_orientation():
+    """Bearings must use the y coordinate, not the row index (#3719).
+
+    A north-up raster has y descending with the row index. Reading the row
+    index as if it were y mirrors every bearing across the east-west axis.
+    """
+    source = np.zeros((3, 3), dtype=np.float64)
+    source[1, 1] = 1.0
+    elev = np.zeros((3, 3), dtype=np.float64)
+
+    def _build(y_coords):
+        return (
+            xr.DataArray(source, dims=['y', 'x'],
+                         coords={'y': y_coords,
+                                 'x': np.arange(3, dtype=np.float64)}),
+            xr.DataArray(elev, dims=['y', 'x'],
+                         coords={'y': y_coords,
+                                 'x': np.arange(3, dtype=np.float64)}),
+        )
+
+    # Descending y: row 0 is north of the source, so it points north (360).
+    north_up = _compute(surface_direction(
+        *_build(np.array([2.0, 1.0, 0.0]))))
+    assert north_up[0, 1] == pytest.approx(360.0, abs=1.0)
+    assert north_up[2, 1] == pytest.approx(180.0, abs=1.0)
+
+    # Ascending y: row 0 is south of the source, so it points south (180).
+    south_up = _compute(surface_direction(
+        *_build(np.array([0.0, 1.0, 2.0]))))
+    assert south_up[0, 1] == pytest.approx(180.0, abs=1.0)
+    assert south_up[2, 1] == pytest.approx(360.0, abs=1.0)
+
+    # East/west is unaffected by the y orientation.
+    for out in (north_up, south_up):
+        assert out[1, 0] == pytest.approx(90.0, abs=1.0)
+        assert out[1, 2] == pytest.approx(270.0, abs=1.0)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+@pytest.mark.parametrize("func", [surface_distance, surface_allocation])
+def test_distance_and_allocation_ignore_axis_direction(backend, func):
+    """Only the bearing cares which way the axes run (#3719).
+
+    _compute hands the backends signed cell sizes so _finalize_direction
+    can build a compass bearing. That is safe only because every
+    edge-cost consumer squares or abs()es the cell size first. Flipping
+    the y axis must leave distance and allocation untouched apart from
+    the row order.
+    """
+    source = np.zeros((6, 6), dtype=np.float64)
+    source[1, 1] = 1.0
+    source[4, 4] = 2.0
+    elev = np.random.default_rng(11).uniform(0, 20, (6, 6))
+
+    def _build(y_coords):
+        arrays = []
+        for data in (source, elev):
+            arr = xr.DataArray(
+                data.copy(), dims=['y', 'x'],
+                coords={'y': y_coords, 'x': np.arange(6, dtype=np.float64)},
+                attrs={'res': (1.0, 1.0)})
+            if backend == 'dask+numpy':
+                if da is None:
+                    pytest.skip("dask not installed")
+                arr.data = da.from_array(arr.data, chunks=(3, 3))
+            arrays.append(arr)
+        return arrays
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", UserWarning)
+        up = _compute(func(*_build(np.arange(6, dtype=np.float64))))
+        down = _compute(func(*_build(np.arange(6, dtype=np.float64)[::-1])))
+
+    np.testing.assert_allclose(down, up, rtol=1e-6, equal_nan=True)
+
+
+def test_direction_matches_proximity_direction():
+    """surface_direction on flat terrain agrees with proximity.direction.
+
+    Both document the same compass convention, so on zero relief with one
+    source they must return the same bearings (#3719).
+    """
+    from xrspatial.proximity import direction
+
+    source = np.zeros((5, 5), dtype=np.float64)
+    source[3, 1] = 1.0
+    elev = np.zeros((5, 5), dtype=np.float64)
+    y_coords = np.array([40.0, 39.0, 38.0, 37.0, 36.0])  # descending
+    x_coords = np.arange(5, dtype=np.float64)
+
+    raster = xr.DataArray(source, dims=['y', 'x'],
+                          coords={'y': y_coords, 'x': x_coords})
+    elevation = xr.DataArray(elev, dims=['y', 'x'],
+                             coords={'y': y_coords, 'x': x_coords})
+
+    np.testing.assert_allclose(
+        _compute(surface_direction(raster, elevation)),
+        _compute(direction(raster)),
+        rtol=1e-5,
+    )
+
+
+@pytest.mark.skipif(da is None, reason="dask not installed")
+@pytest.mark.parametrize("max_distance", [np.inf, 4.0])
+def test_dask_direction_matches_numpy(max_distance):
+    """Dask direction must match numpy on every chunk (#3719).
+
+    The unbounded path re-runs each tile with global source indices; those
+    have to be rebased before the bearing is measured against the block's
+    own index grid, or every chunk but the first is wrong.
+    """
+    source = np.zeros((8, 10), dtype=np.float64)
+    source[2, 3] = 1.0
+    elev = np.random.default_rng(7).uniform(0, 5, (8, 10))
+
+    raster_np = _make_raster(source, backend='numpy')
+    elev_np = _make_raster(elev, backend='numpy')
+    raster_dask = _make_raster(source, backend='dask+numpy', chunks=(4, 5))
+    elev_dask = _make_raster(elev, backend='dask+numpy', chunks=(4, 5))
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", UserWarning)
+        np_result = _compute(surface_direction(
+            raster_np, elev_np, max_distance=max_distance))
+        dask_result = _compute(surface_direction(
+            raster_dask, elev_dask, max_distance=max_distance))
+
+    np.testing.assert_allclose(dask_result, np_result, rtol=1e-5,
+                               equal_nan=True)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy'])
+@pytest.mark.parametrize("max_distance", [np.inf, 4.0])
+def test_direction_nan_mask_matches_distance(backend, max_distance):
+    """Direction is NaN exactly where distance is NaN (#3719).
+
+    Guards the tiled dask path, where a chunk whose source lives in a
+    neighbouring chunk must still get a bearing.
+    """
+    source = np.zeros((8, 10), dtype=np.float64)
+    source[2, 3] = 1.0
+    elev = np.random.default_rng(7).uniform(0, 5, (8, 10))
+
+    raster = _make_raster(source, backend=backend, chunks=(4, 5))
+    elevation = _make_raster(elev, backend=backend, chunks=(4, 5))
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", UserWarning)
+        dist = _compute(surface_distance(raster, elevation,
+                                         max_distance=max_distance))
+        bearing = _compute(surface_direction(raster, elevation,
+                                             max_distance=max_distance))
+
+    np.testing.assert_array_equal(np.isnan(bearing), np.isnan(dist))
+
+
+@pytest.mark.skipif(not has_cuda_and_cupy(), reason="cupy/cuda not available")
+def test_cupy_direction_matches_numpy():
+    """CuPy direction must match the numpy baseline (#3719)."""
+    source = np.zeros((8, 10), dtype=np.float64)
+    source[2, 3] = 1.0
+    elev = np.random.default_rng(7).uniform(0, 5, (8, 10))
+
+    np_result = _compute(surface_direction(
+        _make_raster(source), _make_raster(elev)))
+    cp_result = _compute(surface_direction(
+        _make_raster(source, backend='cupy'),
+        _make_raster(elev, backend='cupy')))
+
+    np.testing.assert_allclose(cp_result, np_result, rtol=1e-5,
+                               equal_nan=True)
+
+
 # ---------------------------------------------------------------------------
 # Tests — max_distance clipping
 # ---------------------------------------------------------------------------