使用 Cartopy 和 Geopandas 遮罩形状文件外部的区域

Question

我有一组用 Cartopy 绘制的数据（经度、纬度、温度）。我可以给出的最小示例是下面的代码（只有 30 个数据点）

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.colors as clr
import pandas as pd
import numpy as np

from metpy.interpolate import interpolate_to_grid, remove_nan_observations
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

canada_east = -95
canada_west = -101.8
canada_north = 52.8
canada_south = 48.85

central_lon = (canada_east + canada_west)/2
central_lat = (canada_north + canada_south)/2

crs = ccrs.LambertConformal(central_longitude = central_lon, central_latitude = central_lat)
lat = np.array([49.8134 50.904  50.698  49.095  49.436  49.9607 49.9601 49.356  50.116
49.402  52.3472 50.411  49.24   49.876  49.591  49.905  49.498  49.088
49.118  50.5947 49.3776 49.148  49.1631 51.358  49.826  50.4324 49.96
49.68   49.875  50.829  51.572])
lon = np.array([-100.3721  -97.273   -99.068   -97.528  -100.308   -98.9054  -98.6367
-99.248   -96.434  -100.93   -101.1099 -100.893  -100.055   -99.909
-97.518   -99.354   -98.03    -99.325   -99.054   -98.0035 -100.5387
-100.491   -97.1454 -100.361   -96.776   -99.4392  -97.7463  -97.984
-95.92    -98.111  -100.488])
tem = np.array([-8.45   -4.026  -5.993  -3.68   -7.35   -7.421  -6.477  -8.03   -3.834
-13.04   -4.057  -8.79   -6.619 -10.89   -4.465  -8.41   -4.861  -9.93
-7.125  -4.424 -11.95   -9.56   -3.86   -7.17   -4.193  -7.653  -4.883
-5.631  -3.004  -4.738  -8.81])

xp, yp, _ = crs.transform_points(ccrs.PlateCarree(), lon, lat ).T
xp, yp, tem = remove_nan_observations(xp, yp, tem)

alt_x, alt_y, data = interpolate_to_grid( xp, yp, tem, minimum_neighbors=2, search_radius=240000, interp_type = 'barnes', hres = 1000)

# Create the figure and grid for subplots
fig = plt.figure(figsize=(17, 12))

# Main ax
ax = plt.subplot(111, projection=crs)
ax.set_extent([canada_west, canada_east, canada_south, canada_north], ccrs.PlateCarree())

# Ading province borders and country borders
provinces_bdr = cfeature.NaturalEarthFeature(category = 'cultural',
                                                name = 'admin_1_states_provinces_lines',
                                                scale = '50m',
                                                linewidth = 0.6,
                                                facecolor='none',
                                                )  # variable to add provinces border

country_bdr = cfeature.NaturalEarthFeature(category= 'cultural',
                                            name = 'admin_0_boundary_lines_land', 
                                            scale = '50m', 
                                            linewidth = 1.5,
                                            facecolor = 'none',
                                            edgecolor = 'k')
                                                
ax.add_feature(provinces_bdr, linestyle='--')
ax.add_feature(country_bdr, linestyle='--')

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.BORDERS)

cf = ax.pcolormesh(alt_x, alt_y, data, cmap=plt.cm.rainbow)

# Read the shape file and add it

shape_feature = ShapelyFeature(Reader('MB_AGregion_Perim_South.shp').geometries(), ccrs.epsg(26914), linewidth = 1, facecolor = (1, 1, 1, 0), edgecolor = (0.5, 0.5, 0.5, 1))

ax.add_feature(shape_feature)
plt.show()

结果如下：

其中内部的灰线是由形状文件生成的。现在我想将着色限制为仅在形状文件内部（因此灰线之外的区域不应由 pcolormesh 着色），但我找不到可行的方法。我已经阅读了这个例子和这个例子，但我无法理解它们。有没有一种简单的方法可以单独使用 geopandas 和/或 cartopy 来做到这一点？

抱歉，我无法在这里上传形状文件，这是我可以做的最好的最小示例。如果我应该做任何改进，请告诉我。我是堆栈溢出的新手，并且乐于接受批评。

Edit1：为了澄清，我希望颜色限制为的形状文件是我用 ShapelyFeature 读取的“MB_AGregion_Perim_South.shp”（我的代码的最后 4 行），它绘制了一条灰线，限制了我的大部分内容染色。

编辑2：正如@Michael Delgado所建议的，我添加了这行代码：

cat_gdf = geopandas.read_file('MB_AGregion_Perim_South.shp')
cat_gdf = cat_gdf.to_crs(epsg = 4326)
mask = shapely.vectorized.contains(cat_gdf.dissolve().geometry.item(), alt_x, alt_y)

其中 alt_x 和 alt_y 是插值结果（请查看上面的示例）。shape文件原来的epsg = 26914，所以我将其转换为4326。

问题是掩码包含所有错误值（这意味着它掩盖了所有内容）。我怀疑这是因为 alt_x 和 alt_y 是用 crs.transform_points(ccrs.PlateCarree(), lon, lat ).T 转换的坐标（如我的代码上面所示）。我四处搜索并尝试将形状文件转换为不同的 epsg 值，但它不起作用。另外，cat_gdf.geometry 是一个多重多边形。会不会是这里的原因呢？

对于将来遇到此问题的任何人，这里有解决方案的详细说明

Answer 1

快速 MRE：

import numpy as np, pandas as pd, geopandas as gpd
import matplotlib.pyplot as plt

x = np.arange(-126, -105, 0.1)
y = np.arange(25, 46, 0.1)
xx, yy = np.meshgrid(x, y)
xnorm = (xx - xx.min()) / (xx.max() - xx.min())
ynorm = (yy - yy.min()) / (yy.max() - yy.min())
v = np.cos((xnorm * 2 - 1) * np.pi) + np.sin((ynorm * 2 - 1) * np.pi)

gdf = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

fig, ax = plt.subplots()
ax.pcolormesh(xx, yy, v)
xlim, ylim = ax.get_xlim(), ax.get_ylim()
gdf.plot(ax=ax, color='none', edgecolor='k')
ax.set_xlim(*xlim)
ax.set_ylim(*ylim)

您可以使用 shapely.vectorized 使用 shapely.geometry 对象来屏蔽一组 x, y 点：

import shapely.vectorized
mask = shapely.vectorized.contains(gdf.dissolve().geometry.item(), xx, yy)

fig, ax = plt.subplots()
ax.pcolormesh(xx, yy, np.where(mask, v, np.nan))
xlim, ylim = ax.get_xlim(), ax.get_ylim()
gdf.plot(ax=ax, color='none', edgecolor='k')
ax.set_xlim(*xlim)
ax.set_ylim(*ylim)