diff --git a/xrspatial/contour.py b/xrspatial/contour.py index d1e8489a..a4c11c52 100644 --- a/xrspatial/contour.py +++ b/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)) diff --git a/xrspatial/tests/test_contour.py b/xrspatial/tests/test_contour.py index 18b06787..0066b66d 100644 --- a/xrspatial/tests/test_contour.py +++ b/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 # ---------------------------------------------------------------------------