基于Cartopy与Geopandas掩膜Shapefile范围外区域

• 现有包含经度（lon）、纬度（lat）、温度（temperature）字段的观测数据集，已通过Cartopy结合MetPy的interpolate_to_grid方法完成Barnes空间插值绘图，配套有仅含30个样点的最小可复现代码
• 绘图时已读取坐标系为EPSG:26914的MB_AGregion_Perim_South.shp文件，通过ShapelyFeature以灰色线条绘制了目标区域边界
• 核心需求：仅在Shapefile划定的边界范围内显示pcolormesh生成的插值色斑，边界外区域不显示填色
• 已尝试操作与异常：
  • 参考Basemap相关的set_clipped_path裁剪方案未能复现效果
  • 尝试通过geopandas读取Shapefile后转为EPSG:4326坐标系，调用shapely.vectorized.contains判断插值网格坐标alt_x、alt_y是否落在面内，得到的掩膜数组全为False，怀疑问题与插值网格为LambertConformal投影坐标、Shapefile为MultiPolygon类型有关
  • 诉求：仅通过geopandas和/或cartopy实现该掩膜效果的简便方案

之前掩膜计算全为False和MultiPolygon类型无关，核心原因是坐标参考系（CRS）不匹配：插值输出的alt_x/alt_y是LambertConformal投影下的平面坐标，而转换为EPSG:4326的Shapefile是经纬度地理坐标，二者数值参考基准、量级完全不同，做包含判断自然会返回全False结果。以下两种方案均可直接实现需求，无额外依赖。

方案1：统一坐标系手动生成掩膜（稳定性最高） #

核心逻辑是把Shapefile转换到和插值网格完全一致的投影坐标系下，再做包含判断生成掩膜，将区域外的插值结果设为NaN即可自动不渲染：

1. 读取Shapefile并转换CRS，和插值使用的LambertConformal投影保持完全一致

   import geopandas as gpd
   import numpy as np
   from shapely.ops import unary_union
   from shapely.vectorized import contains
   # 替换为你实际插值、绘图用的LambertConformal投影参数，必须完全一致
   from cartopy.crs import LambertConformal
   lcc_proj = LambertConformal(central_longitude=105, central_latitude=30)
   
   # 读取目标区域shp
   gdf = gpd.read_file("MB_AGregion_Perim_South.shp")
   # 关键：不要转EPSG:4326，转到和插值网格相同的LCC投影
   gdf_lcc = gdf.to_crs(lcc_proj)
   

2. 合并所有面要素为单个几何对象，简化后续判断逻辑

   region_geom = unary_union(gdf_lcc.geometry)
   

3. 生成掩膜并应用到插值温度网格

   # alt_x、alt_y为interpolate_to_grid输出的二维网格坐标，interp_temp为对应插值温度结果
   in_region_mask = contains(region_geom, alt_x, alt_y)
   # 区域外的值设为掩膜，pcolormesh会自动跳过掩膜区域不渲染
   temp_plot = np.ma.masked_where(~in_region_mask, interp_temp)
   

4. 正常调用pcolormesh绘制temp_plot即可，边界外不会出现填色。

方案2：调用Cartopy内置裁剪接口（代码量最少） #

不需要手动做坐标转换和掩膜计算，直接给绘制完成的色斑图对象设置裁剪路径即可：

import geopandas as gpd
from cartopy.feature import ShapelyFeature

# 读取shp，保留原始EPSG:26914坐标系即可
gdf = gpd.read_file("MB_AGregion_Perim_South.shp")
# 生成区域边界要素
region_feat = ShapelyFeature(gdf.geometry, crs=gdf.crs, facecolor="none", edgecolor="gray")

# 按原有逻辑绘制pcolormesh，拿到返回的网格对象
mesh = ax.pcolormesh(alt_x, alt_y, interp_temp, transform=lcc_proj, cmap="jet")
# 核心：按shp边界裁剪色斑图
mesh.set_clip_path(region_feat, ax.transData)
# 最后绘制区域边界线
ax.add_feature(region_feat, linewidth=1)

注意：如果出现裁剪偏移问题，检查pcolormesh传入的transform参数是否和生成alt_x/alt_y使用的投影完全一致，不要错传为PlateCarree()。

内容的提问来源于stack exchange，提问作者QHoang