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()
Post a Comment for "Overwrite Points From Pcolormesh If They Aren't Contained In A Polygon"