Extend method="multilinear" to grids of any number of dimensions and to scattered points

README.md

@@ -39,7 +39,9 @@ pip install regridding
 
 * [`regrid()`](https://regridding.readthedocs.io/en/latest/_autosummary/regridding.regrid.html),
   which resamples an array onto a new grid using either of two methods:
-  * `"multilinear"`, linear interpolation along one axis;
+  * `"multilinear"`, linear interpolation on logically-rectangular grids of
+    any number of dimensions, rectilinear or curvilinear, at points which may
+    be arranged along any number of axes;
   * `"conservative"`, first-order conservative resampling of 1D grids and of 2D
     logically-rectangular curvilinear grids, using the algorithm described in
     [Ramshaw (1985)](https://doi.org/10.1016/0021-9991(85)90141-X).
@@ -86,6 +88,9 @@ The remaining axes are orthogonal to the operation, and the resampling is
 repeated independently for every position along them.
 This is how a stack of images, or a spectrum for each pixel, is resampled in one
 call.
+The `multilinear` method interpolates at points rather than resampling cells,
+so an explicit `axis_output` may select any number of axes, such as a single
+axis of scattered points to interpolate a 2D grid at.
 
 **Degenerate grids are perturbed.**
 Where a vertex of the output grid lands exactly on an edge of the input grid,

benchmarks/weights_multilinear.py

@@ -0,0 +1,67 @@
+import numpy as np
+import regridding
+
+coordinates_input: tuple[np.ndarray, ...] = ()
+coordinates_output: tuple[np.ndarray, ...] = ()
+
+
+def setup_weights_multilinear_2d(
+    num_point: int,
+    grid: str,
+) -> None:
+    """
+    A grid of about a million vertices, and random points over it.
+
+    The ``rectilinear`` grid is searched along each axis on its own, and the
+    ``curvilinear`` one, which is the same grid rotated and distorted,
+    through its hierarchy of bounding boxes.
+    """
+    global coordinates_input
+    global coordinates_output
+
+    x, y = np.meshgrid(
+        np.linspace(-1, 1, num=1001),
+        np.linspace(-1, 1, num=1000),
+        indexing="ij",
+    )
+
+    if grid == "curvilinear":
+        angle = 0.4
+        x, y = (
+            x * np.cos(angle) - y * np.sin(angle) + 0.05 * x * x,
+            x * np.sin(angle) + y * np.cos(angle) + 0.05 * y * y,
+        )
+
+    rng = np.random.default_rng(0)
+    coordinates_input = (x, y)
+    coordinates_output = (
+        rng.uniform(x.min(), x.max(), size=num_point),
+        rng.uniform(y.min(), y.max(), size=num_point),
+    )
+
+    time_weights_multilinear_2d(num_point, grid)
+
+
+def time_weights_multilinear_2d(
+    num_point: int,
+    grid: str,
+) -> None:
+    regridding.weights(
+        coordinates_input=coordinates_input,
+        coordinates_output=coordinates_output,
+        axis_output=0,
+        method="multilinear",
+        coalesce=False,
+    )
+
+
+time_weights_multilinear_2d.setup = setup_weights_multilinear_2d
+time_weights_multilinear_2d.params = (
+    [1, 10_000, 1_000_000],
+    ["rectilinear", "curvilinear"],
+)
+time_weights_multilinear_2d.param_names = (
+    "number of points",
+    "grid",
+)
+time_weights_multilinear_2d.timeout = 240

docs/index.rst

@@ -36,7 +36,9 @@ Features
 *   :func:`regridding.regrid`, which resamples an array onto a new grid using
     either of two methods:
 
-    *   ``"multilinear"``, linear interpolation along one axis;
+    *   ``"multilinear"``, linear interpolation on logically-rectangular grids
+        of any number of dimensions, rectilinear or curvilinear, at points
+        which may be arranged along any number of axes;
     *   ``"conservative"``, first-order conservative resampling of 1D grids and
         of 2D logically-rectangular curvilinear grids, by one of two
         algorithms described below.
@@ -98,6 +100,9 @@ The remaining axes are orthogonal to the operation, and the resampling is
 repeated independently for every position along them.
 This is how a stack of images, or a spectrum for each pixel, is resampled in
 one call.
+The ``multilinear`` method interpolates at points rather than resampling cells,
+so an explicit ``axis_output`` may select any number of axes, such as a single
+axis of scattered points to interpolate a 2D grid at.
 
 **The conservative method has two algorithms, and picks between them.**
 Both compute the same weights: the fraction of each input cell shared with

regridding/_regrid/_regrid.py

@@ -47,12 +47,18 @@ def regrid(
         of the input grid.
         If :obj:`None`, all the axes of the output grid correspond to resampled
         axes in the input grid.
-        The number of axes should be equal to the number of
-        coordinates in the output grid.
+        For the ``conservative`` method, the number of axes should be equal
+        to the number of coordinates in the output grid, as it should be for
+        the ``multilinear`` method if `axis_output` is :obj:`None`.
+        The ``multilinear`` method interpolates at points rather than
+        resampling cells, so if `axis_output` is given, the points may be
+        arranged along any number of axes, or none at all for a single point.
     method
         The type of regridding to use.
         The ``multilinear`` method interprets `coordinates_input` as the
-        points where `values_input` is sampled.
+        points where `values_input` is sampled, and interpolates between the
+        corners of the cell which contains each output point, on a grid of
+        any number of dimensions, rectilinear or curvilinear.
         The ``conservative`` method interprets `coordinates_input` as the edges
         of the cells containing `values_input`, so `values_input` has one fewer
         element along each resampled axis, and the sum of the result matches the
@@ -69,6 +75,12 @@ def regrid(
         ``extrapolate`` (the default) extrapolates linearly from the nearest
         cell of the input grid, ``nan`` sets those points to
         :obj:`numpy.nan`, and ``raise`` raises a :class:`ValueError`.
+        On a grid which is not rectilinear, the nearest cell is the one along
+        the edge of the grid whose continuation reaches the point, and a
+        point which no continuation reaches, as can happen beyond a strongly
+        curved edge, is set to :obj:`numpy.nan`.
+        A point with a coordinate which is not finite is set to
+        :obj:`numpy.nan`, or raises, whichever `bounds` is.
     perturb
         Whether to perturb `coordinates_output` by a small value to avoid degenerate
         grids. This is helpful for some methods, like ``conservative``, which
@@ -139,6 +151,43 @@ def regrid(
         ax.stairs(values_output, x_output, label="output");
         ax.legend();
 
+    Interpolate an array sampled on a distorted 2D grid at scattered points.
+
+    .. jupyter-execute::
+
+        # Define the vertices of the input grid
+        u, v = np.meshgrid(
+            np.linspace(-1, 1, num=21),
+            np.linspace(-1, 1, num=21),
+            indexing="ij",
+        )
+        x_input = u + 0.2 * v * v
+        y_input = v - 0.2 * u * u
+
+        # Sample a function on the vertices
+        values_input = np.cos(3 * x_input) * np.sin(2 * y_input)
+
+        # Define a set of random points, some of them outside the grid
+        rng = np.random.default_rng(0)
+        x_output = rng.uniform(-1.3, 1.3, size=500)
+        y_output = rng.uniform(-1.3, 1.3, size=500)
+
+        # Interpolate the function at the points
+        values_output = regridding.regrid(
+            coordinates_input=(x_input, y_input),
+            coordinates_output=(x_output, y_output),
+            values_input=values_input,
+            axis_output=0,
+            method="multilinear",
+            bounds="nan",
+        )
+
+        # Plot the input grid and the interpolated points
+        fig, ax = plt.subplots(constrained_layout=True);
+        ax.pcolormesh(x_input, y_input, values_input, shading="gouraud", vmin=-1, vmax=1);
+        ax.scatter(x_output, y_output, c=values_output, vmin=-1, vmax=1, edgecolors="black");
+        ax.set_aspect("equal");
+
     References
     ----------
     .. footbibliography::

regridding/_regrid/_regrid_from_weights.py

@@ -51,8 +51,9 @@ def regrid_from_weights(
         of the input array.
         If :obj:`None`, all the axes of the output array correspond to resampled
         axes in the input grid.
-        The number of axes should be equal to the original number of
-        coordinates in the output grid passed to :func:`regridding.weights`.
+        These should be the axes passed as `axis_output` to
+        :func:`regridding.weights`, which for the ``multilinear`` method may
+        be any number of axes, or none for a single point.
 
     Notes
     -----
@@ -137,8 +138,17 @@ def regrid_from_weights(
 
     shape_output_tmp = values_output.shape
 
-    values_input = values_input.reshape(-1, *shape_input_numba)
-    values_output = values_output.reshape(-1, *shape_output_numba)
+    # The weights address each element by a flat index, so its resampled
+    # axes are flattened together.  That also covers a single point, which
+    # has no resampled axes at all.  The number of elements is given rather
+    # than inferred, since it cannot be when there are no points.
+    num_orthogonal = int(np.prod(shape_orthogonal, dtype=int))
+    values_input = values_input.reshape(
+        num_orthogonal, np.prod(shape_input_numba, dtype=int)
+    )
+    values_output = values_output.reshape(
+        num_orthogonal, np.prod(shape_output_numba, dtype=int)
+    )
 
     flat = weights.reshape(-1)
     unit_weights = getattr(flat[0][2], "unit", None) if flat.size else None
@@ -152,13 +162,16 @@ def regrid_from_weights(
     values_input = np.ascontiguousarray(values_input)
     values_output = np.ascontiguousarray(values_output)
 
-    _regrid_from_weights(
-        weights=weights_numba,
-        values_input=values_input,
-        values_output=values_output,
-    )
+    # an empty typed list has no type for the kernel to be compiled for, and
+    # with no elements there is nothing to add to the zeros of the output
+    if flat.size:
+        _regrid_from_weights(
+            weights=weights_numba,
+            values_input=values_input,
+            values_output=values_output,
+        )
 
-    values_output = values_output.reshape(*shape_output_tmp)
+    values_output = values_output.reshape(shape_output_tmp)
 
     values_output = np.moveaxis(values_output, axis_output_numba, axis_output)
 

regridding/_util.py

@@ -110,7 +110,7 @@ def _normalize_axis(
     if axis is None:
         axis = tuple(range(ndim))
     axis = np.lib.array_utils.normalize_axis_tuple(axis, ndim=ndim)
-    axis = tuple(~(~np.array(axis) % ndim))
+    axis = tuple(~(~np.array(axis, dtype=int) % ndim))
     return axis
 
 
@@ -121,6 +121,7 @@ def _normalize_input_output_coordinates(
     axis_output: None | int | Sequence[int] = None,
     perturb: bool = False,
     seed: "None | int | np.random.Generator" = _seed_default,
+    match_axes: bool = True,
 ) -> tuple[
     tuple[np.ndarray, ...],
     tuple[np.ndarray, ...],
@@ -130,6 +131,30 @@ def _normalize_input_output_coordinates(
     tuple[int, ...],
     tuple[int, ...],
 ]:
+    """
+    Broadcast the coordinates of an input and an output grid against each
+    other, apart from the axes which are resampled.
+
+    Parameters
+    ----------
+    coordinates_input
+        The coordinates of the input grid.
+    coordinates_output
+        The coordinates of the output grid.
+    axis_input
+        The resampled axes of the input grid.
+    axis_output
+        The resampled axes of the output grid.
+    perturb
+        Whether to perturb the output coordinates by a small random amount.
+    seed
+        The seed of the random perturbation.
+    match_axes
+        Whether `axis_output` must have as many axes as `axis_input`.
+        Multilinear interpolation does not need it, since its outputs are
+        points rather than cells, and they can be arranged in any number of
+        axes.
+    """
     if isinstance(coordinates_input, np.ndarray):
         coordinates_input = (coordinates_input,)
 
@@ -159,10 +184,12 @@ def _normalize_input_output_coordinates(
     axis_input = tuple(sorted(axis_input, reverse=True))
     axis_output = tuple(sorted(axis_output, reverse=True))
 
-    if len(axis_output) != len(axis_input):
+    if match_axes and len(axis_output) != len(axis_input):
         raise ValueError(
-            f"The number of axes in `axis_output`, {axis_output}, "
-            f"must match the number of axes in `axis_input`, {axis_input}"
+            f"The number of axes in `axis_output`, "
+            f"{tuple(int(ax) for ax in axis_output)}, "
+            f"must match the number of axes in `axis_input`, "
+            f"{tuple(int(ax) for ax in axis_input)}"
         )
 
     if len(coordinates_input) != len(axis_input):

regridding/_weights/_weights.py

@@ -1,7 +1,7 @@
 from typing import Sequence, Literal
 import numpy as np
 from regridding import _util
-from ._weights_multilinear import _weights_multilinear
+from ._weights_multilinear import weights_multilinear
 from ._weights_conservative import _weights_conservative
 from ._weights_arrays import _weights_to_arrays
 
@@ -141,19 +141,34 @@ def weights(
         of the input grid.
         If :obj:`None`, all the axes of the output grid correspond to resampled
         axes in the input grid.
-        The number of axes should be equal to the number of
-        coordinates in the output grid.
+        For the ``conservative`` method, the number of axes should be equal
+        to the number of coordinates in the output grid, as it should be for
+        the ``multilinear`` method if `axis_output` is :obj:`None`.
+        The ``multilinear`` method interpolates at points rather than
+        resampling cells, so if `axis_output` is given, the points may be
+        arranged along any number of axes, or none at all for a single point.
     weights_input
         Weights applied to the values of the input grid before resampling.
     method
         The type of regridding to use.
+        The ``multilinear`` method interpolates between the corners of the
+        cell of the input grid which contains each output point, on a grid of
+        any number of dimensions, rectilinear or curvilinear.
+        The ``conservative`` method resamples the cells of the input grid
+        onto the cells of the output grid, in one or two dimensions.
     bounds
         How to treat output points that fall outside the input grid.
         Only applies when `method` is ``multilinear``; the ``conservative``
         method assigns no weight to regions outside the input grid.
         ``extrapolate`` (the default) extrapolates linearly from the nearest
         cell of the input grid, ``nan`` sets those points to
         :obj:`numpy.nan`, and ``raise`` raises a :class:`ValueError`.
+        On a grid which is not rectilinear, the nearest cell is the one along
+        the edge of the grid whose continuation reaches the point, and a
+        point which no continuation reaches, as can happen beyond a strongly
+        curved edge, is set to :obj:`numpy.nan`.
+        A point with a coordinate which is not finite is set to
+        :obj:`numpy.nan`, or raises, whichever `bounds` is.
     perturb
         Whether to perturb `coordinates_output` by a small value to avoid degenerate
         grids. This is helpful for some methods, like ``conservative``, which
@@ -339,7 +354,7 @@ def weights(
     clipped = False
 
     if method == "multilinear":
-        result = _weights_multilinear(
+        result = weights_multilinear(
             coordinates_input=coordinates_input,
             coordinates_output=coordinates_output,
             axis_input=axis_input,