diff --git a/xrspatial/kde.py b/xrspatial/kde.py index ff9e8a2bb..0d58d4511 100644 --- a/xrspatial/kde.py +++ b/xrspatial/kde.py @@ -519,6 +519,8 @@ def kde( x, y : array-like 1-D arrays of point coordinates. Alternatively, pass a GeoDataFrame of Point geometries as the first argument and leave *y* unset. + Points with non-finite (NaN or infinite) coordinates or weights + are dropped; a ``ValueError`` is raised if no finite points remain. weights : array-like, optional Per-point weights. Defaults to uniform weights of 1. column : str, optional @@ -589,6 +591,20 @@ def kde( else: w_arr = np.ones(n, dtype=np.float64) + # Drop points with non-finite coordinates or weights (#3628). Same + # policy as xrspatial.interpolate: without this, a single NaN poisons + # the auto-computed extent (all backends) or the whole grid (cupy + # gaussian, which has no cutoff), while the other backends silently + # skip the point. + valid = np.isfinite(x_arr) & np.isfinite(y_arr) & np.isfinite(w_arr) + if not valid.all(): + x_arr, y_arr, w_arr = x_arr[valid], y_arr[valid], w_arr[valid] + if x_arr.shape[0] == 0: + raise ValueError( + "kde(): no valid (finite) points remain after filtering " + "non-finite coordinates and weights" + ) + kid = _kernel_id(kernel) # -- Bandwidth --------------------------------------------------------- @@ -712,7 +728,9 @@ def line_density( Parameters ---------- x1, y1, x2, y2 : array-like - Start and end coordinates of each line segment. + Start and end coordinates of each line segment. Segments with + non-finite (NaN or infinite) endpoints or weights are dropped; + a ``ValueError`` is raised if no finite segments remain. weights : array-like, optional Per-segment weights. Defaults to uniform weights of 1. bandwidth : float or ``'silverman'`` @@ -754,6 +772,20 @@ def line_density( else: w_arr = np.ones(n, dtype=np.float64) + # Drop segments with non-finite endpoints or weights (#3628); a + # single NaN endpoint otherwise poisons the auto-computed extent + # and the output collapses to zeros with NaN coordinates. + valid = (np.isfinite(x1a) & np.isfinite(y1a) & + np.isfinite(x2a) & np.isfinite(y2a) & np.isfinite(w_arr)) + if not valid.all(): + x1a, y1a, x2a, y2a, w_arr = ( + x1a[valid], y1a[valid], x2a[valid], y2a[valid], w_arr[valid]) + if x1a.shape[0] == 0: + raise ValueError( + "line_density(): no valid (finite) segments remain after " + "filtering non-finite endpoints and weights" + ) + kid = _kernel_id(kernel) # Bandwidth from all endpoints diff --git a/xrspatial/tests/test_kde.py b/xrspatial/tests/test_kde.py index 9da1c0cd1..16830595a 100644 --- a/xrspatial/tests/test_kde.py +++ b/xrspatial/tests/test_kde.py @@ -198,6 +198,120 @@ def test_template_not_2d_raises(self): kde([0], [0], bandwidth=1.0, template=t) +# --------------------------------------------------------------------------- +# Non-finite inputs (#3628) +# --------------------------------------------------------------------------- + +class TestNonFiniteInputs: + """NaN/Inf points and weights are dropped, identically on all backends.""" + + def test_nan_point_dropped_matches_clean_input(self, simple_grid): + clean = kde([0.0, 1.0], [0.0, 1.0], bandwidth=1.0, + template=simple_grid) + with_nan = kde([0.0, 1.0, np.nan], [0.0, 1.0, 0.5], bandwidth=1.0, + template=simple_grid) + np.testing.assert_allclose(with_nan.values, clean.values, rtol=1e-12) + + def test_inf_point_dropped(self, simple_grid): + clean = kde([0.0, 1.0], [0.0, 1.0], bandwidth=1.0, + template=simple_grid) + with_inf = kde([0.0, 1.0, np.inf], [0.0, 1.0, 0.5], bandwidth=1.0, + template=simple_grid) + np.testing.assert_allclose(with_inf.values, clean.values, rtol=1e-12) + + def test_nan_weight_dropped(self, simple_grid): + clean = kde([0.0, 1.0], [0.0, 1.0], weights=[1.0, 2.0], + bandwidth=1.0, template=simple_grid) + with_nan = kde([0.0, 1.0, 0.5], [0.0, 1.0, 0.5], + weights=[1.0, 2.0, np.nan], + bandwidth=1.0, template=simple_grid) + np.testing.assert_allclose(with_nan.values, clean.values, rtol=1e-12) + assert not np.isnan(with_nan.values).any() + + def test_nan_point_auto_extent_not_poisoned(self): + """Without a template, the extent must come from finite points only.""" + result = kde([0.0, 1.0, np.nan], [0.0, 1.0, 0.5], bandwidth=1.0, + width=8, height=8) + assert np.isfinite(result.x.values).all() + assert np.isfinite(result.y.values).all() + # The filter runs before the extent is derived, so the result must + # be identical to calling kde on the finite points only. + clean = kde([0.0, 1.0], [0.0, 1.0], bandwidth=1.0, width=8, height=8) + np.testing.assert_array_equal(result.values, clean.values) + np.testing.assert_array_equal(result.x.values, clean.x.values) + np.testing.assert_array_equal(result.y.values, clean.y.values) + + def test_nan_point_silverman_bandwidth_not_poisoned(self): + """Silverman's rule must see only the finite points (#3628).""" + result = kde([0.0, 1.0, 2.0, np.nan], [0.0, 1.0, 2.0, 0.5], + bandwidth='silverman', width=8, height=8) + clean = kde([0.0, 1.0, 2.0], [0.0, 1.0, 2.0], + bandwidth='silverman', width=8, height=8) + np.testing.assert_array_equal(result.values, clean.values) + assert float(result.sum()) > 0.0 + + def test_geodataframe_nan_point_dropped(self, simple_grid): + gpd = pytest.importorskip("geopandas") + shapely_geometry = pytest.importorskip("shapely.geometry") + Point = shapely_geometry.Point + gdf = gpd.GeoDataFrame( + geometry=[Point(0.0, 0.0), Point(1.0, 1.0), Point(np.nan, 0.5)]) + clean = kde([0.0, 1.0], [0.0, 1.0], bandwidth=1.0, + template=simple_grid) + result = kde(gdf, bandwidth=1.0, template=simple_grid) + np.testing.assert_allclose(result.values, clean.values, rtol=1e-12) + + def test_all_nan_points_raises(self, simple_grid): + with pytest.raises(ValueError, match='finite'): + kde([np.nan, np.nan], [0.0, 1.0], bandwidth=1.0, + template=simple_grid) + + @cuda_and_cupy_available + def test_cupy_gaussian_nan_point_matches_numpy(self, simple_grid): + """Eager cupy gaussian previously returned an all-NaN grid (#3628).""" + import cupy + x = [0.0, 1.0, np.nan] + y = [0.0, 1.0, 0.5] + np_result = kde(x, y, bandwidth=1.0, kernel='gaussian', + template=simple_grid) + cupy_template = simple_grid.copy(data=cupy.asarray(simple_grid.values)) + cupy_result = kde(x, y, bandwidth=1.0, kernel='gaussian', + template=cupy_template) + result_np = cupy_result.data.get() + assert not np.isnan(result_np).any() + # rtol matches TestCuPyParity; atol covers fringe pixels past the + # CPU kernel's 4*bw cutoff, which the GPU kernel does not have. + np.testing.assert_allclose(result_np, np_result.values, + rtol=1e-2, atol=1e-6) + + def test_line_density_nan_endpoint_dropped(self): + clean = line_density([0.0], [0.0], [1.0], [1.0], bandwidth=0.5, + x_range=(-1, 2), y_range=(-1, 2), + width=16, height=16) + with_nan = line_density([0.0, np.nan], [0.0, 0.0], [1.0, 1.0], + [1.0, 1.0], bandwidth=0.5, + x_range=(-1, 2), y_range=(-1, 2), + width=16, height=16) + np.testing.assert_allclose(with_nan.values, clean.values, rtol=1e-12) + + def test_line_density_nan_endpoint_auto_extent(self): + result = line_density([0.0, np.nan], [0.0, 0.0], [1.0, 1.0], + [1.0, 1.0], bandwidth=0.5, width=16, height=16) + assert np.isfinite(result.x.values).all() + assert np.isfinite(result.y.values).all() + # Identical to the finite-only call, extent included. + clean = line_density([0.0], [0.0], [1.0], [1.0], + bandwidth=0.5, width=16, height=16) + np.testing.assert_array_equal(result.values, clean.values) + np.testing.assert_array_equal(result.x.values, clean.x.values) + np.testing.assert_array_equal(result.y.values, clean.y.values) + + def test_line_density_all_nan_raises(self): + with pytest.raises(ValueError, match='finite'): + line_density([np.nan], [0.0], [1.0], [1.0], bandwidth=0.5, + width=8, height=8) + + # --------------------------------------------------------------------------- # Memory guard (#1287) # ---------------------------------------------------------------------------