diff --git a/xrspatial/tests/test_geodesic_slope.py b/xrspatial/tests/test_geodesic_slope.py index 64bb4240e..ad965917c 100644 --- a/xrspatial/tests/test_geodesic_slope.py +++ b/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.""" diff --git a/xrspatial/utils.py b/xrspatial/utils.py index 36f050f30..265a912a9 100644 --- a/xrspatial/utils.py +++ b/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]