How to mask data that appears in the ocean using cartopy and matplotlib

Question

Not at all sure what I'm doing wrong besides perhaps the order that I am plotting the ocean in. I am trying to get the ocean feature in to mask the data in the ocean. I am trying to get data to not appear in the ocean and to get the ax.add_feature(cfeature.OCEAN) to be on top of the temperature data I am plotting so I see ocean and no data. Similar to what is happening in the great lakes region where you see lakes and no temperature data.

proj_map = ccrs.Mercator(central_longitude=cLon)
proj_data = ccrs.PlateCarree()

fig = plt.figure(figsize=(30,20))
ax = fig.add_subplot(1,1,1, projection=proj_map)
ax.set_extent([-84,-66,37,47.5])

CT = ax.contourf(Tlat, Tlon, tempF, transform=temp.metpy.cartopy_crs, levels=clevs, 
cmap=cmap)

ax.add_feature(cfeature.COASTLINE.with_scale('10m'), linewidth=0.5)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.BORDERS, linewidth=0.5)
ax.add_feature(cfeature.STATES.with_scale('10m'), linewidth=0.5)
ax.add_feature(USCOUNTIES.with_scale('20m'), linewidth=0.25)

cbar = fig.colorbar(CT, orientation='horizontal', shrink=0.5, pad=0.05)
cbar.ax.tick_params(labelsize=14)
cbar.set_ticks([-50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 
110, 120])
cbar.ax.set_xlabel("Temp ($^\circ$F)",fontsize=20)

Here is what the image looks like

Answer 1

You need to use zorder option to specify proper orders of the plot on the map. Features with largers values of zorder will be plotted on top of those with lower values. In your case, you need zorder of the OCEAN larger than the filled-contour.

Here is a runnable demo code and its sample plot. Read comments in the code for explanation.

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

fig, ax = plt.subplots(figsize=(8, 8), subplot_kw=dict(projection=ccrs.PlateCarree()))
extent = [-84, -66, 37, 47.5]

# generate (x, y), centered at the middle of the `extent`
mean = [(extent[0]+extent[1])/2, (extent[2]+extent[3])/2] #mean
cov = [[7, 3.5], [3.5, 6]]  #co-variance matrix
x, y = np.random.multivariate_normal(mean, cov, 4000).T

# make a 2D histogram
# set the edges of the bins in x and y directions
bin_size = 40
lonrange = np.linspace(extent[0], extent[1], bin_size)
latrange = np.linspace(extent[2], extent[3], bin_size)

# the cell sizes of the bins:
dx = (lonrange[1]- lonrange[0])/2
dy = (latrange[3]- latrange[2])/2

# compute array of center points of the bins' grid
# the dimensions of mesh-grid < the edges by 1
lonrange2 = np.linspace(extent[0]+dx, extent[1]-dx, bin_size-1)
latrange2 = np.linspace(extent[2]+dy, extent[3]-dy, bin_size-1)
x2d, y2d = np.meshgrid(lonrange2, latrange2)

# create 2d-histogram
# zorder is set = 10
h = ax.hist2d(x, y, bins=[lonrange, latrange], zorder=10, alpha=0.75)
#h: (counts, xedges, yedges, image)

ax.add_feature(cfeature.OCEAN, zorder=12)  #zorder > 10
ax.add_feature(cfeature.BORDERS, linewidth=0.5)

ax.gridlines(draw_labels=True, xlocs=list(range(-85, -60, 5)), ylocs=list(range(35, 50, 5)), 
             linewidth=1.8, color='gray', linestyle='--', alpha=0.8, zorder=20)


# plot colorbar, using image from hist2d's result
plt.colorbar(h[3], ax=ax, shrink=0.45)
# finally, show the plot.
plt.show()

The output plot:

If zorder option is not specified:

ax.add_feature(cfeature.OCEAN)

the plot will be: