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!

Corkwood

Have you try to use conturf, I'm finding it producing much smoother plots – Dews

Hi, yes I've tried contourf too, but it doesn't really solve the problem. With the high resolution data, the boundaries between contours are not smooth. It would be ok if I could smooth the data and then contour, but I didn't find a good way to do that yet. – Corkwood

How about passing exact values to contourf on which contours have to bee drawn, with really small intervals that should smooth data on map. – Dews

Answer 1

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)

Esdras

I found that scipy.interpolate.interp2d() is faster than basemap. – Forsberg

The new kid in town is xESMF - ESMF interface to game-changer xarray. – Esdras

xESMF requires Python3.5 or above. – Km

Answer 2

I usually run my data through a Laplace filter for smoothing. Perhaps you could try the function below and see if it helps with your data. The function can be called with or without a mask (e.g land/ocean mask for ocean data points). Hope this helps. T

# Laplace filter for 2D field with/without mask
# M = 1 on - cells used
# M = 0 off - grid cells not used
# Default is without masking

import numpy as np
def laplace_X(F,M):
    jmax, imax = F.shape
    # Add strips of land
    F2 = np.zeros((jmax, imax+2), dtype=F.dtype)
    F2[:, 1:-1] = F
    M2 = np.zeros((jmax, imax+2), dtype=M.dtype)
    M2[:, 1:-1] = M

    MS = M2[:, 2:] + M2[:, :-2]
    FS = F2[:, 2:]*M2[:, 2:] + F2[:, :-2]*M2[:, :-2]

    return np.where(M > 0.5, (1-0.25*MS)*F + 0.25*FS, F)

def laplace_Y(F,M):
    jmax, imax = F.shape

    # Add strips of land
    F2 = np.zeros((jmax+2, imax), dtype=F.dtype)
    F2[1:-1, :] = F
    M2 = np.zeros((jmax+2, imax), dtype=M.dtype)
    M2[1:-1, :] = M

    MS = M2[2:, :] + M2[:-2, :]
    FS = F2[2:, :]*M2[2:, :] + F2[:-2, :]*M2[:-2, :]

    return np.where(M > 0.5, (1-0.25*MS)*F + 0.25*FS, F)


# The mask may cause laplace_X and laplace_Y to not commute
# Take average of both directions

def laplace_filter(F, M=None):
    if M == None:
        M = np.ones_like(F)
    return 0.5*(laplace_X(laplace_Y(F, M), M) +
                laplace_Y(laplace_X(F, M), M))

Komatik

Answer 3

To answer your original question regarding scipy.interpolate.griddata, too:

Have a close look at the parameter specs for that function (e.g. in the SciPy documentation) and make sure that your input arrays have the right shapes. You might need to do something like

import numpy as np
points = np.vstack([a.flat for a in np.meshgrid(lons,lats)]).T # (n,D)
values = sst.ravel() # (n)

etc.

Esdras

Answer 4

If you are working on Linux, you can achieve this using nctoolkit (https://nctoolkit.readthedocs.io/en/latest/).

You have not stated the latlon extent of your data, so I will assume it is a global dataset. Regridding to 1 degree resolution would require the following:

import nctoolkit as nc
filename = '/Users/Nick/Desktop/SST/SST.nc'
data = nc.open_data(filename)

data.to_latlon(lon = [-179.5, 179.5], lat = [-89.5, 89.5], res = [1,1])
# visualize the data
data.plot()

Padnag

Answer 5

Look at this example with xarray... use the ds.interp method and specify the new latitude and longitude values.

http://xarray.pydata.org/en/stable/interpolation.html#example

Selfsacrifice