Interpolation example

Evaltools also has interpolation features. In interpolation module, you can find tools to perform bilinear interpolation at lat/lon locations from grid values. However, interpolation can be easily replaced by scipy or xarray.

First, we import a list of stations where we want interpolated values.

In [1]: import numpy as np

In [2]: import evaltools as evt

In [3]: listing_path = "../doc/sample_data/listing"

In [4]: stations = evt.utils.read_listing(listing_path)

Let’s define our domain boundaries and create arbitrarily simulated data.

In [5]: min_lat = 45.

In [6]: min_lon = 6.

In [7]: d_lon = .5

In [8]: d_lat = .5

In [9]: nb_lon = 28

In [10]: nb_lat = 14

In [11]: sim_a = np.array(
   ....:     [
   ....:         [lat*lon for lon in range(nb_lon)]
   ....:         for lat in range(nb_lat)
   ....:     ]
   ....: )
   ....: 

In [12]: sim_b = np.array(
   ....:     [
   ....:         [lat + lon for lon in range(nb_lon)]
   ....:         for lat in range(nb_lat)
   ....:     ]
   ....: )
   ....: 

We now have a grid of values for two species A and B.

Note

Such grids could have been loaded from a netcdf or a grib file using libraries like netCDF4 or xarray.

With xarray

In [13]: import xarray as xr

In [14]: interp_lat = xr.DataArray(stations.lat, dims='station')

In [15]: interp_lon = xr.DataArray(stations.lon, dims='station')

In [16]: data = xr.DataArray(
   ....:     np.flipud(sim_a),
   ....:     dims=('lat', 'lon'),
   ....:     coords={
   ....:         'lat': np.arange(min_lat, min_lat + d_lat*nb_lat, d_lat),
   ....:         'lon': np.arange(min_lon, min_lon + d_lon*nb_lon, d_lon),
   ....:     },
   ....: )
   ....: 

In [17]: data.interp(lat=interp_lat, lon=interp_lon)
Out[17]: 
<xarray.DataArray (station: 14)> Size: 112B
array([         nan, 134.42659828, 154.110888  , 140.68272232,
       139.47795   ,  17.70654577,  42.82960844,  43.67600354,
        41.72561363,  50.42877972,  51.38245376,  69.97642392,
        52.33722192,  80.67082356])
Coordinates:
    lat      (station) float64 112B 42.51 46.68 47.84 ... 49.57 49.73 49.86
    lon      (station) float64 112B 1.539 12.97 16.53 16.3 ... 15.08 13.4 18.27
  * station  (station) object 112B 'AD0942A' 'AT0VOR1' ... 'CZ0PPLA' 'CZ0TOPR'

With Scipy

In [18]: from scipy.interpolate import RegularGridInterpolator

In [19]: interpolator = RegularGridInterpolator(
   ....:     (
   ....:         np.arange(min_lat, min_lat + d_lat*nb_lat, d_lat),
   ....:         np.arange(min_lon, min_lon + d_lon*nb_lon, d_lon),
   ....:     ),
   ....:     np.flipud(sim_a),
   ....:     method='linear',
   ....:     bounds_error=False,
   ....:     fill_value=np.nan,
   ....: )
   ....: 

In [20]: interpolator(
   ....:     stations[['lat', 'lon']].values
   ....: )
   ....: 
Out[20]: 
array([         nan, 134.42659828, 154.110888  , 140.68272232,
       139.47795   ,  17.70654577,  42.82960844,  43.67600354,
        41.72561363,  50.42877972,  51.38245376,  69.97642392,
        52.33722192,  80.67082356])

With evaltools

We need to create a Grid object, load metadata from a station list file and load one list of stations per species where interpolation will be performed (in our case we will take the first five stations of the stations variable defined above for species A and the last five for species B).

In [21]: interpolator = evt.interpolation.Grid(
   ....:     min_lat, min_lon, d_lon, d_lat, nb_lon, nb_lat,
   ....: )
   ....: 

In [22]: interpolator.load_stations_metadata(listing_path)

In [23]: interpolator.load_station_lists(
   ....:     {'A': stations.index[:5], 'B': stations.index[-5:]}
   ....: )
   ....: 
*** 1 stations dropped since they do not appear in self.stations.index (these are stations not found in the metadata file or located outside the domain).

All that remains to be done is to load the input arrays and perform the interpolation.

In [24]: interpolator.set_grids({'A': sim_a, 'B': sim_b})

In [25]: interpolator.interpolate()
Out[25]: 
{'A': AT0VOR1    134.426627
 AT10001    154.110868
 AT31401    140.682664
 AT31402    139.477942
 dtype: float64,
 'B': CZ0ALIB    19.877197
 CZ0HHKB    22.302000
 CZ0JKOS    22.013803
 CZ0PPLA    18.339802
 CZ0TOPR    27.826797
 dtype: float64}

Note

Before performing the interpolation, you could use view method to visualize input array and check that values are not upside down.

Note

If you then want to perform interpolation from other input arrays for the same species you don’t have to load station lists again but only updating the arrays by calling set_grids method again.