Fix contours() auto-levels silently dropping all contours with +/-inf

xrspatial/contour.py

@@ -8,6 +8,7 @@
 # The algorithm is embarrassingly parallel across quads and across contour
 # levels, making it well suited to Dask chunking and GPU execution.
 
+import warnings
 from typing import TYPE_CHECKING, List, Optional, Sequence, Tuple, Union
 
 import numpy as np
@@ -622,20 +623,34 @@ def contours(
 
     # Determine contour levels.
     if levels is None:
-        if da is not None and isinstance(agg.data, da.Array):
-            vmin, vmax = dask.compute(
-                da.nanmin(agg.data), da.nanmax(agg.data)
+        # Reduce over finite values only.  +/-inf cells would otherwise
+        # poison the range and make np.linspace emit non-finite levels,
+        # silently dropping every contour for the finite terrain.  An
+        # all-non-finite raster yields nan here, caught by the guard below;
+        # ignore the resulting empty-slice warning.
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                "ignore", "All-NaN slice encountered", RuntimeWarning
             )
-            vmin = float(vmin)
-            vmax = float(vmax)
-        elif cupy is not None and hasattr(agg.data, 'get'):
-            vmin = float(cupy.nanmin(agg.data))
-            vmax = float(cupy.nanmax(agg.data))
-        else:
-            vmin = float(np.nanmin(agg.values))
-            vmax = float(np.nanmax(agg.values))
+            if da is not None and isinstance(agg.data, da.Array):
+                finite = da.where(da.isfinite(agg.data), agg.data, np.nan)
+                vmin, vmax = dask.compute(
+                    da.nanmin(finite), da.nanmax(finite)
+                )
+                vmin = float(vmin)
+                vmax = float(vmax)
+            elif cupy is not None and hasattr(agg.data, 'get'):
+                finite = cupy.where(
+                    cupy.isfinite(agg.data), agg.data, cupy.nan
+                )
+                vmin = float(cupy.nanmin(finite))
+                vmax = float(cupy.nanmax(finite))
+            else:
+                finite = np.where(np.isfinite(agg.values), agg.values, np.nan)
+                vmin = float(np.nanmin(finite))
+                vmax = float(np.nanmax(finite))
 
-        if np.isnan(vmin) or np.isnan(vmax):
+        if not np.isfinite(vmin) or not np.isfinite(vmax):
             if return_type == "numpy":
                 return []
             return _to_geopandas([], crs=agg.attrs.get('crs', None))

xrspatial/tests/test_contour.py

@@ -132,6 +132,95 @@ def test_partial_nan(self):
             assert np.all(coords[:, 0] < nan_row_y + 1e-10)
 
 
+# ---------------------------------------------------------------------------
+# Non-finite (inf) handling in automatic level generation (issue #2797)
+# ---------------------------------------------------------------------------
+
+def _make_ramp_with_inf():
+    """Left-to-right ramp with one +inf and one -inf cell in the corner."""
+    data = _make_ramp(ny=5, nx=10)
+    data[0, 0] = np.inf
+    data[1, 0] = -np.inf
+    return data
+
+
+class TestAutoLevelsInf:
+
+    def test_auto_levels_ignore_inf(self):
+        """+/-inf must not poison auto-generated levels (#2797)."""
+        data = _make_ramp_with_inf()
+        agg = create_test_raster(data, backend='numpy')
+        result = contours(agg, n_levels=5)
+        # The finite ramp spans 0..9, so auto-levels still produce contours.
+        assert len(result) > 0
+        for level, _ in result:
+            assert np.isfinite(level)
+
+    def test_auto_levels_match_finite_range(self):
+        """Levels come from the finite min/max, not the inf extremes."""
+        data = _make_ramp_with_inf()
+        finite = _make_ramp(ny=5, nx=10)
+        agg = create_test_raster(data, backend='numpy')
+        finite_agg = create_test_raster(finite, backend='numpy')
+        inf_levels = sorted({lvl for lvl, _ in contours(agg, n_levels=5)})
+        finite_levels = sorted(
+            {lvl for lvl, _ in contours(finite_agg, n_levels=5)}
+        )
+        assert inf_levels == finite_levels
+
+    def test_all_inf_returns_empty(self):
+        """An entirely non-finite raster yields no contours, no crash."""
+        data = np.full((4, 4), np.inf, dtype=np.float64)
+        agg = create_test_raster(data, backend='numpy')
+        result = contours(agg)
+        assert result == []
+
+    def test_explicit_levels_unaffected_by_inf(self):
+        """Explicit levels bypass the range computation entirely."""
+        data = _make_ramp_with_inf()
+        agg = create_test_raster(data, backend='numpy')
+        result = contours(agg, levels=[4.5])
+        assert len(result) > 0
+
+    @dask_array_available
+    def test_auto_levels_ignore_inf_dask(self):
+        """Dask backend ignores inf in the lazy nanmin/nanmax path."""
+        data = _make_ramp_with_inf()
+        np_agg = create_test_raster(data, backend='numpy')
+        dask_agg = create_test_raster(
+            data, backend='dask+numpy', chunks=(3, 4)
+        )
+        np_levels = sorted({lvl for lvl, _ in contours(np_agg, n_levels=5)})
+        dk_levels = sorted({lvl for lvl, _ in contours(dask_agg, n_levels=5)})
+        assert len(dk_levels) > 0
+        assert np_levels == dk_levels
+
+    @cuda_and_cupy_available
+    def test_auto_levels_ignore_inf_cupy(self):
+        """CuPy backend ignores inf in the nanmin/nanmax path."""
+        data = _make_ramp_with_inf()
+        np_agg = create_test_raster(data, backend='numpy')
+        cupy_agg = create_test_raster(data, backend='cupy')
+        np_levels = sorted({lvl for lvl, _ in contours(np_agg, n_levels=5)})
+        cp_levels = sorted({lvl for lvl, _ in contours(cupy_agg, n_levels=5)})
+        assert len(cp_levels) > 0
+        assert np_levels == cp_levels
+
+    @dask_array_available
+    @cuda_and_cupy_available
+    def test_auto_levels_ignore_inf_dask_cupy(self):
+        """Dask+CuPy backend ignores inf in the lazy nanmin/nanmax path."""
+        data = _make_ramp_with_inf()
+        np_agg = create_test_raster(data, backend='numpy')
+        dc_agg = create_test_raster(
+            data, backend='dask+cupy', chunks=(3, 4)
+        )
+        np_levels = sorted({lvl for lvl, _ in contours(np_agg, n_levels=5)})
+        dc_levels = sorted({lvl for lvl, _ in contours(dc_agg, n_levels=5)})
+        assert len(dc_levels) > 0
+        assert np_levels == dc_levels
+
+
 # ---------------------------------------------------------------------------
 # Edge cases
 # ---------------------------------------------------------------------------