如何用Python将栅格降水数据裁剪至国家边界多边形？

完全可以用国家边界多边形精准裁剪目标区域的栅格数据，解决陆海掩码在高分辨率下丢失内陆数据的问题。以下是具体实现步骤和修改后的代码：

核心思路 #

• 加载目标区域的边界矢量数据（本地SHP文件或公开矢量数据集均可）
• 统一矢量数据与栅格数据的坐标系（你的数据使用ccrs.PlateCarree()，对应EPSG:4326）
• 用Xarray的clip方法对栅格数据进行裁剪，仅保留多边形范围内的网格点

修改后的完整代码 #

import geopandas as gpd
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf

# 1. 加载国家边界矢量数据
# 方式1：加载本地SHP文件
# country_boundary = gpd.read_file("path/to/your/country_shapefile.shp")
# 方式2：调用geopandas内置的低分辨率全球边界数据，筛选目标国家
country_boundary = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
country_boundary = country_boundary[country_boundary.name == 'South Africa']  # 匹配你的坐标区域（南非）

# 确保坐标系与栅格数据一致（EPSG:4326对应PlateCarree投影）
country_boundary = country_boundary.to_crs(epsg=4326)

# 2. 裁剪栅格数据
clipped_data = ds1am.rnd24[16][0].clip(country_boundary.geometry)

# 3. 绘制裁剪后的空间分布图
fig, axes = plt.subplots(nrows=1, ncols=1, subplot_kw={'projection': ccrs.PlateCarree()}, figsize=[10,6])

lev_prp=(0,100,200,300,400,500,600,700,800,900,1000,1100,1200)

im0 = clipped_data.plot.contourf(cmap=colourmap_prec, extend='max', add_colorbar=False, levels=lev_prp)

plt.colorbar(im0, fraction=0.035)

axes.set_xticks([26, 28], crs=ccrs.PlateCarree())
axes.set_yticks([-32, -34], crs=ccrs.PlateCarree())
axes.add_feature(cf.BORDERS, alpha=0.5)
axes.add_feature(cf.COASTLINE, alpha=0.5)

# 可选：叠加裁剪边界，验证裁剪范围
country_boundary.boundary.plot(ax=axes, edgecolor='darkred', linewidth=1.2, transform=ccrs.PlateCarree())

plt.tight_layout()
plt.show()

注意事项 #

• 若使用自定义SHP文件，需确认文件坐标系，不一致时用to_crs转换为EPSG:4326
• 如需更高精度的边界，建议使用1:10m分辨率的Natural Earth矢量数据替换默认低分辨率数据
• clip方法会自动保留多边形内部的所有栅格点，避免陆海掩码因分辨率不匹配导致的内陆数据丢失问题

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