Interpolator Objects

This example explains how to create the fundamental objects required for interpolation with pyinterp: the axes and the grid. These objects are the building blocks for defining the space and data on which interpolation will be performed.

First, let’s import the necessary libraries and load a sample dataset.

import numpy as np

import pyinterp
import pyinterp.backends.xarray
import pyinterp.tests

# Load a 3D test grid from the library's test data
ds = pyinterp.tests.load_grid3d()
lon, lat, time, tcw = (
    ds["longitude"].values,
    ds["latitude"].values,
    ds["time"].values,
    ds["tcw"].values,
)

Defining Axes

An axis defines the coordinates for one dimension of the grid. pyinterp provides specialized axis objects for performance and to handle specific cases like circular (e.g., longitude) and temporal axes.

Regular Axis

A regular axis is a simple, monotonically increasing or decreasing array of coordinates. Here, we create an axis for latitude.

y_axis = pyinterp.Axis(lat)
print("Latitude axis:")
print(y_axis)

Latitude axis:
Axis(regular)
  range: [-90, 90]
  step: -2.5
  size: 73

You can use the pyinterp.Axis.find_index() method to find the nearest grid index for a given coordinate.

print(f"Index for latitude 0.12°: {y_axis.find_index(np.array([0.12]))}")

Index for latitude 0.12°: [36]

Periodic Axis (Longitude)

For axes that wrap around, like longitude, you can create a “periodic” axis. This ensures that coordinates are correctly handled at the boundary (e.g., -180° and 180° are treated as the same point).

x_axis = pyinterp.Axis(lon, period=360.0)
print("Longitude axis (periodic):")
print(x_axis)

Longitude axis (periodic):
Axis(regular, period=360)
  range: [0, 357.5]
  step: 2.5
  size: 144

With a periodic axis, boundary points are correctly identified as identical.

idx_1 = x_axis.find_index(np.array([-180]))
idx_2 = x_axis.find_index(np.array([180]))
print(f"Are indices for -180° and 180° the same? {idx_1 == idx_2}")

Are indices for -180° and 180° the same? [ True]

Temporal Axis

For time coordinates, pyinterp provides a highly optimized pyinterp.TemporalAxis class.

t_axis = pyinterp.TemporalAxis(time)
print("Time axis:")
print(t_axis)

Time axis:
TemporalAxis[datetime64[ns]](regular)
  range: [2002-07-01T12:00:00.000000000, 2002-07-31T12:00:00.000000000]
  step: 86400000000000 nanoseconds
  size: 31

Creating a Grid

Once the axes are defined, you can create a grid object that holds the data. The grid object takes the axes and the data array as input. The dimensions of the data array must follow the order of the axes passed to the constructor: here (longitude, latitude, time).

Here, we create a pyinterp.backends.xarray.Grid3D object for the total column water vapor (TCW) data. The TCW values are stored in the dataset as (time, latitude, longitude), so we transpose them to obtain an array ordered as (longitude, latitude, time).

tcw = tcw.T
grid = pyinterp.Grid(x_axis, y_axis, t_axis, tcw)
print("Grid object:")
print(grid)

Grid object:
TemporalGrid3D(shape=(144, 73, 31), dtype=float64, nbytes=2 MB)

Using the XArray Backend

For convenience, pyinterp provides a backend that can directly create a grid interpolator from an xarray.DataArray. This avoids the need to manually create the axes and grid objects.

The backend automatically detects the axis types (regular, circular, temporal) and creates the appropriate grid interpolator.

interpolator = pyinterp.backends.xarray.Grid3D(ds["tcw"])
print("Interpolator from xarray.DataArray:")
print(interpolator)

Interpolator from xarray.DataArray:
TemporalGrid3D(shape=(144, 73, 31), dtype=float64, nbytes=2 MB)