Overwrite points from pcolormesh if they aren't contained in a polygon

Question

I'm trying to plot a map whereby a spatial pattern is plotted over the land using pcolormesh or contourf. A shapefile/polygon that describes the border of the UK is then overlayed onto this plot. My problem is how to directly access the points that fall outside the polygon to set them as 0 or directly colour them a single colour e.g. white. See the following minimal working example

import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np

# Load polygon
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
UK = world[world.iso_a3 == "GBR"]
UK.boundary.plot()

# Simulate a spatial pattern
xlims = (-8, 3)
ylims = (49, 60)
resolution = 0.05
y, x = np.mgrid[slice(ylims[0], ylims[1] + resolution, resolution),
               slice(xlims[0], xlims[1] + resolution, resolution)]
z = 0.5*x+np.sin(x)**2+ np.cos(y)

# Plot
fig, ax=plt.subplots(figsize=(6, 10))
im = ax.pcolormesh(x, y, z, cmap='viridis')
fig.colorbar(im, ax=ax)
UK.boundary.plot(ax=ax, color='black')

I have tried excluding any points in the original dataset and then generating the pcolormesh. However, pcolormesh interpolates between points. This results in a series of points being generated from Northern Ireland down to Cornwall. Just to be clear, what I would like is to fill outside the polygon. Thanks for any help.

Answer 1

Rather than what you request (modifying the values of z), I plot another layer of pcolormesh() on top to get the desired effect. In this process, a new_z array is created with individual values obtained from point-within_polygon operation. A custom colormap, new_binary is created to use with new_z to plot this layer and get the final plot.

import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
from shapely.geometry import Point

# Load polygon
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
UK = world[world.iso_a3 == "GBR"]

# plot 1
#UK.boundary.plot()

# Simulate a spatial pattern
xlims = (-8, 3)
ylims = (49, 60)
resolution = 0.05  # 0.05

# slice()
y, x = np.mgrid[slice(ylims[0], ylims[1] + resolution, resolution),
               slice(xlims[0], xlims[1] + resolution, resolution)]
z = 0.5*x+np.sin(x)**2+ np.cos(y)

# Target geometry, for point inside/outside polygon operation
ukgeom = UK['geometry'].values[0]

def prep_z(Xs,Ys,Zs):
    # Xs,Ys: result of np.meshgrid(lon, lat)
    # Zs: function of(Xs,Ys) to be manipulated; here hard-coded as `new_z`
    for ro,(arow,acol) in enumerate(zip(Xs,Ys)):
        # need 2 level loop to handle each grid point
        for co,xiyi in enumerate(zip(arow,acol)):
            pnt1 = Point(xiyi)
            if pnt1.within(ukgeom):
                new_z[ro][co]=1  #0:white, 1:black with cm='binary'
            else:
                new_z[ro][co]=0
            pass
        pass
# Usage of the function above: prep_z(x,y,z)
# Result: new_z is modified.
# New z for locations in/outside-polygon operation
new_z = np.zeros(z.shape)
prep_z(x,y,z)

# create custom colormap to use later
new_binary = mpl.colors.ListedColormap(np.array([[1., 1., 1., 1.],
       [0., 0., 0., 0.]]))
# 0:white, 1:transparent with cm='new_binary'
# do not use alpha option to get the intended result

# Plot 2
fig, ax = plt.subplots(figsize=(6, 10))
im = ax.pcolormesh(x, y, z, cmap='viridis', zorder=1)

im2 = ax.pcolormesh(x, y, new_z, cmap=new_binary, zorder=2)  # do not use alpha to get transparent mask

UK.boundary.plot(ax=ax, color='black', zorder=10)
fig.colorbar(im, ax=ax, shrink=0.5)

plt.show()