Is it possible to plot contours on top of data in matplotlib (python)?

Question

Below is a figure I have created for a single set of data, with contours plotted over the top. This is exactly what I need, I just need it plotting as a series of subplots for my 7 model runs, with the colourful salinity data, and the contours plotted over the top.

However, when I try to do the subplots, I can get in the salinity data, but I cannot plot the contours on top (or the colorbar, but that's another issue). Ideally, this should look like 7 mini versions of figure 1. However, I get the following error:

TypeError: inner() got multiple values for argument 'ax'

which I don't get when I plot the same thing as a single figure, and can't find a way around.

Unfortunately, the datafiles are 5GB each, so I can't upload them. As an example, I have created a dummy dataset and a basic plotting script to try and plot the contours over the data. I have also created some contours, although they are not exactly like mine, as mine are bathymetry, so needed transforming.

Sorry for any mistakes - pretty new to python! Any advice is massively appreciated, and comments on how to improve my question-asking!

import xarray as xa
import numpy as np
import matplotlib.pyplot as plt
import cmocean.cm as cm
import cartopy.crs as ccrs

#-----------Example Salinity Data---------------

data1 = 15 + 8 * np.random.randn(4,3)
salinity1 = xa.DataArray(data1)
salinity1 = xa.DataArray(data1, dims=['lat', 'lon'])
lons = np.linspace(-2, 6, 3)
lats = np.linspace(54, 62, 4)
salinity1 = xa.DataArray(data1, coords=[lats, lons], dims=['lat', 'lon'])
salt_ocean1 = salinity1

#-----------Example Contour Data---------------

xp = np.arange(-8, 10, 2)
yp = np.arange(-8, 10, 2)
zp = np.ndarray((9,9))
for x in range(0, len(xp)):
    for y in range(0, len(yp)):
        zp[x][y] = xp[x]**2 + yp[y]**2


#-----------Plotting (single figure)---------------

map_proj = ccrs.PlateCarree(central_longitude = 6)


ax = salt_ocean1.plot(transform = ccrs.PlateCarree(),
                      subplot_kws={"projection":map_proj},
                      cmap = cm.haline,
                      add_colorbar=False)

ax2 = plt.contour(xp, yp, zp)

Answer 1

Your example contour data is way over off the coast of Africa, so I changed that a little to make it in the same area your example salinity data. (maybe you used the example "lons" twice, also for the "lats"?)

xp = np.arange(-3, 7, 2)
yp = np.arange(52, 62, 2)
xx, yy = np.meshgrid(xp, yp)
zp = (xx-xx.mean())**2 + (yy-yy.mean())**2

You can probably get your result with only using the plotting methods from xarray (which itself uses Matplotlib). See for example the "facetting" section in the docs:
http://xarray.pydata.org/en/stable/user-guide/plotting.html#faceting

I personally prefer creating the figure and axes using Matplotlib directly, and then passing the axes to plot onto to xarray, using ax=ax.

map_proj = ccrs.PlateCarree(central_longitude=6)

fig, axs = plt.subplots(3,3, figsize=(10,9.5), constrained_layout=True, facecolor='w',
                        subplot_kw=dict(projection=map_proj))
axs = axs.ravel()

# loop over each axes
for i, ax in enumerate(axs, start=1):
    
    # set map extent
    ax.set_extent([-5, 10, 50, 65], crs=ccrs.PlateCarree())
    
    # add some cartopy features
    ax.coastlines(resolution='50m', alpha=.4)
    
    # Plot model data
    qm = salt_ocean1.plot(ax=ax, cmap=cm.haline, add_colorbar=False, transform=ccrs.PlateCarree())
    c = ax.contour(xp, yp, zp, levels=8, cmap='Blues', linewidths=2, transform=ccrs.PlateCarree())
    
    # do this after using xarray (it overrides the title)
    ax.set_title(f'model run: {i}')

# create 1 colorbar for all subplots using (ax=axs)
cb = fig.colorbar(qm, ax=axs, label='Salinity [g/kg]', shrink=.5)

Note that the use of transform=ccrs.PlateCarree() assign that geographical projection to the data. For the Numpy arrays used with contour that's always necessary, since those don't have any metadata. I suspect that for the xarray dataset it depends on the available metadata, but with this example data it's needed. But this projection would need to change if your model or bathymetry data for example has a different projection.

The cartopy gallery has a lot of examples showing how to add any features like the coastlines, land, countries etc:
https://scitools.org.uk/cartopy/docs/v0.13/gallery.html