**Task:** Interpolate data from regular to curvilinear grid

**Solution:** pyresample

Following two excellent contributions on interpolation between grids by Nikolay Koldunov and Oleksandr Huziy I would like to introduce a solution using the pyresample package. I feel it is timely since pyresample does encapsulate the strategy presented by Oleksandr (which I totally support) in fewer function calls. There might also be a speed-up factor to consider for big datasets, since pyresample comes with its own implementation of KD-Trees which was tested faster than the scipy.spatial.cKDTree.

The same data as in Nikolay's and Oleksandr's post will be used for easing comparison.

Some necessary imports:

```
%pylab inline
```

```
#for netcdf
from netCDF4 import Dataset
#for plotting
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
#for array manipulation
import numpy as np
#for interpolation (you will have to install pyresample first)
import pyresample
#for downloading files
import urllib2
```

I repeat Oleksandr's routine for downloading data, and download of the NCEP dataset, as well as the target curvilinear grid:

```
def download_from_link(path):
#download if it does not exist yet
import os
f_name = os.path.basename(path)
if not os.path.isfile(f_name):
remote_con = urllib2.urlopen(path)
with open(f_name, "wb") as f:
f.write(remote_con.read())
```

Load NCEP reanalysis data.

```
download_from_link("ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc")
```

We also need coordinates of the curvilinear grid:

```
download_from_link("https://dl.dropboxusercontent.com/u/107639/curv_grid.nc")
```

Now I'll skip exploring the data (take a look at Oleksandr's post) and concentrate on the reprojection of the NCEP values and the model's curvilinear grid.

```
# load lat-lon of the target grid
fc = Dataset('curv_grid.nc')
lon_curv = fc.variables['xc'][0,:,:]
lat_curv = fc.variables['yc'][0,:,:]
fc.close()
# load lat-lon-value of the origin data
fr = Dataset('air.mon.mean.nc')
air = fr.variables['air'][0,:,:]
lat = fr.variables['lat'][:]
lon = fr.variables['lon'][:]
fr.close()
# get 2D versions of the lat and lon variables (note the -180 here!)
lon2d, lat2d = np.meshgrid(lon - 180, lat)
```

When all the necessary inputs are read in memory, that's where we start using pyresample's magic.

Create a pyresample object holding the origin (NCEP) grid:

```
orig_def = pyresample.geometry.SwathDefinition(lons=lon2d, lats=lat2d)
```

Create another pyresample object for the target (curvilinear) grid:

```
targ_def = pyresample.geometry.SwathDefinition(lons=lon_curv, lats=lat_curv)
```

Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...

```
air_nearest = pyresample.kd_tree.resample_nearest(orig_def, air, \
targ_def, radius_of_influence=500000, fill_value=None)
```

... then with the custom weight function as inverse of the distance (same as Oleksandr):

```
wf = lambda r: 1/r**2
```

```
air_idw = pyresample.kd_tree.resample_custom(orig_def, air, \
targ_def, radius_of_influence=500000, neighbours=10,\
weight_funcs=wf, fill_value=None)
```

(the warning indicates we might loose information by limiting to 10 neighbours...)

Pyresample also natively support gaussian-shaped weighting:

```
air_gauss = pyresample.kd_tree.resample_gauss(orig_def, air, \
targ_def, radius_of_influence=500000, neighbours=10,\
sigmas=250000, fill_value=None)
```

Now let us plot the results from these three approaches. We use imshow() that does not introduce additional interpolation.

```
fig = plt.figure(figsize=(5,15))
ax = fig.add_subplot(311)
ax.imshow(air_nearest,interpolation='nearest')
ax.set_title("Nearest neighbor")
ax = fig.add_subplot(312)
ax.imshow(air_idw,interpolation='nearest')
plt.title("IDW of square distance \n using 10 neighbors");
ax = fig.add_subplot(313)
ax.imshow(air_gauss,interpolation='nearest')
plt.title("Gauss-shape of distance (sigma=25km)\n using 10 neighbors");
```

Now let's say we do not want to resample to the curvilinear coordinates, but to a predefined Lambert Azimuthal Equal Area covering Southern Hemisphere (something based on NSIDC's EASE grids).

```
area_id = 'ease_nh'
name = 'Arctic EASE grid'
proj_id = 'ease_nh'
x_size = 425
y_size = 425
area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
proj_dict = {'a': '6371228.0', 'units': 'm', 'lon_0': '0','proj': 'laea', 'lat_0': '+90'}
targ_def = pyresample.geometry.AreaDefinition(area_id, name, proj_id, proj_dict, x_size, y_size, area_extent)
print targ_def
```

Same function call as above:

```
air_gauss = pyresample.kd_tree.resample_gauss(orig_def, air, \
targ_def, radius_of_influence=500000, neighbours=10,\
sigmas=250000, fill_value=None)
```

And visualization:

```
fig = plt.figure()
ax = fig.add_subplot(111)
im = ax.imshow(air_gauss.transpose(),interpolation='nearest')
fig.colorbar(im).set_label('NCEP air temp')
plt.title("Gauss-shape of distance (sigma=25km)\n using 10 neighbors");
```

Check the pyressample for more cool re-projection stuff, including resampling of satellite data!

## Comments !