Create interpolator objects#

In this example, we are going to build the basic objects allowing to carry out interpolations.

Before starting, we will examine the properties of a Cartesian grid and the different classes associated with these objects.

The first step is to open the NetCDF file and load the data. We use here the NetCDF4 library to detail the different steps, but we will see that we can automate the steps described below using the xarray objects library.

Step-by-step creation of grids#

import timeit

import netCDF4
import numpy
import pandas
import xarray

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

with netCDF4.Dataset(pyinterp.tests.grid3d_path()) as ds:
    lon, lat, time, time_units, tcw = ds.variables[
        "longitude"][:], ds.variables["latitude"][:], ds.variables[
            "time"][:], ds.variables["time"].units, ds.variables["tcw"][:]
    time = numpy.array(netCDF4.num2date(time, time_units),
                       dtype="datetime64[us]")

This regular 3-dimensional grid is associated with three axes:

  • longitudes,

  • latitudes and

  • time.

To perform the calculations quickly, we will build three objects that will be used by the interpolator to search for the data to be used. Let’s start with the y-axis representing the latitude axis.

y_axis = pyinterp.Axis(lat)
y_axis

Out:

<pyinterp.core.Axis>
  min_value: -90
  max_value: 90
  step     : -2.5
  is_circle: false

For example, you can search for the closest point to 0.12 degrees north latitude.

y_axis.find_index([0.12])

Out:

array([36])

Then, the x-axis representing the longitudinal axis. In this case, the axis is an axis representing a 360 degree circle.

x_axis = pyinterp.Axis(lon, is_circle=True)
x_axis

Out:

<pyinterp.core.Axis>
  min_value: 0
  max_value: 357.5
  step     : 2.5
  is_circle: true

The values -180 and 180 degrees represent the same point on the axis.

x_axis.find_index([-180]) == x_axis.find_index([180])

Out:

array([ True])

Finally, we create the time axis

t_axis = pyinterp.TemporalAxis(time)
t_axis

Out:

<pyinterp.core.TemporalAxis>
  min_value: 2002-07-01T12:00:00.000000
  max_value: 2002-07-31T12:00:00.000000
  step     : 86400000000 microseconds

As these objects must communicate in C++ memory space, we use objects specific to the library much faster than other data models and manage the axes representing a circle. For example if we compare these objects to Pandas indexes:

values = lon[10:20] + 1 / 3
index = pandas.Index(lon)
print("pandas.Index: %f" % timeit.timeit(
    "index.searchsorted(values)", globals=dict(index=index, values=values)))
print("pyinterp.Axis %f" % timeit.timeit(
    "x_axis.find_index(values)", globals=dict(x_axis=x_axis, values=values)))

Out:

pandas.Index: 6.921722
pyinterp.Axis 1.510240

This time axis is also very efficient compared to the pandas index.

index = pandas.Index(time)
values = time + numpy.timedelta64(1, "ns")
print("pandas.Index: %f" % timeit.timeit(
    "index.searchsorted(values)", globals=dict(index=index, values=values)))
print("pyinterp.Axis %f" % timeit.timeit(
    "t_axis.find_index(values)", globals=dict(t_axis=t_axis, values=values)))

Out:

pandas.Index: 43.797473
<timeit-src>:6: UserWarning: implicit conversion turns datetime64[ns] into datetime64[us]
pyinterp.Axis 59.356697

Before constructing the tensor for pyinterp, we must begin to organize the tensor data so that it is properly stored in memory for pyinterp.

  • The shape of the tensor must be (len(x_axis), len(y_axis), len(t_axis))

  • The undefined values must be set to nan.

tcw[tcw.mask] = float("nan")

Now we can build the object handling the regular 3-dimensional grid.

Note

Grid data are not copied, the Grid3D class just keeps a reference on the handled array. Axis data are copied for non-uniform axes, and only examined for regular axes.

grid_3d = pyinterp.Grid3D(x_axis, y_axis, t_axis, tcw)
grid_3d

Out:

<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]]])
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.000000
       max_value: 2002-07-31T12:00:00.000000
       step     : 86400000000 microseconds

xarray backend#

The construction of these objects manipulating the regular grids can be done more easily using the xarray library and CF convention usually found in NetCDF files.

interpolator = pyinterp.backends.xarray.RegularGridInterpolator(
    xarray.open_dataset(pyinterp.tests.grid3d_path()).tcw)
interpolator.grid

Out:

<pyinterp.backends.xarray.Grid3D>
array([[[ 0.3334465 ,  0.28214264,  0.35639954, ...,  0.22274017,
          0.23759079,  0.20113754],
        [ 0.3942032 ,  0.35910034,  0.571064  , ...,  0.18763733,
          0.37395096,  0.27269363],
        [ 0.38745117,  0.44010544,  0.6574707 , ...,  0.19438934,
          0.5737648 ,  0.36045074],
        ...,
        [15.77317   , 10.124359  , 11.597317  , ...,  9.539768  ,
         11.003273  , 21.592094  ],
        [12.084702  ,  9.643726  , 12.967667  , ..., 10.471336  ,
         12.0685005 , 24.255836  ],
        [10.15271   , 10.667099  , 11.292194  , ..., 10.606346  ,
         15.883879  , 19.573696  ]],

       [[ 0.3334465 ,  0.28214264,  0.35639954, ...,  0.22274017,
          0.23759079,  0.20113754],
        [ 0.38880157,  0.35504913,  0.5656662 , ...,  0.18763733,
          0.36854935,  0.26729202],
        [ 0.38339996,  0.43065643,  0.6493721 , ...,  0.19573593,
          0.5602646 ,  0.34695053],
        ...,
        [15.729967  , 10.185116  , 11.417755  , ...,  9.635624  ,
         11.898388  , 21.96202   ],
        [12.195412  ,  9.657227  , 12.946064  , ..., 10.348476  ,
         12.546436  , 24.195084  ],
        [10.15271   , 10.667099  , 11.292194  , ..., 10.606346  ,
         15.883879  , 19.573696  ]],

       [[ 0.3334465 ,  0.28214264,  0.35639954, ...,  0.22274017,
          0.23759079,  0.20113754],
        [ 0.38475037,  0.34829712,  0.5602646 , ...,  0.18763733,
          0.36315155,  0.2632408 ],
        [ 0.37935257,  0.4198532 ,  0.6399193 , ...,  0.19978714,
          0.54406357,  0.3334465 ],
        ...,
        [15.689465  , 10.243168  , 11.23954   , ...,  9.73283   ,
         12.794853  , 22.331947  ],
        [12.306118  ,  9.670727  , 12.923113  , ..., 10.225616  ,
         13.027071  , 24.134329  ],
        [10.15271   , 10.667099  , 11.292194  , ..., 10.606346  ,
         15.883879  , 19.573696  ]],

       ...,

       [[ 0.3334465 ,  0.28214264,  0.35639954, ...,  0.22274017,
          0.23759079,  0.20113754],
        [ 0.406353  ,  0.37395096,  0.5737648 , ...,  0.21328735,
          0.38880157,  0.28889465],
        [ 0.4050026 ,  0.47250748,  0.7141762 , ...,  0.21598816,
          0.62641907,  0.4077034 ],
        ...,
        [15.458597  , 10.302574  , 11.629719  , ...,  9.365604  ,
          9.727432  , 19.622301  ],
        [11.724224  ,  9.73148   , 12.9447155 , ..., 10.887165  ,
         11.015423  , 24.116777  ],
        [10.15271   , 10.667099  , 11.292194  , ..., 10.606346  ,
         15.883879  , 19.573696  ]],

       [[ 0.3334465 ,  0.28214264,  0.35639954, ...,  0.22274017,
          0.23759079,  0.20113754],
        [ 0.4023018 ,  0.36854935,  0.5737648 , ...,  0.20383835,
          0.38475037,  0.28484344],
        [ 0.39825058,  0.46170807,  0.69527435, ...,  0.20788956,
          0.60886765,  0.39285278],
        ...,
        [15.563906  , 10.243168  , 11.618919  , ...,  9.425007  ,
         10.15406   , 20.278448  ],
        [11.844383  ,  9.701778  , 12.952816  , ..., 10.748104  ,
         11.366451  , 24.16403   ],
        [10.15271   , 10.667099  , 11.292194  , ..., 10.606346  ,
         15.883879  , 19.573696  ]],

       [[ 0.3334465 ,  0.28214264,  0.35639954, ...,  0.22274017,
          0.23759079,  0.20113754],
        [ 0.39825058,  0.36449814,  0.5737648 , ...,  0.19573593,
          0.37800217,  0.27944183],
        [ 0.39285278,  0.45225525,  0.6750221 , ...,  0.20248795,
          0.5913162 ,  0.37665176],
        ...,
        [15.667862  , 10.182415  , 11.606766  , ...,  9.483063  ,
         10.5779915 , 20.934595  ],
        [11.964542  ,  9.673428  , 12.9609165 , ..., 10.610394  ,
         11.7161255 , 24.208584  ],
        [10.15271   , 10.667099  , 11.292194  , ..., 10.606346  ,
         15.883879  , 19.573696  ]]], dtype=float32)
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

Total running time of the script: ( 1 minutes 51.745 seconds)

Gallery generated by Sphinx-Gallery