nd interpolants added

pydata · Oct 22, 2024 · 60f5afb · 60f5afb
1 parent 5632c8e
commit 60f5afb
Show file tree

Hide file tree

Showing 6 changed files with 367 additions and 148 deletions.
diff --git a/doc/user-guide/interpolation.rst b/doc/user-guide/interpolation.rst
@@ -132,10 +132,11 @@ It is now possible to safely compute the difference ``other - interpolated``.
 Interpolation methods
 ---------------------
 
-We use :py:class:`scipy.interpolate.interp1d` for 1-dimensional interpolation.
+We use either :py:class:`scipy.interpolate.interp1d` or special interpolants from
+:py:class:`scipy.interpolate` for 1-dimensional interpolation (see :py:meth:`~xarray.Dataset.interp`).
 For multi-dimensional interpolation, an attempt is first made to decompose the
 interpolation in a series of 1-dimensional interpolations, in which case
-:py:class:`scipy.interpolate.interp1d` is used. If a decomposition cannot be
+the relevant 1-dimensional interpolator is used. If a decomposition cannot be
 made (e.g. with advanced interpolation), :py:func:`scipy.interpolate.interpn` is
 used.
 

diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
@@ -2217,19 +2217,44 @@ def interp(
         coords: Mapping[Any, Any] | None = None,
         method: InterpOptions = "linear",
         assume_sorted: bool = False,
+        reduce: bool = True,
         kwargs: Mapping[str, Any] | None = None,
         **coords_kwargs: Any,
     ) -> Self:
-        """Interpolate a DataArray onto new coordinates
+        """
+        Interpolate a DataArray onto new coordinates.
+
+        Performs univariate or multivariate interpolation of a Dataset onto new coordinates,
+        utilizing either NumPy or SciPy interpolation routines.
+
+        When interpolating along multiple dimensions, the process attempts to decompose the
+        interpolation into independent interpolations along one dimension at a time, unless
+        `reduce=False` is passed. The specific interpolation method and dimensionality
+        determine which interpolant is used:
+
+        1. **Interpolation along one dimension of 1D data (`method='linear'`)**
+            - Uses :py:class:`numpy.interp`, unless `fill_value='extrapolate'` is provided via `kwargs`.
+
+        2. **Interpolation along one dimension of N-dimensional data (N ≥ 1)**
+            - Methods {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "quintic", "polynomial"}
+                use :py:class:`scipy.interpolate.interp1d`, unless conditions permit the use of :py:class:`numpy.interp`
+                (as in the case of `method='linear'` for 1D data).
+            - If `method='polynomial'`, the `order` keyword argument must also be provided. In this case,
+                :py:class:`scipy.interpolate.interp1d` is called with `kind=order`.
+
+        3. **Special interpolants for interpolation along one dimension of N-dimensional data (N ≥ 1)**
+            - Depending on the `method`, the following interpolants from :py:class:`scipy.interpolate` are used:
+                - `"pchip"`: :py:class:`scipy.interpolate.PchipInterpolator`
+                - `"barycentric"`: :py:class:`scipy.interpolate.BarycentricInterpolator`
+                - `"krogh"`: :py:class:`scipy.interpolate.KroghInterpolator`
+                - `"akima"` or `"makima"`: :py:class:`scipy.interpolate.Akima1dInterpolator`
+                    (`makima` is handled by passing `makima` to `method`).
 
-        Performs univariate or multivariate interpolation of a DataArray onto
-        new coordinates using scipy's interpolation routines. If interpolating
-        along an existing dimension,  either :py:class:`scipy.interpolate.interp1d`
-        or a 1-dimensional scipy interpolator (e.g. :py:class:`scipy.interpolate.KroghInterpolator`)
-        is called. When interpolating along multiple existing dimensions, an
-        attempt is made to decompose the interpolation into multiple
-        1-dimensional interpolations. If this is possible, the 1-dimensional interpolator is called.
-        Otherwise, :py:func:`scipy.interpolate.interpn` is called.
+        4. **Interpolation along multiple dimensions of multi-dimensional data**
+            - Uses :py:func:`scipy.interpolate.interpn` for methods {"linear", "nearest", "slinear",
+                "cubic", "quintic", "pchip"}.
+
+        Out-of-range values are filled with NaN, unless specified otherwise via `kwargs` to the numpy/scipy interpolant.
 
         Parameters
         ----------
@@ -2238,20 +2263,16 @@ def interp(
             New coordinate can be a scalar, array-like or DataArray.
             If DataArrays are passed as new coordinates, their dimensions are
             used for the broadcasting. Missing values are skipped.
-        method : {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "polynomial"}, default: "linear"
-            The method used to interpolate. The method should be supported by
-            the scipy interpolator:
-
-            - ``interp1d``: {"linear", "nearest", "zero", "slinear",
-              "quadratic", "cubic", "polynomial"}
-            - ``interpn``: {"linear", "nearest"}
-
-            If ``"polynomial"`` is passed, the ``order`` keyword argument must
-            also be provided.
+        method : str
+            Interpolation method to use (see descriptions above).
         assume_sorted : bool, default: False
             If False, values of x can be in any order and they are sorted
             first. If True, x has to be an array of monotonically increasing
             values.
+        reduce : bool, default: True
+            If True, the interpolation is decomposed into independent interpolations along one dimension at a time,
+            where the interpolation coordinates are independent. Setting this to be True alters the behavior of certain
+            multi-dimensional interpolants compared to the default SciPy output.
         kwargs : dict-like or None, default: None
             Additional keyword arguments passed to scipy's interpolator. Valid
             options and their behavior depend whether ``interp1d`` or
@@ -2267,12 +2288,14 @@ def interp(
 
         Notes
         -----
-        scipy is required.
+        - SciPy is required for certain interpolation methods.
+        - Allowing `reduce=True` (the default) may alter the behavior of interpolation along multiple dimensions
+            compared to the default behavior in SciPy.
 
         See Also
         --------
-        scipy.interpolate.interp1d
-        scipy.interpolate.interpn
+        Dataset.interp
+        Dataset.reindex_like
 
         :doc:`xarray-tutorial:fundamentals/02.2_manipulating_dimensions`
             Tutorial material on manipulating data resolution using :py:func:`~xarray.DataArray.interp`
@@ -2354,6 +2377,7 @@ def interp(
             method=method,
             kwargs=kwargs,
             assume_sorted=assume_sorted,
+            reduce=reduce,
             **coords_kwargs,
         )
         return self._from_temp_dataset(ds)
@@ -2363,48 +2387,75 @@ def interp_like(
         other: T_Xarray,
         method: InterpOptions = "linear",
         assume_sorted: bool = False,
+        reduce: bool = True,
         kwargs: Mapping[str, Any] | None = None,
     ) -> Self:
         """Interpolate this object onto the coordinates of another object,
         filling out of range values with NaN.
 
-        If interpolating along a single existing dimension,
-        :py:class:`scipy.interpolate.interp1d` is called. When interpolating
-        along multiple existing dimensions, an attempt is made to decompose the
-        interpolation into multiple 1-dimensional interpolations. If this is
-        possible, :py:class:`scipy.interpolate.interp1d` is called. Otherwise,
-        :py:func:`scipy.interpolate.interpn` is called.
+        When interpolating along multiple dimensions, the process attempts to decompose the
+        interpolation into independent interpolations along one dimension at a time, unless
+        `reduce=False` is passed. The specific interpolation method and dimensionality
+        determine which interpolant is used:
+
+        1. **Interpolation along one dimension of 1D data (`method='linear'`)**
+            - Uses :py:class:`numpy.interp`, unless `fill_value='extrapolate'` is provided via `kwargs`.
+
+        2. **Interpolation along one dimension of N-dimensional data (N ≥ 1)**
+            - Methods {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "quintic", "polynomial"}
+                use :py:class:`scipy.interpolate.interp1d`, unless conditions permit the use of :py:class:`numpy.interp`
+                (as in the case of `method='linear'` for 1D data).
+            - If `method='polynomial'`, the `order` keyword argument must also be provided.
+
+        3. **Special interpolants for interpolation along one dimension of N-dimensional data (N ≥ 1)**
+            - Depending on the `method`, the following interpolants from :py:class:`scipy.interpolate` are used:
+                - `"pchip"`: :py:class:`scipy.interpolate.PchipInterpolator`
+                - `"barycentric"`: :py:class:`scipy.interpolate.BarycentricInterpolator`
+                - `"krogh"`: :py:class:`scipy.interpolate.KroghInterpolator`
+                - `"akima"` or `"makima"`: :py:class:`scipy.interpolate.Akima1dInterpolator`
+                    (`makima` is handled by passing the `makima` flag).
+
+        4. **Interpolation along multiple dimensions of multi-dimensional data**
+            - Uses :py:func:`scipy.interpolate.interpn` for methods {"linear", "nearest", "slinear",
+                "cubic", "quintic", "pchip"}.
 
         Parameters
         ----------
         other : Dataset or DataArray
             Object with an 'indexes' attribute giving a mapping from dimension
             names to an 1d array-like, which provides coordinates upon
             which to index the variables in this dataset. Missing values are skipped.
-        method : {"linear", "nearest", "zero", "slinear", "quadratic", "cubic", "polynomial"}, default: "linear"
-            The method used to interpolate. The method should be supported by
-            the scipy interpolator:
-
-            - {"linear", "nearest", "zero", "slinear", "quadratic", "cubic",
-              "polynomial"} when ``interp1d`` is called.
-            - {"linear", "nearest"} when ``interpn`` is called.
-
-            If ``"polynomial"`` is passed, the ``order`` keyword argument must
-            also be provided.
+        method : str
+            Interpolation method to use (see descriptions above).
         assume_sorted : bool, default: False
             If False, values of coordinates that are interpolated over can be
             in any order and they are sorted first. If True, interpolated
             coordinates are assumed to be an array of monotonically increasing
             values.
+        reduce : bool, default: True
+            If True, the interpolation is decomposed into independent interpolations along one dimension at a time,
+            where the interpolation coordinates are independent. Setting this to be True alters the behavior of certain
+            multi-dimensional interpolants compared to the default SciPy output.
         kwargs : dict, optional
-            Additional keyword passed to scipy's interpolator.
+            Additional keyword arguments passed to the interpolant.
 
         Returns
         -------
         interpolated : DataArray
             Another dataarray by interpolating this dataarray's data along the
             coordinates of the other object.
 
+        Notes
+        -----
+        scipy is required.
+        If the dataarray has object-type coordinates, reindex is used for these
+        coordinates instead of the interpolation.
+
+        See Also
+        --------
+        DataArray.interp
+        DataArray.reindex_like
+
         Examples
         --------
         >>> data = np.arange(12).reshape(4, 3)
@@ -2460,24 +2511,18 @@ def interp_like(
         Coordinates:
           * x        (x) int64 32B 10 20 30 40
           * y        (y) int64 24B 70 80 90
-
-        Notes
-        -----
-        scipy is required.
-        If the dataarray has object-type coordinates, reindex is used for these
-        coordinates instead of the interpolation.
-
-        See Also
-        --------
-        DataArray.interp
-        DataArray.reindex_like
         """
+
         if self.dtype.kind not in "uifc":
             raise TypeError(
                 f"interp only works for a numeric type array. Given {self.dtype}."
             )
         ds = self._to_temp_dataset().interp_like(
-            other, method=method, kwargs=kwargs, assume_sorted=assume_sorted
+            other,
+            method=method,
+            kwargs=kwargs,
+            assume_sorted=assume_sorted,
+            reduce=reduce,
         )
         return self._from_temp_dataset(ds)