How to plot POP2 grid (POP_gx1v7) data on map with xarray and cartopy?

Question

I am working with ocean variables of the LENS2 ensemble of CESM2 and am now trying to plot the data on a geographical grid. The ocean model POP2 is based on a displaced-pole grid with poles in Greenland and Antarctica. The netCDF file contains the related WGS84 latitudes and longitudes but I do not manage to plot the data, especially north of 60 degree North. I am conducting most of my work with xarray and would like to stay within this library, but if necessary I am also prepared to switch to other tools.

Ideally I would like to still be able to see the original grid cells instead of having the code interpolated. (Interpolation with plt.contour() works if I plot data below 60 degrees North. However, I need to plot data within the Nordic Seas and there the plot collapses.)

To show the problem, here is the simplest case:

field = ds_temp.TEMP.isel(time=0, z_t=0)
extent = (ds_temp.TLONG.min(), ds_temp.TLONG.max(), 
          ds_temp.TLAT.min(), ds_temp.TLAT.max())

fig = plt.figure(figsize=(12, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
cs = ax.imshow(field, extent=extent, 
               origin='lower', 
               transform=ccrs.PlateCarree())
ax.coastlines()
plt.show()

This results in the following plot:

As you can see the cartopy coastlines and the data’s land mask are shifted. How can I plot my data with cartopy and if possible staying with xarray functions?

Answer 1

Instead of calling ax.coastlines() you could try to draw your own coastlines around continents appearing as white areas on the map. Most likely, your orginal data in field has some invalid values for mainland areas of the map. If you print it, it could probably look something like this:

[[-1.7999999523162842 -1.7999999523162842 -1.7999999523162842 ...
  -1.7999999523162842 -1.7999999523162842 -1.7999999523162842]
 [-1.7999999523162842 -1.7999999523162842 -1.7999999523162842 ...
  -1.7999999523162842 -1.7999999523162842 -1.7999999523162842]
 [-1.7999999523162842 -1.7999999523162842 -1.7999999523162842 ...
  -1.7999999523162842 -1.7999999523162842 -1.7999999523162842]
 ...
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]]

Those double dashes look like missing data in masked arrays. So one way to plot coastlines could be using matplotlib.pyplot.contour with a slightly modified data, where those "missing values" are replaced with some other ones, allowing you to compare them against a suitable threshold for drawing contours:

import os
import matplotlib.pyplot as plt
from netCDF4 import Dataset as netcdf_dataset

from cartopy import config
import cartopy.crs as ccrs


# get the path of the file -i t can be found in the repo data directory.
fname = os.path.join(config["repo_data_dir"],
                     'netcdf', 'HadISST1_SST_update.nc'
                     )

# get the data, latitudes and longtitudes
dataset = netcdf_dataset(fname)
sst = dataset.variables['sst'][0, :, :]
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]

ax = plt.axes(projection=ccrs.PlateCarree())

# plot data
ax.contourf(lons,
             lats,
             sst,
             60,
             transform=ccrs.PlateCarree()
             )

# replace invalid values (corresponding to land) with some unrealistically high values
very_high_value = 100
fixed_sst = sst.filled(very_high_value)

# use a slightly smaller threshold value to plot contour lines around land
ax.contour(lons,
           lats,
           fixed_sst,
           levels=[very_high_value-1],
           colors='black',
           transform=ccrs.PlateCarree()
           )

plt.show()

The code above is based on this example, showing the data with matplotlib.pyplot.contourf, and resulting in the following picture:

It doesn't look particularly nice, but it seems to do what you want. You will need to choose a suitable projection. If you use imshow() for plotting the original data, then you may need to adjust the extent parameters in imshow() and contour() separately.