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

[brendancol]

Closes #3719

surface_direction() was returning wrong bearings in two unrelated ways, and
fixing the second one uncovered a third. surface_distance() and
surface_allocation() were correct throughout.

What changed

• _finalize_direction() measured the offset to the source in row/column
  indices scaled by an absolute cell size, so the row index stood in for y.
  On a north-up raster (y descending as the row index rises) every bearing
  came out mirrored across the east-west axis, and surface_direction()
  disagreed with proximity.direction() on identical input. _compute() now
  reads the sign of each axis off the coordinate and passes the backends
  signed cell sizes. Every edge-cost site either squares the cell size or
  takes abs() of it first, so distances and allocations are byte-identical
  to before.

• _assemble_sd() handed _finalize_direction() the global source indices
  that _run_tile() records, while the bearing was measured against the
  block's own np.arange grid. The chunk offset stayed in the answer and
  every chunk except (0, 0) reported the bearing it would have had as a
  standalone raster. It now rebases to block-local indices. This is the
  default path: max_distance defaults to np.inf, which routes to
  _sd_dask_iterative().

• Rebasing broke the no-source mask, which tested src_row < 0 — a source
  in the chunk above rebases to a legitimately negative local index. A finite
  in-budget distance always comes with an assigned source, so the mask now
  keys on dist alone. Direction ends up with the same NaN mask as distance
  and allocation, which is what the docstring already promised.

Backends

numpy, cupy, dask+numpy, dask+cupy. All four checked against the numpy
baseline on a descending-y raster, bounded and unbounded, with zero NaN-mask
mismatches and zero numeric difference. The cupy eager path shares
_finalize_direction's index math and was affected by the orientation bug
too; dask+cupy inherits the tiled fix through the CPU iterative path.

Test plan

• pytest xrspatial/tests/test_surface_distance.py — 46 passed (37 before)
• pytest xrspatial/tests/test_proximity.py xrspatial/tests/test_cost_distance.py — 686 passed
• flake8 clean on both changed files (the F841 at
  test_surface_distance.py:234 predates this branch and is untouched)
• New: test_direction_follows_y_axis_orientation pins both axis
  directions; test_direction_matches_proximity_direction pins agreement
  with the sibling API; test_dask_direction_matches_numpy and
  test_direction_nan_mask_matches_distance cover bounded and unbounded
  dask; test_cupy_direction_matches_numpy closes the GPU gap
• Checked out main's surface_distance.py under the new tests:
  test_direction_follows_y_axis_orientation,
  test_direction_matches_proximity_direction, and both
  test_dask_direction_matches_numpy cases fail there and pass here.
  test_direction_nan_mask_matches_distance and
  test_cupy_direction_matches_numpy pass on main as well — they pin
  invariants this change had to preserve, and the mask one is what caught
  the third defect mid-fix.

Found by /sweep-accuracy against xrspatial/surface_distance.py. The state
CSV row for the module is in the second commit.

[brendancol]

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

Reviewed the full surface_distance.py and test_surface_distance.py at the
head of deep-sweep-accuracy-surface_distance-2026-08-16-01, not just the
diff. All three defects are real and the reproductions in #3719 replay against
main. Distances and allocations are untouched, which the existing 37 tests
confirm.

Blockers (must fix before merge)

None.

Suggestions (should fix, not blocking)

• xrspatial/surface_distance.py:642 — the cupy DIRECTION branch still
  masks on srow < 0. _finalize_direction dropped that clause at line
  399-406 because a rebased block-local index can legitimately be negative.
  The eager cupy path never rebases, so srow is always non-negative wherever
  dist is finite and the two masks agree today. But leaving one backend
  keying on the source index and the other on dist is the asymmetry that
  produced this bug in the first place. Make it cp.isinf(dist) | (dist > max_distance) so both paths state the same rule.

• xrspatial/surface_distance.py:1355-1359 — the signed cell size relies on
  an unwritten invariant. Passing signed values through the existing
  cellsize_x / cellsize_y parameters is correct only because every other
  consumer squares them (_compute_seeds_sd builds diag_dist from
  cellsize_x ** 2 + cellsize_y ** 2, _edge_seeds uses cardinal_dist ** 2) or takes abs() first (_surface_distance_dask,
  _surface_distance_dask_bounded, _surface_distance_dask_cupy). I walked
  every call site and the invariant holds, but nothing enforces it, and the
  next person to add a cell-size consumer has no way to know. A test asserting
  that a descending-y raster produces the same DISTANCE and ALLOCATION output
  as an ascending-y one would pin it cheaply.

Nits (optional improvements)

• xrspatial/surface_distance.py:344 — _coord_step_sign reads the span,
  not a step. values[-1] - values[0] gives the right answer for any
  monotonic axis and matches how calc_res already derives resolution, but
  surface_distance never checks monotonicity the way
  proximity._check_monotonic_coords does. On a scrambled coordinate the sign
  is whatever the endpoints happen to say. Not worth a validation pass, but
  worth a word in the docstring that the axis is assumed monotonic.

• Geodesic direction is still index-scaled. In method='geodesic' mode
  _finalize_direction scales index offsets by a cell size in degrees, so
  bearings drift away from the equator where a degree of longitude is shorter
  than a degree of latitude. This PR fixes the sign but not the scale. It
  predates the change and is not made worse by it; worth a separate issue
  rather than growing this diff.

What looks good

• The three defects are separated cleanly, and the commit message explains how
  fixing the second exposed the third rather than presenting it as one insight.
• Distances and allocations come out byte-identical, which is the right
  outcome for a direction-only fix and is what makes the signed-cell-size
  approach defensible.
• test_direction_matches_proximity_direction pins agreement with the sibling
  API instead of hard-coding numbers, so the two functions cannot drift apart
  again silently.
• test_direction_nan_mask_matches_distance covers bounded and unbounded dask
  and is what caught the third defect mid-fix. Keeping it in the suite is more
  valuable than the bearing assertions.
• The verification ran all four backends on a descending-y raster rather than
  reasoning about the GPU path from source.

Checklist

• Algorithm matches reference — bearings now agree with
  proximity.direction, the in-repo reference for this convention
• All implemented backends produce consistent results — numpy, cupy,
  dask+numpy, dask+cupy checked bounded and unbounded, zero difference
• NaN handling is correct — direction now shares the distance NaN mask
• Edge cases covered — both axis directions, no-source raster, chunks
  whose source lives in a neighbouring chunk
• Dask chunk boundaries handled correctly — the tiled path was the bug
• No premature materialization — srow - row_offset is a per-block numpy
  op inside an existing map_blocks chunk function
• Benchmark exists — benchmarks/benchmarks/surface_distance.py, and this
  change has no throughput effect
• README feature matrix — not applicable, no new function and no backend
  support change
• Docstrings present and accurate — the surface_direction compass
  convention in the public docstring is now what the function returns

[brendancol]

PR Review: follow-up pass after da98601

Re-read the changed hunks and re-ran the suite. Disposition of the first pass:

Fixed

• The cupy DIRECTION branch no longer masks on srow < 0. Both it and
  _finalize_direction now key on dist, so the two backends state the same
  rule and the GPU side cannot reintroduce the asymmetry.
• test_distance_and_allocation_ignore_axis_direction pins the invariant the
  signed cell size leans on: flip the y axis and distance and allocation come
  back unchanged. Parametrized over numpy and dask+numpy, and over
  surface_distance and surface_allocation. 50 tests pass (37 on main).
• _coord_step_sign's docstring now says the axis is assumed monotonic and
  points at calc_res making the same assumption.

Deferred, with a linked issue

• Geodesic bearings still scale index offsets by a cell size in degrees, so
  they tilt with cos(lat) away from the equator. Filed as #3733. It predates
  this PR, and the fix needs a decision this one should not make: whether to
  scale the longitude leg by cos(lat) at the pixel, or compute a forward
  azimuth from the two lat/lon pairs, which changes what "bearing" means over
  a long path since a great circle's initial azimuth is not its final one.

No new findings. Nothing outstanding from my side.

[brendancol]

Overlap notice: PR #3729 (test-coverage sweep) and PR #3731 (accuracy sweep) both fix the dask surface_direction bug — global source indices measured against chunk-local index grids — found independently in the same /deep-sweep run, along with #3713 filed by the documentation sweep. Both rewrite _finalize_direction and _tile_fn, so they will conflict.

They are not equivalent, and the difference matters:

• Fix surface_direction bearings on north-up rasters and on chunked dask input #3731 fixes three defects in surface_direction, not one. Besides the chunk-local index bug, it fixes bearings built from row/column indices scaled by an absolute cell size, which mirrors every bearing north/south on a standard north-up (y-descending) raster and disagrees with proximity.direction() on identical input — on all backends, eager included. It then fixes a no-source mask testing src_row < 0, which a rebased source from the chunk above can satisfy.
• Cover the surface_distance backends that had no tests, and fix dask surface_direction (#3725, #3713) #3729 fixes the chunk-local bug and brings the broad backend test coverage (37 to 106 tests, including the previously uncovered _surface_distance_dask_cupy).

The trap in merging #3729 alone: its direction tests assert that dask matches numpy. If numpy is itself mirrored, that parity holds while both backends are wrong, so the suite would go green over a live defect. #3731's test_direction_cardinal_points pins bearings against known compass directions instead, which is what actually catches it.

Suggested resolution: take #3731's source changes as the base (superset of the direction fixes), keep #3729's test additions, and reconcile the overlapping hunks in _finalize_direction / _tile_fn by hand. Then run both PRs' full test files against the merged tree — neither PR's own green run proves the combination.

Note for sequencing: PR #3724 adds a docstring caveat describing this bug as unfixed, so it should land after whichever fix merges, with that paragraph dropped.