Python中不同坐标系shapefile与雷达NetCDF影像匹配及裁剪方法

前置依赖 #

先安装需要的第三方库：

pip install geopandas xarray rioxarray pyproj numpy

实现逻辑 #

核心思路是先统一两个数据的坐标系，再用矢量范围裁剪栅格，你可以选择将矢量转到雷达的投影坐标系下裁剪，也可以选择将雷达重投影到矢量的坐标系下再裁剪，以下是可运行的完整代码：

完整代码示例 #

import geopandas as gpd
import xarray as xr
import rioxarray
from pyproj import CRS

# ---------------------- 1. 定义两个数据的坐标系 ----------------------
# Shapefile的GDA_1994_Lambert_Conformal_Conic坐标系，按给出的参数构建PROJ字符串
crs_shape = CRS.from_proj4("+proj=lcc +lat_1=-36 +lat_2=-38 +lat_0=-37 +lon_0=145 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +units=m +no_defs")

# 雷达影像的Albers等面积投影坐标系，按给出的参数构建PROJ字符串
crs_radar = CRS.from_proj4("+proj=aea +lat_1=-36.3 +lat_2=-39.4 +lat_0=-37.852 +lon_0=144.752 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")

# ---------------------- 2. 读取并处理Shapefile ----------------------
# 读取Shapefile
shp = gpd.read_file("你的shapefile路径.shp")
# 给Shapefile设置坐标系
shp = shp.set_crs(crs_shape)
# 转换到雷达坐标系，方便后续裁剪
shp_radar_crs = shp.to_crs(crs_radar)

# ---------------------- 3. 读取并处理NetCDF雷达数据 ----------------------
# 读取NetCDF文件
ds = xr.open_dataset("你的雷达nc文件路径.nc")
# 雷达的x、y单位是km，转换为米（和投影单位匹配）
ds['x'] = ds['x'] * 1000
ds['y'] = ds['y'] * 1000
# 给数据集设置坐标系和空间维度
ds = ds.rio.set_crs(crs_radar)
ds = ds.rio.set_spatial_dims(x_dim="x", y_dim="y")

# ---------------------- 4. 裁剪雷达影像 ----------------------
# 用转换后的Shapefile裁剪雷达数据，all_touched设为True可以保留边界所有触及的栅格
ds_clipped = ds.rio.clip(shp_radar_crs.geometry.values, shp_radar_crs.crs, all_touched=True)

# ---------------------- 可选：将裁剪后的雷达转到Shapefile的坐标系 ----------------------
ds_clipped_reproj = ds_clipped.rio.reproject(crs_shape)

# ---------------------- 保存结果 ----------------------
# 保存为NetCDF
ds_clipped.to_netcdf("裁剪后的雷达数据.nc")
# 也可以保存为GTiff格式
ds_clipped.rio.to_raster("裁剪后的雷达数据.tif")

注意事项 #

• 如果你的Shapefile本身已经自带了正确的坐标系定义，可以不用手动构建crs_shape，直接用shp.crs获取即可
• 若NetCDF里的雷达数据有多个变量，裁剪操作会自动作用于所有带空间坐标的变量
• 可以根据你的需求调整all_touched参数，设为False时只有栅格中心落在矢量范围内才会保留

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