Validate planar cell size in slope() (#2897)

xrspatial/slope.py

@@ -360,7 +360,9 @@ def slope(agg: xr.DataArray,
     size. If the coordinates are in degrees (lat/lon) but the elevation values
     are in meters, the result is wrong by orders of magnitude. When this
     mismatch is detected, a ``UserWarning`` is emitted suggesting you reproject
-    to a projected CRS or use ``method='geodesic'``.
+    to a projected CRS or use ``method='geodesic'``. The ``'planar'`` method also
+    raises a ``ValueError`` if the cell size on either axis is zero, negative, or
+    non-finite, since those values cannot produce a valid slope.
 
     References
     ----------
@@ -401,6 +403,17 @@ def slope(agg: xr.DataArray,
     if method == 'planar':
         warn_if_unit_mismatch(agg)
         cellsize_x, cellsize_y = get_dataarray_resolution(agg)
+        # Reject negatives too (curvature() only checks == 0): a negative cell
+        # size would silently flip the slope sign rather than error out.
+        if (not np.isfinite(cellsize_x) or not np.isfinite(cellsize_y)
+                or cellsize_x <= 0 or cellsize_y <= 0):
+            raise ValueError(
+                "slope(method='planar') requires a positive, finite cell size "
+                f"on both axes; got cellsize_x={cellsize_x!r}, "
+                f"cellsize_y={cellsize_y!r}. "
+                "Set agg.attrs['res'] to a (x, y) tuple of positive floats, "
+                "or attach numeric x/y coordinates to the DataArray."
+            )
         mapper = ArrayTypeFunctionMapping(
             numpy_func=_run_numpy,
             cupy_func=_run_cupy,

xrspatial/tests/test_slope.py

@@ -473,6 +473,94 @@ def test_center_nan_propagates_dask_cupy():
                                   verify_attrs=False)
 
 
+# ---------------------------------------------------------------------------
+# Planar cell-size validation (issue #2897)
+# ---------------------------------------------------------------------------
+
+INVALID_RES = [
+    (0, 0),         # zero on both axes
+    (0, 1),         # zero on x
+    (1, 0),         # zero on y
+    (-1, 1),        # negative on x
+    (1, -1),        # negative on y
+    (np.inf, 1),    # inf on x
+    (1, np.inf),    # inf on y
+    (np.nan, 1),    # nan on x
+    (1, np.nan),    # nan on y
+]
+
+
+@pytest.mark.parametrize("res", INVALID_RES)
+def test_planar_invalid_cellsize_numpy(res):
+    """slope(method='planar') must reject zero, negative, or non-finite cell
+    sizes with a clear ValueError instead of a raw ZeroDivisionError, silent
+    NaN interior, or plausible-but-wrong values.
+    """
+    data = np.zeros((5, 5), dtype=np.float32)
+    agg = create_test_raster(data, backend='numpy', attrs={'res': res})
+    with pytest.raises(ValueError, match="positive, finite cell size"):
+        slope(agg)
+
+
+@dask_array_available
+@pytest.mark.parametrize("res", INVALID_RES)
+def test_planar_invalid_cellsize_dask_numpy(res):
+    data = np.zeros((5, 5), dtype=np.float32)
+    agg = create_test_raster(data, backend='dask+numpy',
+                             attrs={'res': res}, chunks=(3, 3))
+    with pytest.raises(ValueError, match="positive, finite cell size"):
+        slope(agg)
+
+
+@cuda_and_cupy_available
+@pytest.mark.parametrize("res", INVALID_RES)
+def test_planar_invalid_cellsize_cupy(res):
+    data = np.zeros((5, 5), dtype=np.float32)
+    agg = create_test_raster(data, backend='cupy', attrs={'res': res})
+    with pytest.raises(ValueError, match="positive, finite cell size"):
+        slope(agg)
+
+
+@dask_array_available
+@cuda_and_cupy_available
+@pytest.mark.parametrize("res", INVALID_RES)
+def test_planar_invalid_cellsize_dask_cupy(res):
+    data = np.zeros((5, 5), dtype=np.float32)
+    agg = create_test_raster(data, backend='dask+cupy',
+                             attrs={'res': res}, chunks=(3, 3))
+    with pytest.raises(ValueError, match="positive, finite cell size"):
+        slope(agg)
+
+
+def test_planar_valid_cellsize_still_works():
+    """A normal positive cell size must not trip the new guard."""
+    data = np.zeros((5, 5), dtype=np.float32)
+    agg = create_test_raster(data, backend='numpy', attrs={'res': (10, 10)})
+    result = slope(agg)
+    # slope() sets attrs['units'] = 'degrees' (#2894), so the output attrs
+    # legitimately differ from the input; the units contract is covered
+    # separately by the units tests below.
+    general_output_checks(agg, result, verify_attrs=False)
+
+
+def test_geodesic_unaffected_by_invalid_res_attr():
+    """The geodesic path takes its scale from lat/lon coords, so a bogus 'res'
+    attribute must not trigger the planar cell-size guard.
+    """
+    lat = np.linspace(40.0, 40.4, 5)
+    lon = np.linspace(-105.0, -104.6, 5)
+    data = np.zeros((5, 5), dtype=np.float64)
+    agg = xr.DataArray(
+        data,
+        coords={'y': lat, 'x': lon},
+        dims=['y', 'x'],
+        attrs={'res': (0, 0)},
+    )
+    # Must not raise the planar guard; geodesic uses lat/lon coords.
+    result = slope(agg, method='geodesic')
+    general_output_checks(agg, result, verify_attrs=False)
+
+
 # slope is reported in degrees, so the output must advertise units='degrees'
 # rather than leaking the input elevation units (e.g. 'm'). Regression for #2894.
 def test_units_overridden_to_degrees_numpy():