Working with NEMO files

Here are some examples for plotting outputs on a NEMO grid.

Building variable

Let’s first build a bathymetric field using the vertical scale factors.

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.feature as cfeature
import as ccrs

data = xr.open_dataset('data/')
data = data.isel(t=0)
tmask = data['tmask'].values
e3t = data['e3t_0'].values
lon = data['glamt'].values
lat = data['gphit'].values

bathy = np.sum(e3t * tmask, axis=0)
bathy = == 0, bathy)

First try

If we first try to use the pcolormesh as we learned, here is what comes out:

fig = plt.figure(figsize=(12, 15))
ax = plt.axes(projection=ccrs.PlateCarree())
cs = ax.pcolormesh(lon, lat, bathy, transform=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND, zorder=50)
ax.add_feature(cfeature.COASTLINE, zorder=51)
cb = plt.colorbar(cs, shrink=0.3)
/home/barrier/Softwares/anaconda3/lib/python3.7/site-packages/cartopy/mpl/ MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  result = matplotlib.axes.Axes.pcolormesh(self, *args, **kwargs)

We have an error message saying that the longitudes and latitudes are not monotonic. Let’s improve our figure.

Better way

As indicated in the documentation, the pcolormesh method assumes the following layout:

(X[i+1, j], Y[i+1, j])          (X[i+1, j+1], Y[i+1, j+1])
                      | C[i,j] |
    (X[i, j], Y[i, j])          (X[i, j+1], Y[i, j+1])

Therefore, the good way to draw is to provide the coordinates of the F points (upper-right corners), and to give a sub-array of T points.

Grid layout of NEMO outputs:

lonf = data['glamf'].data
latf = data['gphif'].data
fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=0))
cs = ax.pcolormesh(lonf, latf, bathy[1:, 1:], transform=ccrs.PlateCarree(), shading='flat')
ax.set_extent([-20, 20, -20, 20], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE, zorder=2)
ax.add_feature(cfeature.LAND, zorder=1)
cb = plt.colorbar(cs, shrink=0.5)

Contour plots

However, the drawing of contour plots is not simple on irregular grid. Instead, we need to use the tricontour method, as indicated here.

First, we recover the coordinates on the T points, not on the F points as for pcolormesh.

lont = data['glamt'].data
latt = data['gphit'].data

Then, we extract the data mask.

mask = (

Now, we extract the bathy, longitudes and latitudes on wet points and we convert into 1D arrays:

lon1d = np.ravel(lont[~mask])
lat1d = np.ravel(latt[~mask])
bat1d = np.ravel(bathy[~mask])
masked_array(data=[121.15876481013115, 64.95892756174798,
                   515.0322590608843, ..., 24.57716311091656,
                   24.57716311091656, 24.57716311091656],
             mask=[False, False, False, ..., False, False, False],

The next step is to convert our 1D geographical coordinates (lon/lat) into the coordinates of the output map. If we want to draw our contours on a Mollweide projection:

projin = ccrs.PlateCarree()
projout = ccrs.Mollweide(central_longitude=180)
#projout = ccrs.PlateCarree(central_longitude=0)

output = projout.transform_points(projin, lon1d, lat1d)
lonout = output[..., 0]
latout = output[..., 1]

Now, we can add contours using the tricontour method:

fig = plt.figure(figsize=(12, 12))
ax = plt.axes(projection=projout)
cs = ax.pcolormesh(lonf, latf, bathy[1:, 1:], transform=projin)
cl = ax.tricontour(lonout, latout, bat1d, levels=np.arange(0, 6000 + 1000, 1000), colors='k', linewidths=0.5)
ax.add_feature(cfeature.LAND, zorder=100)
l = ax.add_feature(cfeature.COASTLINE, zorder=101, linewidth=2)