如果点不包含在多边形中，则覆盖 pcolormesh 中的点 [英] Overwrite points from pcolormesh if they aren't contained in a polygon

问题描述

我正在尝试绘制地图，其中使用 pcolormesh 或 contourf 在土地上绘制空间模式.然后将描述英国边界的 shapefile/多边形叠加到该图上.我的问题是如何直接访问落在多边形之外的点以将它们设置为 0 或直接将它们着色为单一颜色，例如白色的.请参阅以下最小工作示例

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')

我尝试排除原始数据集中的任何点，然后生成 pcolormesh.但是， pcolormesh 在点之间进行插值.这导致从北爱尔兰到康沃尔产生一系列点.为了清楚起见，我想要的是填充多边形之外.感谢您的帮助.

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.

推荐答案

不是你要求的(修改 z 的值)，我绘制另一层 pcolormesh()在顶部以获得所需的效果.在此过程中，使用从 point-within_polygon 操作获得的各个值创建了一个 new_z 数组.创建自定义颜色图 new_binary 以与 new_z 一起使用以绘制该图层并获得最终图.

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()

这篇关于如果点不包含在多边形中，则覆盖 pcolormesh 中的点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！