Task:
Interpolate data from regular to curvilinear grid
Solution:
Basemap.interp functionUnfortunately geophysical data distributed on a large variety of grids, and from time to time we have to compare our variables to each other. Often plotting a simple map is enough, but if you want to go a bit beyond qualitative comparison then you have to interpolate data from one grid to another. One of the easiest way to do this is to use basemap.interp function from Matplotlib Basemap library. Here I will show how to prepare your data and how to perform interpolation.
Some necessary imports:
%pylab inline
from scipy.io import netcdf
As usual we are going to use NCEP reanalysis data. This time monthly mean air temperature. Load with wget:
!wget ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc
We also need coordinates of the curvilinear grid:
!wget https://www.dropbox.com/s/9xzgyjs08zyuwzw/curv_grid.nc
Let's have a look at our curvilinear grid:
fc = netcdf.netcdf_file('curv_grid.nc')
lon_curv = fc.variables['xc'][0,:,:]
lat_curv = fc.variables['yc'][0,:,:]
Distribution of longitudes, for example, looks like this:
imshow(lon_curv)
This is grid from the regional Arctic Ocean model with 40 km resolution:
Now load NCEP data (only first time step), that have 2.5 degree resolution:
fr = netcdf.netcdf_file('air.mon.mean.nc')
air = fr.variables['air'][0,:,:]
lat = fr.variables['lat'][:]
lon = fr.variables['lon'][:]
First look at the data:
imshow(air)
#colorbar()
and longitudes (1d):
lon
There are couple of problems with this longitude. First - Basemap wants it to be in the -180:180 format, not 0:360. This can be fixed:
lon1 = lon.copy()
for n, l in enumerate(lon1):
if l >= 180:
lon1[n]=lon1[n]-360.
lon = lon1
lon
Second - the order of both lat and lon should be increasing. So we have to reshape our longitudes and data in a way that the first values on the left correspond to -180 longitude. We split arrays in the middle:
air1 = air[:,0:72]
air2 = air[:,72:]
lon1 = lon[0:72]
lon2 = lon[72:]
And then put this halves back together, but switching left and right:
air_new = np.hstack((air2, air1))
lon_new = np.hstack((lon2, lon1))
Now lons look OK:
lon_new
And this is our new data array:
imshow(air_new)
What about lats?
lat
They are not increasing, and we have to flip them up side down:
lat_new = np.flipud(lat)
lat_new
The data should be flipped accordingly:
air_new = np.flipud(air_new)
Now import mpl_toolkits.basemap:
import mpl_toolkits.basemap
and interpolate:
result = mpl_toolkits.basemap.interp(air_new, lon_new, lat_new, \
lon_curv, lat_curv, checkbounds=False, masked=False, order=1)
Here as an input we use our modified 1d coordinate variables and data, as well as two 2d arrays with coordinates of curvilinear grid we interpolate to. You can change type of interpolation by setting the order argument. It is 0 for nearest-neighbor interpolation, 1 for bilinear interpolation, 3 for cubic spline (default 1).
Quick look at the result:
imshow(result)
colorbar()
Not very clear what we get here. Let's put this data on the map:
from mpl_toolkits.basemap import Basemap
m = Basemap(projection='npstere',boundinglat=40,lon_0=0,resolution='l')
x, y = m(lon_curv, lat_curv)
lon_new2, lat_new2 = meshgrid(lon_new, lat_new)
x2, y2 = m(lon_new2, lat_new2)
fig = plt.figure(figsize=(10,10))
subplot(121)
m.drawcoastlines(linewidth=1.)
cs = m.contourf(x,y,result, levels = np.arange(-40,10,5), extend='both')
colorbar(orientation='horizontal')
subplot(122)
m.drawcoastlines(linewidth=1.)
cs = m.contourf(x2,y2,air_new, levels = np.arange(-40,10,5), extend='both')
colorbar(orientation='horizontal')
Look quite similar to me :)
But if we look at pcolor versions of figures, were grid is visible, you will see what basemap.interp did.
fig = plt.figure(figsize=(10,20))
subplot(211)
m.drawcoastlines(linewidth=1.)
cs = m.pcolormesh(x,y,result, vmin =-40, vmax=10)
subplot(212)
m.drawcoastlines(linewidth=1.)
cs = m.pcolormesh(x2,y2,air_new, vmin =-40, vmax=10)
Comments !