Note
Click here to download the full example code or to run this example in your browser via Binder
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.
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)
