Regridding regular netcdf data

Question

I have a netcdf file containing global sea-surface temperatures. Using matplotlib and Basemap, I've managed to make a map of this data, with the following code:

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

filename = '/Users/Nick/Desktop/SST/SST.nc'
fh = Dataset(filename, mode='r')

lons = fh.variables['LON'][:]
lats = fh.variables['LAT'][:]
sst = fh.variables['SST'][:].squeeze()

fig = plt.figure()

m = Basemap(projection='merc', llcrnrlon=80.,llcrnrlat=-25.,urcrnrlon=150.,urcrnrlat=25.,lon_0=115., lat_0=0., resolution='l')

lon, lat = np.meshgrid(lons, lats)
xi, yi = m(lon, lat)

cs = m.pcolormesh(xi,yi,sst, vmin=18, vmax=32)

m.drawmapboundary(fill_color='0.3')
m.fillcontinents(color='0.3', lake_color='0.3')
cbar = m.colorbar(cs, location='bottom', pad="10%", ticks=[18., 20., 22., 24., 26., 28., 30., 32.])
cbar.set_label('January SST (' + u'\u00b0' + 'C)')
plt.savefig('SST.png', dpi=300)

The problem is that the data is very high resolution (9km grid) which makes the resulting image quite noisy. I would like to put the data onto a lower resolution grid (e.g. 1 degree), but I'm struggling to work out how this could be done. I followed a worked solution to try and use the matplotlib griddata function by inserting the code below into my above example, but it resulted in 'ValueError: condition must be a 1-d array'.

xi, yi = np.meshgrid(lons, lats)

X = np.arange(min(x), max(x), 1)
Y = np.arange(min(y), max(y), 1)

Xi, Yi = np.meshgrid(X, Y)

Z = griddata(xi, yi, z, Xi, Yi)

I'm a relative beginner to Python and matplotlib, so I'm not sure what I'm doing wrong (or what a better approach might be). Any advice appreciated!

Solution

If you regrid your data to a coarser lat/lon grid using e.g. bilinear interpolation, this will result in a smoother field.

The NCAR ClimateData guide has a nice introduction to regridding (general, not Python-specific).

The most powerful implementation of regridding routines available for Python is, to my knowledge, the Earth System Modeling Framework (ESMF) Python interface (ESMPy). If this is a bit too involved for your application, you should look into

1. EarthPy tutorials on regridding (e.g. using Pyresample, cKDTree, or Basemap).
2. Turning your data into an Iris cube and using Iris' regridding functions.

Perhaps start by looking at the EarthPy regridding tutorial using Basemap, since you are using it already.

The way to do this in your example would be

from mpl_toolkits import basemap
from netCDF4 import Dataset

filename = '/Users/Nick/Desktop/SST/SST.nc'
with Dataset(filename, mode='r') as fh:
   lons = fh.variables['LON'][:]
   lats = fh.variables['LAT'][:]
   sst = fh.variables['SST'][:].squeeze()

lons_sub, lats_sub = np.meshgrid(lons[::4], lats[::4])

sst_coarse = basemap.interp(sst, lons, lats, lons_sub, lats_sub, order=1)

This performs bilinear interpolation (order=1) on your SST data onto a sub-sampled grid (every fourth point). Your plot will look more coarse-grained afterwards. If you do not like that, interpolate back onto the original grid with e.g.

sst_smooth = basemap.interp(sst_coarse, lons_sub[0,:], lats_sub[:,0], *np.meshgrid(lons, lats), order=1)