diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index cc43f3153..4db6dc5cf 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -24,6 +24,7 @@ mcda,2026-06-10,3146,MEDIUM,5,"Cat5 backend failures, all raise loudly (no wrong morphology,2026-04-30,"1397,1399",HIGH,2;5,HIGH fixed in #1397/PR #1398: morph_erode/dilate seeded centre cell into running min/max even when kernel[centre]==0 (all 4 backends). HIGH fixed in #1399/PR #1400: dask backends raised on 1xN/Nx1 kernels because empty-slice writeback (0:-0). multispectral,2026-03-30T14:00:00Z,1094,,, normalize,2026-05-01,,,,rescale and standardize across all 4 backends. NaN/inf filtered via isfinite mask before min/max/mean/std. Constant input handled (range=0 -> new_min; std=0 -> 0.0). Output dtype float64 consistently. Backend parity covered by test_matches_numpy. No accuracy issues found. +pathfinding,2026-07-03,3629;3630;3631,HIGH,1;3;5,"Cat1/3 HIGH+MEDIUM #3629: _get_pixel_id used int(abs(point-coord0)/cellsize) so out-of-bounds points on the coords[0] side folded to mirrored interior pixels (silent wrong path instead of ValueError) and truncation gave a half-cell floor bias plus fp flip at exact centers (0.3 on 0.1-res grid -> pixel 2); fixed with signed step + round-to-nearest-center. Cat5 HIGH+MEDIUM #3630: multi_stop_search crashed on dask+cupy (np.asarray(seg.values) hits cupy implicit-conversion TypeError) and returned numpy-backed output for cupy input while a_star_search preserves array type; note test_multi_stop_cupy_matches_numpy had a tautological conversion expression masking it (test-coverage sweep: no dask+cupy multi_stop test existed). Cat2/5 MEDIUM #3631: _hpa_star_search returned a partial finite cost trail when refinement failed (89 finite px on a wall-split 200x200 grid with unreachable goal) vs the all-NaN no-path contract elsewhere. Cat6 clean: a_star cost == scipy csgraph dijkstra on 8-connected 20x25 grid with anisotropic cells + random NaN barriers, friction and no-friction, delta 0.0; A->B==B->A symmetric; dask matches. Cat4: planar coordinate-unit distances (no haversine) is the library-wide convention, noted not flagged. CUDA available: cupy + dask+cupy paths executed (cupy is CPU-fallback by design, dask sparse A* verified vs numpy). HPA* suboptimality is inherent/documented, not flagged." perlin,2026-04-10T12:00:00Z,,,,Improved Perlin noise implementation correct. Fade/gradient functions verified. Backend-consistent. Continuous at cell boundaries. polygon_clip,2026-06-10,3186,HIGH,5,"Cat5 backend inconsistency: dask+cupy clip_polygon rasterizes the mask with a uniform chunk size from the raster's first chunk, then feeds raster+mask to da.map_blocks (positional block pairing). Non-uniform raster chunks gave the mask a different block layout -> IndexError/ValueError (or silent mis-stamp). Repro (8,6) rechunk ((3,5),(6,)) on dask+cupy raised ValueError Shapes do not align; dask+numpy was fine via xarray.where rechunk. Fix #3186/PR: rechunk cond to raster.data.chunks[-2:] before map_blocks; added non-uniform regression tests for dask+numpy and dask+cupy. use_cuda->gpu migration in that branch was already landed by #3089/#3122. CUDA available; cupy+dask+cupy verified, 25 tests pass. Cats 1-4 clean: numpy path uses raster.where, cupy path operates on raw arrays, NaN inputs preserved, no neighborhood ops/curvature. Prior fix #1197/#1200 (crop+all_touched) merged and unrelated." polygonize,2026-05-29,2606,HIGH,5,"Cat 5 HIGH: dask connectivity=8 cross-chunk merge filled diagonal notch where same-value regions meet only at a corner across a chunk boundary; total area exceeded raster. Hole ring was dropped because containment tested hole[0] (on exterior at pinch). Fixed via _ring_interior_point in PR for #2606. numpy, dask+numpy, dask+cupy area parity now holds; 4-conn was already correct. cupy + dask+cupy paths validated on GPU host. Other cats clean: NaN masked on numpy/cupy float paths (tested), _is_close handles +/-inf via exact-equality short-circuit, atol/rtol/simplify_tolerance reject NaN/inf, integer GPU CCL matches numpy." diff --git a/xrspatial/pathfinding.py b/xrspatial/pathfinding.py index 8f7a23bd0..9f65f371d 100644 --- a/xrspatial/pathfinding.py +++ b/xrspatial/pathfinding.py @@ -44,8 +44,22 @@ def _get_pixel_id(point, raster, xdim=None, ydim=None): x_coords = raster.coords[xdim].data cellsize_x, cellsize_y = get_dataarray_resolution(raster, xdim, ydim) - py = int(abs(point[0] - y_coords[0]) / cellsize_y) - px = int(abs(point[1] - x_coords[0]) / cellsize_x) + cellsize_x = abs(float(cellsize_x)) + cellsize_y = abs(float(cellsize_y)) + + # coords may be ascending or descending; use a signed step so points + # outside the raster get out-of-range indices (rejected by _is_inside) + # instead of folding back inside through abs(). + sign_y = -1.0 if len(y_coords) > 1 and y_coords[1] < y_coords[0] else 1.0 + sign_x = -1.0 if len(x_coords) > 1 and x_coords[1] < x_coords[0] else 1.0 + + # round to the nearest cell center (not truncate) so a point in the + # upper half of a cell, or float noise at an exact center, does not + # shift the pixel by one. + py = int(round( + (float(point[0]) - float(y_coords[0])) / (sign_y * cellsize_y))) + px = int(round( + (float(point[1]) - float(x_coords[0])) / (sign_x * cellsize_x))) # return index of row and column where the `point` located. return py, px @@ -833,8 +847,12 @@ def a_star_search(surface: xr.DataArray, 2D array of values to bin. start : array-like object of 2 numeric elements (y, x) or (lat, lon) coordinates of the starting point. + The point is mapped to the pixel whose cell center is nearest. + A point outside the raster bounds raises a ``ValueError``. goal : array like object of 2 numeric elements (y, x) or (lat, lon) coordinates of the goal location. + Mapped to the nearest cell center; a point outside the raster + bounds raises a ``ValueError``. barriers : array like object, default=[] List of values inside the surface which are barriers (cannot cross). diff --git a/xrspatial/tests/test_pathfinding.py b/xrspatial/tests/test_pathfinding.py index 4fccd29fc..6ffafb762 100644 --- a/xrspatial/tests/test_pathfinding.py +++ b/xrspatial/tests/test_pathfinding.py @@ -157,6 +157,83 @@ def test_a_star_search_connectivity( np.testing.assert_allclose(path_agg_4, result_4_connectivity, equal_nan=True) +# ----------------------------------------------------------------------- +# Point -> pixel mapping (issue #3629) +# ----------------------------------------------------------------------- + +def _make_unit_raster(h=5, w=5): + import xarray as xr + r = xr.DataArray(np.ones((h, w)), dims=['y', 'x'], + attrs={'res': (1.0, 1.0)}) + r['y'] = np.linspace(h - 1, 0, h) # descending + r['x'] = np.linspace(0, w - 1, w) + return r + + +def test_get_pixel_id_nearest_center(): + from xrspatial.pathfinding import _get_pixel_id + + agg = _make_unit_raster() + # exact centers map to themselves + assert _get_pixel_id((4.0, 0.0), agg, 'x', 'y') == (0, 0) + assert _get_pixel_id((0.0, 4.0), agg, 'x', 'y') == (4, 4) + # a point in the upper half of a cell maps to the nearest center, + # not the truncated lower-index cell + assert _get_pixel_id((4.0, 2.6), agg, 'x', 'y') == (0, 3) + assert _get_pixel_id((1.6, 0.0), agg, 'x', 'y') == (2, 0) + + # float noise at an exact center: on a 0.1-res grid the pixel-3 + # center is 0.30000000000000004; a user-typed 0.3 must still map + # to pixel 3 + import xarray as xr + agg2 = xr.DataArray(np.ones((2, 10)), dims=['y', 'x']) + agg2['y'] = [0.1, 0.0] + agg2['x'] = np.arange(10) * 0.1 + assert _get_pixel_id((0.0, 0.3), agg2, 'x', 'y') == (1, 3) + + +def test_get_pixel_id_ascending_coords(): + from xrspatial.pathfinding import _get_pixel_id + import xarray as xr + + agg = xr.DataArray(np.ones((5, 5)), dims=['y', 'x'], + attrs={'res': (1.0, 1.0)}) + agg['y'] = np.linspace(0, 4, 5) # ascending + agg['x'] = np.linspace(0, 4, 5) + assert _get_pixel_id((0.0, 0.0), agg, 'x', 'y') == (0, 0) + assert _get_pixel_id((4.0, 2.6), agg, 'x', 'y') == (4, 3) + # outside on the low side -> negative index (rejected by _is_inside) + py, px = _get_pixel_id((-1.0, 0.0), agg, 'x', 'y') + assert py < 0 + + +def test_a_star_search_out_of_bounds_raises(): + """Points outside the raster must raise on every side, not fold inside.""" + agg = _make_unit_raster() # y in [0, 4] descending, x in [0, 4] + + inside = (2.0, 2.0) + outside_points = [ + (5.0, 2.0), # above the top edge (y_coords[0] side) + (-1.0, 2.0), # below the bottom edge + (2.0, -1.0), # left of the first column (x_coords[0] side) + (2.0, 5.0), # right of the last column + ] + for pt in outside_points: + with pytest.raises(ValueError, match="start location outside"): + a_star_search(agg, pt, inside) + with pytest.raises(ValueError, match="goal location outside"): + a_star_search(agg, inside, pt) + + +def test_multi_stop_out_of_bounds_waypoint_raises(): + from xrspatial import multi_stop_search + + agg = _make_unit_raster() + # waypoint above the top edge previously folded to an interior row + with pytest.raises(ValueError, match="outside the surface bounds"): + multi_stop_search(agg, [(5.0, 0.0), (0.0, 4.0)]) + + # ----------------------------------------------------------------------- # Helper for multi-backend test rasters # -----------------------------------------------------------------------