Fix geodesic coord resolution to prefer real lat/lon over numeric y/x (#2903)

* Fix geodesic coord resolution to prefer real lat/lon over numeric y/x (#2896)

_find_coord accepted the numeric dimension coord before scanning for a
coord named latitude/longitude, so a curvilinear raster with numeric y/x
index coords plus real 2-D lat/lon coords used the pixel indices as
lat/lon. The range check did not catch it because small indices fall
inside the geographic ranges.

Prefer an explicitly named lat/lon coord over the dimension coord. This
covers slope, aspect, and surface_distance, which share _extract_latlon_coords.

Add curvilinear coverage (dims y/x + 2-D lat/lon + numeric y/x) across the
four backends.

* Address review nit: document explicit-name coord resolution order (#2896)

Author: brendancol

xrspatial/tests/test_geodesic_slope.py

@@ -43,6 +43,40 @@ def _make_geo_raster(elev, lat_start, lat_end, lon_start, lon_end,
     return raster
 
 
+def _make_curvilinear_raster(elev, lat_start, lat_end, lon_start, lon_end,
+                             backend='numpy', chunks=(3, 3)):
+    """Build a curvilinear DataArray: dims ('y', 'x') with numeric y/x index
+    coords AND real 2-D lat/lon coords over the same geographic grid.
+
+    This is the layout that exposed the coordinate-resolution bug: the numeric
+    y/x pixel-index coords must not be used as lat/lon when real lat/lon coords
+    are present.
+    """
+    H, W = elev.shape
+    lat1d = np.linspace(lat_start, lat_end, H)
+    lon1d = np.linspace(lon_start, lon_end, W)
+    lon2d, lat2d = np.meshgrid(lon1d, lat1d)
+    raster = xr.DataArray(
+        elev.astype(np.float64),
+        dims=['y', 'x'],
+        coords={
+            'y': np.arange(H, dtype=np.float64),
+            'x': np.arange(W, dtype=np.float64),
+            'lat': (('y', 'x'), lat2d),
+            'lon': (('y', 'x'), lon2d),
+        },
+    )
+
+    if 'cupy' in backend:
+        import cupy
+        raster.data = cupy.asarray(raster.data)
+
+    if 'dask' in backend and da is not None:
+        raster.data = da.from_array(raster.data, chunks=chunks)
+
+    return raster
+
+
 def _flat_surface(H=6, W=8, elev=500.0):
     """Constant-elevation surface — slope should be 0 everywhere interior."""
     return np.full((H, W), elev, dtype=np.float64)
@@ -108,6 +142,87 @@ def test_north_tilted_has_positive_slope(self):
         assert np.all(interior > 0)
 
 
+class TestGeodesicSlopeCurvilinear:
+    """Curvilinear layout: dims ('y', 'x') with numeric y/x index coords plus
+    real 2-D lat/lon coords. The geodesic path must use the lat/lon coords, not
+    the pixel indices, so the result must match the equivalent 1-D lat/lon grid.
+
+    Pixel indices (0..N) fall inside the accepted geographic ranges, so the
+    range validation does not catch the mistake — only the slope value does.
+    """
+
+    def test_curvilinear_matches_1d_latlon(self):
+        elev = _east_tilted_surface(H=8, W=10, grade=100.0)
+        r_curv = _make_curvilinear_raster(elev, 40.0, 41.0, 10.0, 11.0)
+        r_ref = _make_geo_raster(elev, 40.0, 41.0, 10.0, 11.0)
+        s_curv = slope(r_curv, method='geodesic')
+        s_ref = slope(r_ref, method='geodesic')
+        np.testing.assert_allclose(
+            s_curv.values, s_ref.values, rtol=1e-5, equal_nan=True
+        )
+
+    def test_curvilinear_ignores_pixel_index_coords(self):
+        """Slope must reflect the real geographic grid, not the 0..N indices.
+
+        Using the pixel indices as lat/lon collapses the east tilt to a tiny
+        value (~0.007 vs ~0.067), so a correct interior slope is the signal.
+        """
+        elev = _east_tilted_surface(H=8, W=10, grade=100.0)
+        r_curv = _make_curvilinear_raster(elev, 40.0, 41.0, 10.0, 11.0)
+        interior = slope(r_curv, method='geodesic').values[1:-1, 1:-1]
+        assert np.all(np.isfinite(interior))
+        assert np.all(interior > 0.05)
+
+
+@dask_array_available
+class TestGeodesicSlopeCurvilinearDask:
+
+    def test_curvilinear_numpy_equals_dask(self):
+        elev = _east_tilted_surface(H=8, W=10, grade=100.0)
+        r_np = _make_curvilinear_raster(elev, 40.0, 41.0, 10.0, 11.0,
+                                        backend='numpy')
+        r_da = _make_curvilinear_raster(elev, 40.0, 41.0, 10.0, 11.0,
+                                        backend='dask+numpy', chunks=(4, 5))
+        s_np = slope(r_np, method='geodesic')
+        s_da = slope(r_da, method='geodesic')
+        np.testing.assert_allclose(
+            s_np.values, s_da.values, rtol=1e-5, equal_nan=True
+        )
+
+
+@cuda_and_cupy_available
+class TestGeodesicSlopeCurvilinearCupy:
+
+    def test_curvilinear_numpy_equals_cupy(self):
+        elev = _east_tilted_surface(H=8, W=10, grade=100.0)
+        r_np = _make_curvilinear_raster(elev, 40.0, 41.0, 10.0, 11.0,
+                                        backend='numpy')
+        r_cu = _make_curvilinear_raster(elev, 40.0, 41.0, 10.0, 11.0,
+                                        backend='cupy')
+        s_np = slope(r_np, method='geodesic')
+        s_cu = slope(r_cu, method='geodesic')
+        np.testing.assert_allclose(
+            s_np.values, s_cu.data.get(), rtol=1e-5, equal_nan=True
+        )
+
+
+@dask_array_available
+@cuda_and_cupy_available
+class TestGeodesicSlopeCurvilinearDaskCupy:
+
+    def test_curvilinear_numpy_equals_dask_cupy(self):
+        elev = _east_tilted_surface(H=8, W=10, grade=100.0)
+        r_np = _make_curvilinear_raster(elev, 40.0, 41.0, 10.0, 11.0,
+                                        backend='numpy')
+        r_dc = _make_curvilinear_raster(elev, 40.0, 41.0, 10.0, 11.0,
+                                        backend='dask+cupy', chunks=(4, 5))
+        s_np = slope(r_np, method='geodesic')
+        s_dc = slope(r_dc, method='geodesic')
+        np.testing.assert_allclose(
+            s_np.values, s_dc.data.compute().get(), rtol=1e-5, equal_nan=True
+        )
+
+
 class TestGeodesicSlopeLatitudeInvariance:
     """Same physical slope at equator vs 60N should give similar geodesic slope."""
 

xrspatial/utils.py

@@ -935,6 +935,11 @@ def warn_if_unit_mismatch(agg: xr.DataArray) -> None:
 # Known dimension / coordinate names (lower-cased for matching)
 _LAT_NAMES = {'lat', 'latitude', 'y'}
 _LON_NAMES = {'lon', 'longitude', 'x'}
+# Names that unambiguously mean geographic lat/lon. These take precedence
+# over a numeric dimension coord so a curvilinear raster with numeric y/x
+# index coords plus real lat/lon coords resolves to the lat/lon coords.
+_EXPLICIT_LAT_NAMES = {'lat', 'latitude'}
+_EXPLICIT_LON_NAMES = {'lon', 'longitude'}
 
 
 def _extract_latlon_coords(agg: xr.DataArray):
@@ -965,9 +970,11 @@ def _extract_latlon_coords(agg: xr.DataArray):
     dim_y, dim_x = agg.dims[-2], agg.dims[-1]
 
     # --- locate lat coordinate ---
-    lat_coord = _find_coord(agg, dim_y, _LAT_NAMES, 'latitude')
+    lat_coord = _find_coord(agg, dim_y, _LAT_NAMES, _EXPLICIT_LAT_NAMES,
+                            'latitude')
     # --- locate lon coordinate ---
-    lon_coord = _find_coord(agg, dim_x, _LON_NAMES, 'longitude')
+    lon_coord = _find_coord(agg, dim_x, _LON_NAMES, _EXPLICIT_LON_NAMES,
+                            'longitude')
 
     lat_vals = np.asarray(lat_coord.values, dtype=np.float64)
     lon_vals = np.asarray(lon_coord.values, dtype=np.float64)
@@ -994,15 +1001,30 @@ def _extract_latlon_coords(agg: xr.DataArray):
     return lat_2d, lon_2d
 
 
-def _find_coord(agg, dim_name, known_names, label):
-    """Find a coordinate matching *dim_name* or one of *known_names*."""
-    # 1) Try the dimension name directly
+def _find_coord(agg, dim_name, known_names, explicit_names, label):
+    """Find a coordinate matching *dim_name* or one of *known_names*.
+
+    A coordinate whose name is unambiguously geographic (*explicit_names*,
+    e.g. ``lat``/``longitude``) is preferred over the dimension coord. This
+    keeps a curvilinear raster with numeric ``y``/``x`` index coords plus
+    real lat/lon coords from silently using the pixel indices as lat/lon.
+    If several explicit names are present (e.g. both ``lat`` and
+    ``latitude``), the first one in coord order wins.
+    """
+    # 1) Prefer an explicitly named geographic coordinate (lat/lon).
+    for name in agg.coords:
+        if str(name).lower() in explicit_names:
+            coord = agg.coords[name]
+            if np.issubdtype(coord.dtype, np.number):
+                return coord
+
+    # 2) Fall back to the dimension name directly.
     if dim_name in agg.coords:
         coord = agg.coords[dim_name]
         if np.issubdtype(coord.dtype, np.number):
             return coord
 
-    # 2) Scan all coords for a known name
+    # 3) Scan all coords for any other known name (e.g. y/x).
     for name in agg.coords:
         if str(name).lower() in known_names:
             coord = agg.coords[name]