Fast interpolation of geo-referenced data

pyinterp#

Fast interpolation of geo-referenced data for Python and xarray.

pyinterp interpolates gridded fields in 2, 3 or 4 dimensions and scattered observations, such as model outputs, satellite tracks or in situ measurements. Its core is written in C++ and runs on all the cores of the machine; it knows that longitudes wrap around the globe, that a time axis holds dates, and that a field may have undefined values, near the coasts for instance.

Installation#

pip install "pyinterp[xarray]"

or, with conda:

conda install -c conda-forge pyinterp

Binary wheels are published for Linux, macOS (Apple Silicon) and Windows, for Python 3.12 and later. The xarray extra is needed only for the xarray interface; see the installation guide.

Quick start#

Once pyinterp is imported, xarray data arrays have a pyinterp accessor that works like DataArray.interp:

import numpy as np
import xarray as xr

import pyinterp
import pyinterp.tests

# A global mean sea surface, stored as (lat, lon) with longitudes in [0, 360)
ds = pyinterp.tests.load_grid2d()

# Interpolate on a new grid: the longitude axis is periodic, so the grid can
# be queried in [-180, 180)
mss = ds.mss.pyinterp.interp(
    lon=np.arange(-180, 180, 0.1),
    lat=np.arange(-80, 80, 0.1),
    method="bicubic",
)

# Interpolate point by point, along a track
track = ds.mss.pyinterp.interp(
    lon=xr.DataArray(np.linspace(-170, 170, 1000), dims="obs"),
    lat=xr.DataArray(np.linspace(-60, 60, 1000), dims="obs"),
)

Scattered observations are indexed in an R*-tree:

from pyinterp.geometry import geographic

rng = np.random.default_rng(0)
lon, lat = rng.uniform(-180, 180, 10_000), rng.uniform(-80, 80, 10_000)
values = np.cos(np.radians(lat)) * np.sin(np.radians(lon))

# Index the observations on the WGS-84 ellipsoid
tree = pyinterp.RTree3D(spheroid=geographic.Spheroid())
tree.packing(np.column_stack((lon, lat)), values)

# Estimate the field at new positions from the 11 nearest observations
estimate, neighbors = pyinterp.inverse_distance_weighting(
    tree, np.array([[10.0, 45.0], [-70.0, -30.0]]), k=11
)

The gallery shows these and the other features on real data.

Why pyinterp?#

xarray.DataArray.interp and scipy.interpolate are general-purpose tools. pyinterp is built for geo-referenced data:

Periodic longitudes: a grid stored in [0, 360) is interpolated in [-180, 180), and across the antimeridian, without rolling the data.
Dates: time axes of type datetime64 are handled with integer arithmetic, at any resolution, without conversion to floats.
Speed: the interpolators are written in C++ and process the points in parallel; the axes find the cell of a point by a direct computation when they are regular. Bilinear interpolation runs 5 times faster than SciPy on one core and 38 times faster on all the cores; see the benchmarks.
Undefined values: the gaps of a field, such as the land of an ocean grid, can be filled beforehand, so that the windowed methods (bicubic, splines, Akima, Steffen…) do not spread them; at the edges of the grid, their window can shrink to the available data.
Beyond grids: scattered data, statistics and geometry use the same coordinates and the same conventions.

What is inside#

Regular grids: 1D to 4D interpolation (nearest, bilinear, inverse distance weighting, bicubic, splines, Akima, Steffen, polynomial), with temporal axes, through xarray or NumPy.
Conservative regridding: area-weighted averages between rectilinear grids, which preserve the integral of the field, on the sphere or on the ellipsoid.
Reduced Gaussian grids: interpolation of the ECMWF N and O grids in 2D, 3D and time, reading only the part of the field the points need.
Scattered data: R*-trees in Cartesian or geodetic coordinates, with inverse distance weighting, radial basis functions, kriging and window functions, and optimal interpolation.
Gap filling: LOESS, Gauss-Seidel relaxation, multigrid and FFT inpainting.
Statistics: 1D and 2D binning, histograms, descriptive statistics and T-Digest quantiles, which can be merged across Dask chunks.
Geometry: bindings to Boost.Geometry in Cartesian and geographic coordinates, geohash encoding, and orbit utilities for satellite tracks.

The design notes explain how the axes, the R*-trees and the coordinate conversions work.

Contributing and citing#

pyinterp is developed by CNES and CLS under the BSD-3-Clause license. Bug reports, questions and pull requests are welcome on GitHub; see the contributing guide. If pyinterp helps your research, please cite it using the metadata of the CITATION.cff file.

pyinterp

Contents