Accept GeoDataFrames in interpolation + expose on .xrs accessor with coregister

README.md

@@ -445,9 +445,9 @@ For a broader catalog of spectral indices and sensor-specific band combinations,
 
 | Name | Description | Source | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
 |:----------:|:------------|:------:|:----------------------:|:--------------------:|:-------------------:|:------:|
-| [IDW](xrspatial/interpolate/_idw.py) | Inverse Distance Weighting from scattered points to a raster grid | Standard (IDW) | ✅ | 🔼 | 🔼 | 🔼 |
-| [Kriging](xrspatial/interpolate/_kriging.py) | Ordinary Kriging with automatic variogram fitting (spherical, exponential, gaussian) | Standard (ordinary kriging) | 🔼 | 🔼 | 🔼 | 🔼 |
-| [Spline](xrspatial/interpolate/_spline.py) | Thin Plate Spline interpolation with optional smoothing | Standard (TPS) | 🔼 | 🔼 | 🔼 | 🔼 |
+| [IDW](xrspatial/interpolate/_idw.py) | Inverse Distance Weighting from scattered points (arrays or a GeoDataFrame) to a raster grid | Standard (IDW) | ✅ | 🔼 | 🔼 | 🔼 |
+| [Kriging](xrspatial/interpolate/_kriging.py) | Ordinary Kriging with automatic variogram fitting (spherical, exponential, gaussian); accepts arrays or a GeoDataFrame | Standard (ordinary kriging) | 🔼 | 🔼 | 🔼 | 🔼 |
+| [Spline](xrspatial/interpolate/_spline.py) | Thin Plate Spline interpolation with optional smoothing; accepts arrays or a GeoDataFrame | Standard (TPS) | 🔼 | 🔼 | 🔼 | 🔼 |
 
 -----------
 
@@ -494,7 +494,7 @@ For a broader catalog of spectral indices and sensor-specific band combinations,
 
 | Name | Description | Source | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
 |:-----|:------------|:------:|:------------------:|:-----------------:|:---------------------:|:---------------------:|
-| [KDE](xrspatial/kde.py) | Point-to-raster kernel density estimation (Gaussian, Epanechnikov, quartic) | Silverman 1986 | 🔼 | 🔼 | 🔼 | 🔼 |
+| [KDE](xrspatial/kde.py) | Point-to-raster kernel density estimation (Gaussian, Epanechnikov, quartic); accepts arrays or a GeoDataFrame | Silverman 1986 | 🔼 | 🔼 | 🔼 | 🔼 |
 | [Line Density](xrspatial/kde.py) | Line-segment-to-raster density estimation | Standard | 🔼 |  |  |  |
 
 --------

docs/source/reference/interpolation.rst

@@ -4,6 +4,20 @@
 Interpolation
 *************
 
+.. note::
+
+   Each function accepts either raw ``x``, ``y``, ``z`` arrays or a
+   GeoDataFrame of Point geometries (``column`` selects the value to
+   interpolate, defaulting to the first numeric column).  The same
+   operations are available on the ``.xrs`` accessor, where the caller
+   raster supplies the output grid, chunks, and CRS::
+
+       dem.xrs.kriging(points_gdf, column="elevation")
+       dem.xrs.spline(points_gdf, coregister=True)
+
+   ``coregister=True`` reprojects a GeoDataFrame's points from their CRS
+   into the caller's CRS before interpolation.
+
 Inverse Distance Weighting (IDW)
 =================================
 .. autosummary::

docs/source/reference/kde.rst

@@ -9,6 +9,16 @@ KDE
    Kernel density estimation converts point or line data into
    continuous density surfaces on a raster grid.
 
+   ``kde`` accepts either raw ``x``, ``y`` arrays or a GeoDataFrame of
+   Point geometries (``column`` selects per-point weights).  It is also
+   available on the ``.xrs`` accessor, where the caller raster supplies the
+   output grid and CRS::
+
+       grid.xrs.kde(points_gdf, coregister=True)
+
+   ``coregister=True`` reprojects the points from their CRS into the
+   caller's CRS first.
+
 KDE
 ===
 .. autosummary::

xrspatial/accessor.py

@@ -761,6 +761,59 @@ def _wrap_accessor_attr(accessor, name):
     return _AccessorTool(attr)
 
 
+def _check_point_crs_match(src, target):
+    """Raise if point CRS and caller CRS disagree.
+
+    Both arguments are already-normalized ``pyproj.CRS`` instances (or
+    ``None``).  A ``None`` on either side is a no-op, mirroring how
+    ``rasterize`` treats a missing CRS.
+    """
+    if src is None or target is None:
+        return
+    if not src.equals(target):
+        def _label(crs):
+            epsg = crs.to_epsg()
+            return f"EPSG:{epsg}" if epsg is not None else crs.name
+        raise ValueError(
+            f"CRS mismatch: points are in {_label(src)} but the caller is "
+            f"in {_label(target)}. Pass coregister=True to reproject the "
+            f"points into the caller CRS."
+        )
+
+
+def _interp_on_accessor(obj, func, points, rest, *, coregister, **kwargs):
+    """Run interpolation *func* against caller raster *obj* as the template.
+
+    *obj* supplies the output grid, chunks, CRS, and attrs (passed as
+    ``template``).  *points* is the first argument to *func*: a GeoDataFrame,
+    or raw ``x`` for the legacy array form, in which case *rest* carries the
+    remaining positional coordinates (``y``, ``z``).  When *points* is a
+    GeoDataFrame and ``coregister=True``, its geometries are reprojected into
+    the caller's CRS before interpolation; otherwise a genuine CRS mismatch
+    raises (matching ``rasterize``'s default).  ``coregister`` needs a CRS on
+    both sides and only applies to GeoDataFrame input.
+    """
+    from xrspatial.interpolate._vector import is_geodataframe
+
+    if is_geodataframe(points):
+        target = _to_pyproj_crs(obj.attrs.get('crs'))
+        src = _to_pyproj_crs(points.crs)
+        if coregister:
+            if src is None or target is None:
+                raise ValueError(
+                    "coregister=True needs a CRS on both the points "
+                    "(GeoDataFrame.crs) and the caller (attrs['crs'])"
+                )
+            points = points.to_crs(target)
+        else:
+            _check_point_crs_match(src, target)
+    elif coregister:
+        raise ValueError(
+            "coregister=True requires a GeoDataFrame input carrying a CRS"
+        )
+    return func(points, *rest, template=obj, **kwargs)
+
+
 @xr.register_dataarray_accessor("xrs")
 class XrsSpatialDataArrayAccessor:
     """DataArray accessor exposing xarray-spatial operations."""
@@ -1189,17 +1242,54 @@ def mahalanobis(self, other_bands, **kwargs):
 
     # ---- Interpolation ----
 
-    def idw(self, x, y, z, **kwargs):
+    def kde(self, x, y=None, *, coregister=False, **kwargs):
+        """Kernel density of points onto this DataArray's grid.
+
+        Pass a GeoDataFrame of Point geometries as the first argument
+        (``column`` selects per-point weights), or raw ``x``/``y`` arrays.
+        With ``coregister=True`` a GeoDataFrame's points are reprojected
+        from their CRS into this array's CRS first.  See
+        :func:`xrspatial.kde`.
+        """
+        from .kde import kde
+        return _interp_on_accessor(self._obj, kde, x, (y,),
+                                   coregister=coregister, **kwargs)
+
+    def idw(self, x, y=None, z=None, *, coregister=False, **kwargs):
+        """Inverse-distance-weighted interpolation onto this grid.
+
+        Pass a GeoDataFrame (``column`` selects the value to interpolate)
+        or raw ``x``/``y``/``z`` arrays.  ``coregister=True`` reprojects a
+        GeoDataFrame's points into this array's CRS first.  See
+        :func:`xrspatial.idw`.
+        """
         from .interpolate import idw
-        return idw(x, y, z, self._obj, **kwargs)
+        return _interp_on_accessor(self._obj, idw, x, (y, z),
+                                   coregister=coregister, **kwargs)
+
+    def kriging(self, x, y=None, z=None, *, coregister=False, **kwargs):
+        """Ordinary-kriging interpolation onto this grid.
 
-    def kriging(self, x, y, z, **kwargs):
+        Pass a GeoDataFrame (``column`` selects the value to interpolate)
+        or raw ``x``/``y``/``z`` arrays.  ``coregister=True`` reprojects a
+        GeoDataFrame's points into this array's CRS first.  See
+        :func:`xrspatial.kriging`.
+        """
         from .interpolate import kriging
-        return kriging(x, y, z, self._obj, **kwargs)
+        return _interp_on_accessor(self._obj, kriging, x, (y, z),
+                                   coregister=coregister, **kwargs)
+
+    def spline(self, x, y=None, z=None, *, coregister=False, **kwargs):
+        """Thin-plate-spline interpolation onto this grid.
 
-    def spline(self, x, y, z, **kwargs):
+        Pass a GeoDataFrame (``column`` selects the value to interpolate)
+        or raw ``x``/``y``/``z`` arrays.  ``coregister=True`` reprojects a
+        GeoDataFrame's points into this array's CRS first.  See
+        :func:`xrspatial.spline`.
+        """
         from .interpolate import spline
-        return spline(x, y, z, self._obj, **kwargs)
+        return _interp_on_accessor(self._obj, spline, x, (y, z),
+                                   coregister=coregister, **kwargs)
 
     # ---- Preview ----
 
@@ -1858,6 +1948,56 @@ def rasterize(self, geometries, **kwargs):
             "to use as rasterize template"
         )
 
+    # ---- Interpolation ----
+
+    def _interp_template(self):
+        """First 2-D y/x variable to use as an interpolation template."""
+        ds = self._obj
+        for var in ds.data_vars:
+            da = ds[var]
+            if da.ndim == 2 and 'y' in da.dims and 'x' in da.dims:
+                return da
+        raise ValueError(
+            "Dataset has no 2D variable with 'y' and 'x' dimensions "
+            "to use as an interpolation template"
+        )
+
+    def kde(self, x, y=None, *, coregister=False, **kwargs):
+        """Kernel density of points onto a 2-D variable's grid.
+
+        See :meth:`XrsSpatialDataArrayAccessor.kde`.
+        """
+        from .kde import kde
+        return _interp_on_accessor(self._interp_template(), kde, x, (y,),
+                                   coregister=coregister, **kwargs)
+
+    def idw(self, x, y=None, z=None, *, coregister=False, **kwargs):
+        """Inverse-distance-weighted interpolation onto a 2-D variable's grid.
+
+        See :meth:`XrsSpatialDataArrayAccessor.idw`.
+        """
+        from .interpolate import idw
+        return _interp_on_accessor(self._interp_template(), idw, x, (y, z),
+                                   coregister=coregister, **kwargs)
+
+    def spline(self, x, y=None, z=None, *, coregister=False, **kwargs):
+        """Thin-plate-spline interpolation onto a 2-D variable's grid.
+
+        See :meth:`XrsSpatialDataArrayAccessor.spline`.
+        """
+        from .interpolate import spline
+        return _interp_on_accessor(self._interp_template(), spline, x, (y, z),
+                                   coregister=coregister, **kwargs)
+
+    def kriging(self, x, y=None, z=None, *, coregister=False, **kwargs):
+        """Ordinary-kriging interpolation onto a 2-D variable's grid.
+
+        See :meth:`XrsSpatialDataArrayAccessor.kriging`.
+        """
+        from .interpolate import kriging
+        return _interp_on_accessor(self._interp_template(), kriging, x, (y, z),
+                                   coregister=coregister, **kwargs)
+
     # ---- GeoTIFF I/O ----
 
     def to_geotiff(self, path, var=None, **kwargs):

xrspatial/interpolate/_idw.py

@@ -10,6 +10,7 @@
                              _validate_scalar, cuda_args, ngjit)
 
 from ._validation import extract_grid_coords, validate_points
+from ._vector import resolve_xyz
 
 try:
     import cupy
@@ -281,14 +282,16 @@ def _check_idw_memory(grid_pixels, k):
         )
 
 
-def idw(x, y, z, template, power=2.0, k=None,
-        fill_value=np.nan, name='idw'):
+def idw(x, y=None, z=None, template=None, power=2.0, k=None,
+        fill_value=np.nan, name='idw', *, column=None):
     """Inverse Distance Weighting interpolation.
 
     Parameters
     ----------
     x, y, z : array-like
-        Coordinates and values of scattered input points.
+        Coordinates and values of scattered input points.  Alternatively,
+        pass a GeoDataFrame of Point geometries as the first argument and
+        leave *y*/*z* unset; *template* and *column* must then be keywords.
     template : xr.DataArray
         2-D DataArray whose grid defines the output raster.
     power : float, default 2.0
@@ -301,6 +304,9 @@ def idw(x, y, z, template, power=2.0, k=None,
         Value for pixels with zero total weight.
     name : str, default 'idw'
         Name of the output DataArray.
+    column : str, optional
+        When the first argument is a GeoDataFrame, the column whose values
+        are interpolated.  Defaults to the first numeric column.
 
     Returns
     -------
@@ -313,6 +319,7 @@ def idw(x, y, z, template, power=2.0, k=None,
         ``(grid_pixels, k)`` distance and index arrays from the
         ``cKDTree`` query would exceed 80% of available memory.
     """
+    x, y, z = resolve_xyz(x, y, z, column=column, func_name='idw')
     _validate_raster(template, func_name='idw', name='template')
     x_arr, y_arr, z_arr = validate_points(x, y, z, func_name='idw')
     _validate_scalar(power, func_name='idw', name='power',

xrspatial/interpolate/_kriging.py

@@ -11,6 +11,7 @@
                              _validate_scalar)
 
 from ._validation import extract_grid_coords, validate_points
+from ._vector import resolve_xyz
 
 try:
     import cupy
@@ -479,14 +480,16 @@ def _check_kriging_memory(n_points, grid_pixels, is_dask=False):
         )
 
 
-def kriging(x, y, z, template, variogram_model='spherical', nlags=15,
-            return_variance=False, name='kriging'):
+def kriging(x, y=None, z=None, template=None, variogram_model='spherical',
+            nlags=15, return_variance=False, name='kriging', *, column=None):
     """Ordinary Kriging interpolation.
 
     Parameters
     ----------
     x, y, z : array-like
-        Coordinates and values of scattered input points.
+        Coordinates and values of scattered input points.  Alternatively,
+        pass a GeoDataFrame of Point geometries as the first argument and
+        leave *y*/*z* unset; *template* and *column* must then be keywords.
     template : xr.DataArray
         2-D DataArray whose grid defines the output raster.
     variogram_model : str, default 'spherical'
@@ -498,6 +501,9 @@ def kriging(x, y, z, template, variogram_model='spherical', nlags=15,
         If True, return ``(prediction, variance)`` tuple.
     name : str, default 'kriging'
         Name of the output DataArray.
+    column : str, optional
+        When the first argument is a GeoDataFrame, the column whose values
+        are interpolated.  Defaults to the first numeric column.
 
     Returns
     -------
@@ -512,6 +518,7 @@ def kriging(x, y, z, template, variogram_model='spherical', nlags=15,
         matrix, or prediction matrix) would exceed 80% of available
         memory.
     """
+    x, y, z = resolve_xyz(x, y, z, column=column, func_name='kriging')
     _validate_raster(template, func_name='kriging', name='template')
     x_arr, y_arr, z_arr = validate_points(x, y, z, func_name='kriging')
 

xrspatial/interpolate/_spline.py

@@ -12,6 +12,7 @@
                              _validate_scalar, cuda_args, ngjit)
 
 from ._validation import extract_grid_coords, validate_points
+from ._vector import resolve_xyz
 
 try:
     import cupy
@@ -300,20 +301,26 @@ def _check_spline_memory(n_points, grid_pixels):
         )
 
 
-def spline(x, y, z, template, smoothing=0.0, name='spline'):
+def spline(x, y=None, z=None, template=None, smoothing=0.0, name='spline',
+           *, column=None):
     """Thin Plate Spline interpolation.
 
     Parameters
     ----------
     x, y, z : array-like
-        Coordinates and values of scattered input points.
+        Coordinates and values of scattered input points.  Alternatively,
+        pass a GeoDataFrame of Point geometries as the first argument and
+        leave *y*/*z* unset; *template* and *column* must then be keywords.
     template : xr.DataArray
         2-D DataArray whose grid defines the output raster.
     smoothing : float, default 0.0
         Smoothing parameter.  0 forces exact interpolation through
         the data points.
     name : str, default 'spline'
         Name of the output DataArray.
+    column : str, optional
+        When the first argument is a GeoDataFrame, the column whose values
+        are interpolated.  Defaults to the first numeric column.
 
     Returns
     -------
@@ -325,6 +332,7 @@ def spline(x, y, z, template, smoothing=0.0, name='spline'):
         If the worst-case allocation (kernel block K or augmented
         system A) would exceed 80% of available memory.
     """
+    x, y, z = resolve_xyz(x, y, z, column=column, func_name='spline')
     _validate_raster(template, func_name='spline', name='template')
     x_arr, y_arr, z_arr = validate_points(x, y, z, func_name='spline')
     _validate_scalar(smoothing, func_name='spline', name='smoothing',