Skip to content Skip to sidebar Skip to footer

Overwrite Points From Pcolormesh If They Aren't Contained In A Polygon

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

Solution 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]

defprep_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) inenumerate(zip(Xs,Ys)):
        # need 2 level loop to handle each grid pointfor co,xiyi inenumerate(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]=0passpass# 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()

enter image description here

Post a Comment for "Overwrite Points From Pcolormesh If They Aren't Contained In A Polygon"