带旋转网格的NetCDF与EPSG:3003 shp配准及裁剪问题求助

核心问题本质 #

旋转网格（rlon/rlat）属于旋转极坐标系，这类NetCDF不会自带常规经纬度或标准EPSG编码，直接给数据集赋值CRS会完全偏离真实坐标，这也是你之前尝试转经纬度、裁剪失败的核心原因。必须先精准识别旋转投影参数，再完成坐标匹配。

分步解决方案 #

1. 提取旋转网格的投影参数 #

首先查看NetCDF元数据，找到旋转极的关键参数——通常数据集里会有rotated_pole这类grid_mapping变量，其属性包含grid_north_pole_latitude（极纬度）、grid_north_pole_longitude（极经度）：

import xarray as xr
ds = xr.open_dataset('path/to/your/data.nc')
# 打印旋转极参数，替换为你数据集实际的grid_mapping变量名
print(ds.rotated_pole)

2. 定义旋转网格的自定义CRS #

用pyproj基于提取的参数创建对应CRS，注意椭球体参数（a/b）要和NetCDF保持一致：

from pyproj import CRS

# 替换为你从数据集里拿到的实际参数
pole_lat = ds.rotated_pole.grid_north_pole_latitude.values
pole_lon = ds.rotated_pole.grid_north_pole_longitude.values

rotated_crs = CRS.from_proj4(
    f"+proj=ob_tran +o_proj=longlat +lon_0={pole_lon} +o_lat_p={pole_lat} +a=6371229 +b=6371229"
)

3. 为NetCDF绑定正确的空间属性 #

用rioxarray明确指定空间维度，并写入旋转CRS：

import rioxarray

# 绑定rlon/rlat为x/y空间维度
ds = ds.rio.set_spatial_dims(x_dim="rlon", y_dim="rlat")
# 写入旋转CRS
ds = ds.rio.write_crs(rotated_crs)

4. 统一坐标系 #

为提升裁剪效率，将EPSG:3003的Shapefile转换为旋转网格的CRS：

import geopandas as gpd

gdf = gpd.read_file('path/to/your/shapefile.shp')
# 转换坐标系至旋转CRS
gdf_rotated = gdf.to_crs(rotated_crs)

5. 执行裁剪操作 #

现在可以正常用rioxarray完成裁剪：

clipped_ds = ds.rio.clip(gdf_rotated.geometry, from_disk=True)

可选：将结果转回常规坐标系 #

如果需要把裁剪后的数据转成EPSG:3003或经纬度：

# 转EPSG:3003
clipped_ds_3003 = clipped_ds.rio.reproject("EPSG:3003")
# 转WGS84经纬度
clipped_ds_wgs84 = clipped_ds.rio.reproject("EPSG:4326")

常见错误排查 #

• CRS识别失败：检查proj4参数是否正确，尤其是椭球体a/b值要和NetCDF元数据一致
• 裁剪后无数据：验证Shapefile转换后的范围是否与NetCDF的rlon/rlat范围重叠
• 空间维度不被识别：必须用rio.set_spatial_dims明确指定x/y轴为rlon/rlat

完整修正代码 #

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

# 加载数据
ds = xr.open_dataset('path/to/your/data.nc')
gdf = gpd.read_file('path/to/your/shapefile.shp')

# 提取旋转极参数
rotated_pole = ds.rotated_pole
pole_lat = rotated_pole.grid_north_pole_latitude.values
pole_lon = rotated_pole.grid_north_pole_longitude.values

# 定义旋转CRS
rotated_crs = CRS.from_proj4(
    f"+proj=ob_tran +o_proj=longlat +lon_0={pole_lon} +o_lat_p={pole_lat} +a=6371229 +b=6371229"
)

# 绑定NetCDF空间属性
ds = ds.rio.set_spatial_dims(x_dim="rlon", y_dim="rlat")
ds = ds.rio.write_crs(rotated_crs)

# 转换Shapefile坐标系
gdf_rotated = gdf.to_crs(rotated_crs)

# 裁剪数据
clipped_ds = ds.rio.clip(gdf_rotated.geometry, from_disk=True)

# 可选：转换为EPSG:3003
clipped_ds_3003 = clipped_ds.rio.reproject("EPSG:3003")

内容的提问来源于stack exchange，提问作者Hundesa Siraj Mohammed