Drop non-finite points and segments in kde/line_density (#3628)

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

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)
 # ---------------------------------------------------------------------------