如何在含rlat/rlon维度的NetCDF中使用lon/lat坐标？

我有一个NetCDF文件，对应的xarray.Dataset信息如下：

xarray.Dataset {
dimensions: 
    rlat = 183 ;
    rlon = 182 ;
    time = 1 ;

coordinates:
    float32 lat(rlat, rlon) ;
    float32 lon(rlat, rlon) ;
    datetime64[ns] time (time) ;
    float32 rlat(rlat) ;
    float32 rlon(rlon) ;

variables:
    float32 CaPA_coarse_A_PR_SFC(time, rlat, rlon) ;
    float32 rotated_pole() ;

attributes:
    product = CaPA_coarse ;
    Conventions = CF-1.6 ;
    License = These data are provided by the Canadian Surface Prediction Archive CaSPar. You should have received a copy of the License agreement with the data. Otherwise you can find them under http://caspar-data.ca/doc/caspar_license.txt or email caspar.data@uwaterloo.ca. ;
    Remarks = Variable names are following the convention <Product>_<Type:A=Analysis,P=Prediction>_<ECCC name>_<Level/Tile/Category>. Variables with level '10000' are at surface level. The height [m] of variables with level '0XXXX' needs to be inferrred using the corresponding fields of geopotential height (GZ_0XXXX-GZ_10000). The variables UUC, VVC, UVC, and WDC are not modelled but inferred from UU and VV for convenience of the users. Precipitation (PR) is reported as 6-hr accumulations for CaPA_fine and CaPA_coarse. Precipitation (PR) are accumulations since beginning of the forecast for GEPS, GDPS, REPS, RDPS, HRDPS, and CaLDAS.}

我希望用lon/lat坐标替代rlat/rlon维度，但不知道怎么操作；同时尝试将该文件重投影到EPSG:4326（WGS84），但不清楚应该用什么作为输入投影。

1. 确定输入投影 #

这个数据集采用的是旋转极投影，核心参数存储在rotated_pole变量的属性中。可以通过以下步骤获取并构建输入坐标系：

import xarray as xr
import pyproj

# 打开数据集
ds = xr.open_dataset("你的NetCDF文件路径.nc")
# 查看rotated_pole的属性，里面包含投影的关键参数
print(ds.rotated_pole.attrs)

# 根据属性构建投影坐标系（CaPA数据通常用以下参数，可根据实际输出调整）
proj_params = {
    "proj": "ob_tran",
    "o_proj": "latlon",
    "lon_0": ds.rotated_pole.attrs.get("grid_north_pole_longitude"),
    "o_lat_p": ds.rotated_pole.attrs.get("grid_north_pole_latitude"),
    "a": 6371229,
    "b": 6371229
}
input_crs = pyproj.CRS(proj_params)

2. 重投影到EPSG:4326并替换坐标维度 #

原数据的lon/lat是二维不规则坐标，无法直接替换一维的rlat/rlon维度，最实用的方式是重投影到规则的WGS84网格，推荐两种工具：

方式一：用xesmf做重采样 #

import xesmf as xe

# 创建目标规则网格（可自定义分辨率，这里以0.1°为例）
target_grid = xe.util.grid_2d(
    lon_min=ds.lon.min().item(), lon_max=ds.lon.max().item(),
    lat_min=ds.lat.min().item(), lat_max=ds.lat.max().item(),
    dx=0.1, dy=0.1
)

# 构建重投影器，选择插值方法（bilinear为双线性插值，可按需选nearest等）
regridder = xe.Regridder(ds, target_grid, "bilinear")
# 执行重投影
ds_wgs84 = regridder(ds)
# 此时数据集的维度为time, y, x，对应规则的lat/lon坐标

方式二：用rioxarray处理投影转换 #

import rioxarray

# 打开数据集并设置输入坐标系
ds = xr.open_dataset("你的NetCDF文件路径.nc").rio.set_crs(input_crs)
# 直接重投影到EPSG:4326
ds_wgs84 = ds.rio.reproject("EPSG:4326")
# 重投影后的数据集会自动用规则的lon/lat作为坐标维度

注：仅添加lon/lat坐标（不替换维度） #

如果只是想保留原网格，仅把lon/lat作为附加坐标，可直接用assign_coords，但这种方式无法替换rlat/rlon维度：

ds = ds.assign_coords(lon=ds.lon, lat=ds.lat)

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