Added implementation for 1D linear interpolation to `regridding.regri…

…d()`.

docs/refs.bib

@@ -33,3 +33,16 @@ @article{Kumar2018
  keywords = {Axis-crossing, Computational Geometry, Polygons, Point-in-Polygon test, Winding Number},
  abstract = {This work is an extension of an axis-crossing algorithm to compute winding number for solving point in polygon for an arbitary polygon. Polygons are popular drawings in computer graphics to represent different types of structures with approximations. Solutions for point-in-polygons are many, like even-odd rule, positive-negative number, and winding number. This paper mainly deals with improvements of ‘A winding number and point in polygon algorithm’. Point in polygon is a fundamental problem and has various applications in ray tracing, computer graphics, image processing, gaming applications, robotics, acoustics, geo-science etc. The main focus of this paper explains about winding number for a closed polygon ‘S’, to test whether point ‘P’ lies either inside or outside with respect to positive and negative axis-crossing algorithm method.}
 }
+@article{Ramshaw1985,
+ title = {Conservative rezoning algorithm for generalized two-dimensional meshes},
+ journal = {Journal of Computational Physics},
+ volume = {59},
+ number = {2},
+ pages = {193-199},
+ year = {1985},
+ issn = {0021-9991},
+ doi = {https://doi.org/10.1016/0021-9991(85)90141-X},
+ url = {https://www.sciencedirect.com/science/article/pii/002199918590141X},
+ author = {John D Ramshaw},
+ abstract = {A method is presented for transferring a conserved quantity Q from one generalized mesh to another when the volumetric density of Q is uniform within each cell of the original mesh.}
+}

pyproject.toml

@@ -26,6 +26,7 @@ test = [
 ]
 doc = [
  "pytest",
+ "scipy",
  "matplotlib",
  "graphviz",
  "sphinx-autodoc-typehints",

regridding/__init__.py

@@ -1,5 +1,6 @@
 from . import math
 from . import geometry
 from ._find_indices import *
+from ._weights import *
 from ._interp_ndarray import *
 from ._regrid import *

regridding/_regrid/__init__.py

@@ -0,0 +1,2 @@
+from ._regrid_from_weights import *
+from ._regrid import *

regridding/_regrid/_regrid.py

@@ -0,0 +1,165 @@
+from typing import Sequence, Literal
+import numpy as np
+import regridding
+from . import regrid_from_weights
+
+__all__ = [
+ "regrid",
+]
+
+
+def regrid(
+ coordinates_input: tuple[np.ndarray, ...],
+ coordinates_output: tuple[np.ndarray, ...],
+ values_input: np.ndarray,
+ values_output: None | np.ndarray = None,
+ axis_input: None | int | Sequence[int] = None,
+ axis_output: None | int | Sequence[int] = None,
+ method: Literal["multilinear", "conservative"] = "multilinear",
+) -> np.ndarray:
+ """
+ Regrid an array of values defined on a logically-rectangular curvilinear
+ grid onto a new logically-rectangular curvilinear grid.
+
+ Parameters
+ ----------
+ coordinates_input
+ Coordinates of the input grid.
+ coordinates_output
+ Coordinates of the output grid.
+ Should have the same number of coordinates as the input grid.
+ values_input
+ Input array of values to be resampled.
+ values_output
+ Optional array in which to place the output.
+ axis_input
+ Logical axes of the input grid to resample.
+ If :obj:`None`, resample all the axes of the input grid.
+ The number of axes should be equal to the number of
+ coordinates in the input grid.
+ axis_output
+ Logical axes of the output grid corresponding to the resampled axes
+ 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.
+ method
+ The type of regridding to use.
+ The ``conservative`` method uses the algorithm described in
+ :footcite:t:`Ramshaw1985`
+
+ See Also
+ --------
+ :func:`regridding.regrid`
+ :func:`regridding.regrid_from_weights`
+
+ Examples
+ --------
+
+ Regrid a 1D array using multilinear interpolation.
+
+ .. jupyter-execute::
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+ import regridding
+
+ # Define the input grid
+ x_input = np.linspace(-1, 1, num=11)
+
+ # Define the input array
+ values_input = np.square(x_input)
+
+ # Define the output grid
+ x_output = np.linspace(-1, 1, num=51)
+
+ # Regrid the input array onto the output grid
+ values_output = regridding.regrid(
+ coordinates_input=(x_input,),
+ coordinates_output=(x_output,),
+ values_input=values_input,
+ method="multilinear",
+ )
+
+ # Plot the results
+ plt.figure(figsize=(6, 3));
+ plt.scatter(x_input, values_input, s=100, label="input", zorder=1);
+ plt.scatter(x_output, values_output, label="interpolated", zorder=0);
+ plt.legend();
+
+ |
+
+ Regrid a 2D array using conservative resampling.
+
+ .. jupyter-execute::
+
+ # Define the number of edges in the input grid
+ num_x = 66
+ num_y = 66
+
+ # Define a dummy linear grid
+ x = np.linspace(-5, 5, num=num_x)
+ y = np.linspace(-5, 5, num=num_y)
+ x, y = np.meshgrid(x, y, indexing="ij")
+
+ # Define the curvilinear input grid using the dummy grid
+ angle = 0.4
+ x_input = x * np.cos(angle) - y * np.sin(angle) + 0.05 * x * x
+ y_input = x * np.sin(angle) + y * np.cos(angle) + 0.05 * y * y
+
+ # Define the test pattern
+ pitch = 16
+ a_input = 0 * x[:~0,:~0]
+ a_input[::pitch, :] = 1
+ a_input[:, ::pitch] = 1
+ a_input[pitch//2::pitch, pitch//2::pitch] = 1
+
+ # Define a rectilinear output grid using the limits of the input grid
+ x_output = np.linspace(x_input.min(), x_input.max(), num_x // 2)
+ y_output = np.linspace(y_input.min(), y_input.max(), num_y // 2)
+ x_output, y_output = np.meshgrid(x_output, y_output, indexing="ij")
+
+ # Regrid the test pattern onto the new grid
+ a_output = regridding.regrid(
+ coordinates_input=(x_input, y_input),
+ coordinates_output=(x_output, y_output),
+ values_input=a_input,
+ method="conservative",
+ )
+
+ fig, axs = plt.subplots(
+ ncols=2,
+ sharex=True,
+ sharey=True,
+ figsize=(8, 4),
+ constrained_layout=True,
+ );
+ axs[0].pcolormesh(x_input, y_input, a_input);
+ axs[0].set_title("input array");
+ axs[1].pcolormesh(x_output, y_output, a_output);
+ axs[1].set_title("regridded array");
+
+ |
+
+ References
+ ----------
+ .. footbibliography::
+ """
+ weights, shape_input, shape_output = regridding.weights(
+ coordinates_input=coordinates_input,
+ coordinates_output=coordinates_output,
+ axis_input=axis_input,
+ axis_output=axis_output,
+ method=method,
+ )
+ result = regrid_from_weights(
+ weights=weights,
+ shape_input=shape_input,
+ shape_output=shape_output,
+ values_input=values_input,
+ values_output=values_output,
+ axis_input=axis_input,
+ axis_output=axis_output,
+ )
+ return result

regridding/_regrid/_regrid_from_weights.py

@@ -0,0 +1,112 @@
+from typing import Sequence
+import numpy as np
+import numba
+from regridding import _util
+
+__all__ = [
+ "regrid_from_weights",
+]
+
+
+def regrid_from_weights(
+ weights: np.ndarray,
+ shape_input: tuple[int, ...],
+ shape_output: tuple[int, ...],
+ values_input: np.ndarray,
+ values_output: None | np.ndarray = None,
+ axis_input: None | int | Sequence[int] = None,
+ axis_output: None | int | Sequence[int] = None,
+) -> np.ndarray:
+ """
+ Regrid an array of values using weights computed by
+ :func:`regridding.weights`.
+
+ Parameters
+ ----------
+ weights
+ Ragged array of weights computed by :func:`regridding.weights`.
+ shape_input
+ Broadcasted shape of the input coordinates computed by :func:`regridding.weights`.
+ shape_output
+ Broadcasted shape of the output coordinates computed by :func:`regridding.weights`.
+ values_input
+ Input array of values to be resampled.
+ values_output
+ Optional array in which to place the output.
+ axis_input
+ Logical axes of the input array to resample.
+ If :obj:`None`, resample all the axes of the input array.
+ The number of axes should be equal to the number of
+ coordinates in the original input grid passed to :func:`regridding.weights`.
+ axis_output
+ Logical axes of the output array corresponding to the resampled axes
+ 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`.
+
+ See Also
+ --------
+ :func:`regridding.regrid`
+ :func:`regridding.regrid_from_weights`
+ """
+
+ shape_input = np.broadcast_shapes(shape_input, values_input.shape)
+ values_input = np.broadcast_to(values_input, shape=shape_input, subok=True)
+ ndim_input = len(shape_input)
+ axis_input = _util._normalize_axis(axis_input, ndim=ndim_input)
+
+ if values_output is None:
+ shape_output = np.broadcast_shapes(
+ shape_output,
+ tuple(
+ shape_input[ax] if ax not in axis_input else 1
+ for ax in _util._normalize_axis(None, ndim_input)
+ ),
+ )
+ values_output = np.zeros_like(values_input, shape=shape_output)
+ else:
+ if values_output.shape != shape_output:
+ raise ValueError(f"")
+ values_output.fill(0)
+ values_output_copy = values_output
+ ndim_output = len(shape_output)
+ axis_output = _util._normalize_axis(axis_output, ndim=ndim_output)
+
+ axis_input_numba = ~np.arange(len(axis_input))[::-1]
+ axis_output_numba = ~np.arange(len(axis_output))[::-1]
+
+ shape_input_numba = tuple(shape_input[ax] for ax in axis_input)
+ shape_output_numba = tuple(shape_output[ax] for ax in axis_output)
+
+ values_input = np.moveaxis(values_input, axis_input, axis_input_numba)
+ values_output = np.moveaxis(values_output, axis_output, axis_output_numba)
+
+ weights = numba.typed.List(weights.reshape(-1))
+ values_input = values_input.reshape(-1, *shape_input_numba)
+ values_output = values_output.reshape(-1, *shape_output_numba)
+
+ _regrid_from_weights(
+ weights=weights,
+ values_input=values_input,
+ values_output=values_output,
+ )
+
+ return values_output_copy
+
+
+@numba.njit()
+def _regrid_from_weights(
+ weights: numba.typed.List,
+ values_input: np.ndarray,
+ values_output: np.ndarray,
+) -> None:
+ for d in numba.prange(len(weights)):
+ weights_d = weights[d]
+ values_input_d = values_input[d].reshape(-1)
+ values_output_d = values_output[d].reshape(-1)
+
+ for w in range(len(weights_d)):
+ i_input, i_output, weight = weights_d[w]
+ values_output_d[int(i_output)] += weight * values_input_d[int(i_input)]