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 timeit

import numpy
import pandas

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:
<pyinterp.core.Axis>
  min_value: -90
  max_value: 90
  step     : -2.5
  is_circle: false

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([0.12])}')

Index for latitude 0.12°: [36]

Circular Axis (Longitude)

For axes that wrap around, like longitude, you can create a “circular” 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, is_circle=True)
print('\nLongitude axis (circular):')
print(x_axis)

Longitude axis (circular):
<pyinterp.core.Axis>
  min_value: 0
  max_value: 357.5
  step     : 2.5
  is_circle: true

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

print('Are indices for -180° and 180° the same? '
      f'{x_axis.find_index([-180]) == x_axis.find_index([180])}')

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:
<pyinterp.core.TemporalAxis>
  min_value: 2002-07-01T12:00:00.000000000
  max_value: 2002-07-31T12:00:00.000000000
  step     : 86400000000000 nanoseconds

Performance Comparison

pyinterp axis objects are implemented in C++ and are significantly faster for lookups than equivalent objects in libraries like pandas. Let’s compare the performance of pyinterp.Axis against pandas.Index for finding indices.

Longitude Axis:

values = lon[10:20] + 1 / 3
index = pandas.Index(lon)
print('Performance for Longitude Axis:')
print(f'  pandas.Index:  '
      f"{timeit.timeit('index.searchsorted(values)', globals=globals()):.6f}s")
print(f'  pyinterp.Axis: '
      f"{timeit.timeit('x_axis.find_index(values)', globals=globals()):.6f}s")

Performance for Longitude Axis:
  pandas.Index:  1.364317s
  pyinterp.Axis: 0.848444s

Time Axis:

index = pandas.Index(time)
values = time + numpy.timedelta64(1, 'ns')
print('Performance for Time Axis:')
print(f'  pandas.Index:        '
      f"{timeit.timeit('index.searchsorted(values)', globals=globals()):.6f}s")
print(f'  pyinterp.TemporalAxis: '
      f"{timeit.timeit('t_axis.find_index(values)', globals=globals()):.6f}s")

Performance for Time Axis:
  pandas.Index:        19.485896s
  pyinterp.TemporalAxis: 4.704602s

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 data must be organized to match the order of the axes: (time, latitude, longitude).

Here, we create a pyinterp.Grid3D object for the total column water vapor (TCW) data. Before constructing the tensor for pyinterp, we must begin to organize the data in a grid with the values of the axes in the first dimensions of the tensor.

In our case, the time, latitude, longitude axes must be sorted in this order.

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

Grid object:
<pyinterp.grid.Grid3D>
array([[[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [12.08470387,  9.64372637, 12.96766807, ..., 10.47133655,
         12.06850269, 24.25583891],
        [15.77317208, 10.12436132, 11.59731843, ...,  9.53976881,
         11.00327523, 21.59209521],
        ...,
        [ 0.38745328,  0.44010711,  0.65747292, ...,  0.19438924,
          0.57376683,  0.36045132],
        [ 0.39420377,  0.35910122,  0.57106663, ...,  0.18763875,
          0.3739523 ,  0.27269494],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [12.19541192,  9.65722735, 12.9460665 , ..., 10.34847761,
         12.54643744, 24.19508449],
        [15.72996894, 10.18511573, 11.41775537, ...,  9.63562578,
         11.89839032, 21.96202211],
        ...,
        [ 0.38340299,  0.43065642,  0.64937233, ...,  0.19573934,
          0.56026585,  0.34695034],
        [ 0.38880338,  0.35505093,  0.56566624, ...,  0.18763875,
          0.36855191,  0.26729454],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [12.30611997,  9.67072833, 12.92311483, ..., 10.22561868,
         13.02707239, 24.13433007],
        [15.68946599, 10.24316996, 11.23954241, ...,  9.73283285,
         12.79485551, 22.33194901],
        ...,
        [ 0.37935269,  0.41985564,  0.63992164, ...,  0.19978964,
          0.54406467,  0.33344936],
        [ 0.38475309,  0.34830043,  0.56026585, ...,  0.18763875,
          0.36315151,  0.26324425],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       ...,

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [11.72422765,  9.73148275, 12.9447164 , ..., 10.88716679,
         11.01542611, 24.1167788 ],
        [15.4585992 , 10.30257428, 11.62972078, ...,  9.36560614,
          9.72743245, 19.62230197],
        ...,
        [ 0.40500456,  0.47250947,  0.71417704, ...,  0.21599081,
          0.62642066,  0.40770475],
        [ 0.40635466,  0.3739523 ,  0.57376683, ...,  0.21329062,
          0.38880338,  0.28889612],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [11.84438639,  9.70178059, 12.95281699, ..., 10.74810667,
         11.36645164, 24.16403223],
        [15.56390686, 10.24316996, 11.61892   , ...,  9.42501046,
         10.15406348, 20.27844969],
        ...,
        [ 0.39825407,  0.46170868,  0.69527567, ...,  0.20789022,
          0.60886938,  0.39285367],
        [ 0.40230436,  0.36855191,  0.57376683, ...,  0.20383993,
          0.38475309,  0.28484582],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [11.96454513,  9.67342853, 12.96091758, ..., 10.61039666,
         11.71612706, 24.20858547],
        [15.66786442, 10.18241554, 11.60676911, ...,  9.48306468,
         10.5779943 , 20.9345974 ],
        ...,
        [ 0.39285367,  0.45225799,  0.67502419, ...,  0.20248983,
          0.59131811,  0.3766525 ],
        [ 0.39825407,  0.36450161,  0.57376683, ...,  0.19573934,
          0.37800259,  0.27944543],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]]], shape=(144, 73, 31))
Axis:
 * x: <pyinterp.core.Axis>
       min_value: 0
       max_value: 357.5
       step     : 2.5
       is_circle: true
 * y: <pyinterp.core.Axis>
       min_value: -90
       max_value: 90
       step     : -2.5
       is_circle: false
 * z: <pyinterp.core.TemporalAxis>
       min_value: 2002-07-01T12:00:00.000000000
       max_value: 2002-07-31T12:00:00.000000000
       step     : 86400000000000 nanoseconds

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:
<pyinterp.backends.xarray.Grid3D>
array([[[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [12.08470387,  9.64372637, 12.96766807, ..., 10.47133655,
         12.06850269, 24.25583891],
        [15.77317208, 10.12436132, 11.59731843, ...,  9.53976881,
         11.00327523, 21.59209521],
        ...,
        [ 0.38745328,  0.44010711,  0.65747292, ...,  0.19438924,
          0.57376683,  0.36045132],
        [ 0.39420377,  0.35910122,  0.57106663, ...,  0.18763875,
          0.3739523 ,  0.27269494],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [12.19541192,  9.65722735, 12.9460665 , ..., 10.34847761,
         12.54643744, 24.19508449],
        [15.72996894, 10.18511573, 11.41775537, ...,  9.63562578,
         11.89839032, 21.96202211],
        ...,
        [ 0.38340299,  0.43065642,  0.64937233, ...,  0.19573934,
          0.56026585,  0.34695034],
        [ 0.38880338,  0.35505093,  0.56566624, ...,  0.18763875,
          0.36855191,  0.26729454],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [12.30611997,  9.67072833, 12.92311483, ..., 10.22561868,
         13.02707239, 24.13433007],
        [15.68946599, 10.24316996, 11.23954241, ...,  9.73283285,
         12.79485551, 22.33194901],
        ...,
        [ 0.37935269,  0.41985564,  0.63992164, ...,  0.19978964,
          0.54406467,  0.33344936],
        [ 0.38475309,  0.34830043,  0.56026585, ...,  0.18763875,
          0.36315151,  0.26324425],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       ...,

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [11.72422765,  9.73148275, 12.9447164 , ..., 10.88716679,
         11.01542611, 24.1167788 ],
        [15.4585992 , 10.30257428, 11.62972078, ...,  9.36560614,
          9.72743245, 19.62230197],
        ...,
        [ 0.40500456,  0.47250947,  0.71417704, ...,  0.21599081,
          0.62642066,  0.40770475],
        [ 0.40635466,  0.3739523 ,  0.57376683, ...,  0.21329062,
          0.38880338,  0.28889612],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [11.84438639,  9.70178059, 12.95281699, ..., 10.74810667,
         11.36645164, 24.16403223],
        [15.56390686, 10.24316996, 11.61892   , ...,  9.42501046,
         10.15406348, 20.27844969],
        ...,
        [ 0.39825407,  0.46170868,  0.69527567, ...,  0.20789022,
          0.60886938,  0.39285367],
        [ 0.40230436,  0.36855191,  0.57376683, ...,  0.20383993,
          0.38475309,  0.28484582],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]],

       [[10.15271338, 10.66710078, 11.29219624, ..., 10.60634637,
         15.88388013, 19.57369844],
        [11.96454513,  9.67342853, 12.96091758, ..., 10.61039666,
         11.71612706, 24.20858547],
        [15.66786442, 10.18241554, 11.60676911, ...,  9.48306468,
         10.5779943 , 20.9345974 ],
        ...,
        [ 0.39285367,  0.45225799,  0.67502419, ...,  0.20248983,
          0.59131811,  0.3766525 ],
        [ 0.39825407,  0.36450161,  0.57376683, ...,  0.19573934,
          0.37800259,  0.27944543],
        [ 0.33344936,  0.28214562,  0.35640102, ...,  0.2227413 ,
          0.23759238,  0.20113973]]], shape=(144, 73, 31))
Axis:
 * x: <pyinterp.core.Axis>
       min_value: 0
       max_value: 357.5
       step     : 2.5
       is_circle: true
 * y: <pyinterp.core.Axis>
       min_value: -90
       max_value: 90
       step     : -2.5
       is_circle: false
 * z: <pyinterp.core.TemporalAxis>
       min_value: 2002-07-01T12:00:00.000000000
       max_value: 2002-07-31T12:00:00.000000000
       step     : 86400000000000 nanoseconds