如何用Shapefile裁剪非标准网格坐标的xarray Dataset

我有一个NetCDF数据集，维度是米为单位的xc和yc，CRS为Lambert方位等面积投影，参数如下：

Grid Mapping Name: Lambert Azimuthal Equal Area
Longitude of Projection Origin: 0.0
Latitude of Projection Origin: -90.0
False Easting: 0.0
False Northing: 0.0
Semi-major Axis: 6378137.0
Inverse Flattening: 298.257223563
Proj4 String: +proj=laea +lon_0=0 +datum=WGS84 +ellps=WGS84 +lat_0=-90.0

尝试用Shapefile裁剪时遇到问题，现有代码及问题如下：

读取NetCDF数据集 #

import geopandas as gpd
import rioxarray
import xarray as xr

# Open the netCDF dataset
ds = xr.open_dataset("./sea ice_2021_07.nc")
print(ds)

数据集输出：

<xarray.Dataset> Dimensions: (xc: 432, yc: 432) Coordinates: * xc
(xc) float64 -5.388e+03 -5.362e+03 ... 5.362e+03 5.388e+03 * yc
(yc) float64 5.388e+03 5.362e+03 ... -5.362e+03 -5.388e+03
time datetime64[ns] ...
lat (yc, xc) float32 ...
lon (yc, xc) float32 ... Data variables:
ice_conc (yc, xc) float64 ...

读取Shapefile #

gdf = gpd.read_file('polygon.shp').geometry.to_list()
print(gdf)

输出（经纬度格式）：

[<POLYGON ((-180 -64.4, -180 -64.4, -180 -64.5, -180 -64.5, -180 -64.6, -180 ...>]

裁剪代码及错误结果 #

ds.rio.write_crs("epsg:4326", inplace=True)
clipped = ds.rio.clip(gdf, crs="epsg:4326")
print(clipped)

裁剪结果仅显示1个纬度，且所有值为null：

<xarray.Dataset>
Dimensions: (xc: 13, yc: 1)
Coordinates:

• xc (xc) float64 -137.5 -112.5 -87.5 ... 137.5 162.5
• yc (yc) float64 -62.5
  time datetime64[ns] ...
  lat (yc, xc) float32 ...
  lon (yc, xc) float32 ...
  Lambert_Azimuthal_Grid int64 0
  Data variables:
  ice_conc (yc, xc) float64 nan nan nan nan ... nan nan nan nan

问题根源是错误地将数据集CRS设为WGS84（epsg:4326），但实际数据集的xc/yc是Lambert方位等面积投影下的米制坐标，需正确匹配CRS后再进行投影转换。

步骤1：为NetCDF数据集设置正确的CRS #

数据集的xc/yc是Lambert方位等面积投影（南极点为原点），可用给定的Proj4字符串或匹配的EPSG编码（EPSG:3031）设置CRS：

# 方式1：用Proj4字符串设置
ds.rio.write_crs("+proj=laea +lon_0=0 +datum=WGS84 +ellps=WGS84 +lat_0=-90.0", inplace=True)

# 方式2：用EPSG编码（更简洁，参数匹配）
# ds.rio.write_crs("epsg:3031", inplace=True)

步骤2：将Shapefile转换为数据集的CRS #

Shapefile是WGS84经纬度格式，需转换为与数据集一致的Lambert投影，确保裁剪时坐标系统匹配：

# 读取Shapefile并保留GeoDataFrame格式
gdf = gpd.read_file('polygon.shp')
# 转换到数据集的CRS
gdf_proj = gdf.to_crs(ds.rio.crs)

步骤3：执行裁剪操作 #

用转换后的Shapefile进行裁剪：

# 裁剪数据集
clipped_ds = ds.rio.clip(gdf_proj.geometry, crs=ds.rio.crs)
# 查看结果
print(clipped_ds)

关键说明 #

• 不能直接将数据集CRS设为epsg:4326，因为xc/yc是投影坐标而非经纬度，强制设置会导致坐标匹配完全错误。
• 裁剪前必须确保Shapefile和数据集的CRS一致，推荐将Shapefile转成数据集的投影（避免投影坐标转经纬度的精度损失）。
• 若数据集中已有lat/lon变量，也可指定经纬度作为空间维度，但效率和精度不如原生投影坐标。

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