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 pandas
import pyinterp
import pyinterp.backends.xarray
import pyinterp.tests
import numpy
import xarray

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:

Axis([90, 87.5, 85, ..., -87.5, -90], 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:

Axis([0, 2.5, 5, ..., 355, 357.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:

TemporalAxis(array(['2002-07-01T12:00:00.000000', '2002-07-02T12:00:00.000000',
                    '2002-07-03T12:00:00.000000', '2002-07-04T12:00:00.000000',
                    '2002-07-05T12:00:00.000000', '2002-07-06T12:00:00.000000',
                    '2002-07-07T12:00:00.000000', '2002-07-08T12:00:00.000000',
                    '2002-07-09T12:00:00.000000', '2002-07-10T12:00:00.000000',
                    '2002-07-11T12:00:00.000000', '2002-07-12T12:00:00.000000',
                    '2002-07-13T12:00:00.000000', '2002-07-14T12:00:00.000000',
                    '2002-07-15T12:00:00.000000', '2002-07-16T12:00:00.000000',
                    '2002-07-17T12:00:00.000000', '2002-07-18T12:00:00.000000',
                    '2002-07-19T12:00:00.000000', '2002-07-20T12:00:00.000000',
                    '2002-07-21T12:00:00.000000', '2002-07-22T12:00:00.000000',
                    '2002-07-23T12:00:00.000000', '2002-07-24T12:00:00.000000',
                    '2002-07-25T12:00:00.000000', '2002-07-26T12:00:00.000000',
                    '2002-07-27T12:00:00.000000', '2002-07-28T12:00:00.000000',
                    '2002-07-29T12:00:00.000000', '2002-07-30T12:00:00.000000',
                    '2002-07-31T12:00:00.000000'], dtype='datetime64[us]'))

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: 9.641530
pyinterp.Axis 1.183933

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: 73.936987
pyinterp.Axis 3.178607

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))

tcw = tcw.T

• 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>
Axis:
  x: Axis([0, 2.5, 5, ..., 355, 357.5], is_circle=true)
  y: Axis([90, 87.5, 85, ..., -87.5, -90], is_circle=false)
  z: TemporalAxis(array(['2002-07-01T12:00:00.000000', '2002-07-02T12:00:00.000000',
                    '2002-07-03T12:00:00.000000', '2002-07-04T12:00:00.000000',
                    '2002-07-05T12:00:00.000000', '2002-07-06T12:00:00.000000',
                    '2002-07-07T12:00:00.000000', '2002-07-08T12:00:00.000000',
                    '2002-07-09T12:00:00.000000', '2002-07-10T12:00:00.000000',
                    '2002-07-11T12:00:00.000000', '2002-07-12T12:00:00.000000',
                    '2002-07-13T12:00:00.000000', '2002-07-14T12:00:00.000000',
                    '2002-07-15T12:00:00.000000', '2002-07-16T12:00:00.000000',
                    '2002-07-17T12:00:00.000000', '2002-07-18T12:00:00.000000',
                    '2002-07-19T12:00:00.000000', '2002-07-20T12:00:00.000000',
                    '2002-07-21T12:00:00.000000', '2002-07-22T12:00:00.000000',
                    '2002-07-23T12:00:00.000000', '2002-07-24T12:00:00.000000',
                    '2002-07-25T12:00:00.000000', '2002-07-26T12:00:00.000000',
                    '2002-07-27T12:00:00.000000', '2002-07-28T12:00:00.000000',
                    '2002-07-29T12:00:00.000000', '2002-07-30T12:00:00.000000',
                    '2002-07-31T12:00:00.000000'], dtype='datetime64[us]'))
Data:
  [[[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]]]

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>
Axis:
  x: Axis([0, 2.5, 5, ..., 355, 357.5], is_circle=true)
  y: Axis([-90, -87.5, -85, ..., 87.5, 90], is_circle=false)
  z: TemporalAxis(array(['2002-07-01T12:00:00.000000000', '2002-07-02T12:00:00.000000000',
                    '2002-07-03T12:00:00.000000000', '2002-07-04T12:00:00.000000000',
                    '2002-07-05T12:00:00.000000000', '2002-07-06T12:00:00.000000000',
                    '2002-07-07T12:00:00.000000000', '2002-07-08T12:00:00.000000000',
                    '2002-07-09T12:00:00.000000000', '2002-07-10T12:00:00.000000000',
                    '2002-07-11T12:00:00.000000000', '2002-07-12T12:00:00.000000000',
                    '2002-07-13T12:00:00.000000000', '2002-07-14T12:00:00.000000000',
                    '2002-07-15T12:00:00.000000000', '2002-07-16T12:00:00.000000000',
                    '2002-07-17T12:00:00.000000000', '2002-07-18T12:00:00.000000000',
                    '2002-07-19T12:00:00.000000000', '2002-07-20T12:00:00.000000000',
                    '2002-07-21T12:00:00.000000000', '2002-07-22T12:00:00.000000000',
                    '2002-07-23T12:00:00.000000000', '2002-07-24T12:00:00.000000000',
                    '2002-07-25T12:00:00.000000000', '2002-07-26T12:00:00.000000000',
                    '2002-07-27T12:00:00.000000000', '2002-07-28T12:00:00.000000000',
                    '2002-07-29T12:00:00.000000000', '2002-07-30T12:00:00.000000000',
                    '2002-07-31T12:00:00.000000000'], dtype='datetime64[ns]'))
Data:
  [[[ 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  ]]]