如何用R语言terra包将旋转极NetCDF转换为标准经纬度网格

我有一份覆盖北美的气候投影文件，它的范围不是标准经纬度格式，需要转换为标准经纬度坐标（经度-180至180，纬度-90至90）。我曾考虑过偏移范围，但觉得这个方法不对。请问用R语言的terra包完成这个转换的最佳方式是什么？

data <- terra::rast("snw_NAM-44_CCCma-CanESM2_historical-r1_r1i1p1_CCCma-CanRCM4_r2_mon_195001-195012.nc")
data
#> class       : SpatRaster 
#> dimensions  : 130, 155, 12  (nrow, ncol, nlyr)
#> resolution  : 0.4400001, 0.44  (x, y)
#> extent      : -34.1, 34.1, -28.82, 28.38  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source      : snw_NAM-44_CCCma-CanESM2_historical-r1_r1i1p1_CCCma-CanRCM4_r2_mon_195001-195012.nc:snw 
#> varname     : snw (Surface Snow Amount) 
#> names       :  snw_1,  snw_2,  snw_3,  snw_4,  snw_5,  snw_6, ... 
#> unit        : kg m-2, kg m-2, kg m-2, kg m-2, kg m-2, kg m-2, ... 
#> time (days) : 1954-08-08 to 1955-07-08

以下是这份带旋转极的NetCDF文件详情：

ncdf4::nc_open("snw_NAM-44_CCCma-CanESM2_historical-r1_r1i1p1_CCCma-CanRCM4_r2_mon_195001-195012.nc")
#> File snw_NAM-44_CCCma-CanESM2_historical-r1_r1i1p1_CCCma-CanRCM4_r2_mon_195001-195012.nc (NC_FORMAT_CLASSIC):
#> 
#>      5 variables (excluding dimension variables):
#>         double time_bnds[bnds,time]   
#>         double lon[rlon,rlat]   
#>             long_name: longitude
#>             units: degrees_east
#>         double lat[rlon,rlat]   
#>             long_name: latitude
#>             units: degrees_north
#>         char rotated_pole[]   
#>             grid_mapping_name: rotated_latitude_longitude
#>             grid_north_pole_longitude: 83
#>             grid_north_pole_latitude: 42.5
#>         float snw[rlon,rlat,time]   
#>             long_name: Surface Snow Amount
#>             standard_name: surface_snow_amount
#>             units: kg m-2
#>             cell_methods: time:mean (interval:1200 seconds)
#>             coordinates: lon lat
#>             grid_mapping: rotated_pole
#>             _FillValue: 1.00000002004088e+20
#> 
#>      4 dimensions:
#>         time  Size:12   *** is unlimited *** 
#>             long_name: time
#>             standard_name: time
#>             units: days since 1949-12-01
#>             axis: T
#>             calendar: 365_day
#>             bounds: time_bnds
#>         rlon  Size:155 
#>             long_name: longitude in rotated pole grid
#>             units: degrees
#>             axis: X
#>             standard_name: grid_longitude
#>         rlat  Size:130 
#>             long_name: latitude in rotated pole grid
#>             units: degrees
#>             axis: Y
#>             standard_name: grid_latitude
#>         bnds  Size:2 (no dimvar)
#> 
#>     26 global attributes:
#>         title: CanRCM4 model output prepared for CanSISE Project
#>         institution: CCCma (Canadian Centre for Climate Modelling and Analysis, Victoria, BC, Canada)
#>         institute_id: CCCma
#>         experiment: CanSISE downscaling run driven by CCCma-CanESM2 eia-001
#>         experiment_id: historical-r1
#>         driving_experiment: CCCma-CanESM2, historical-r1, r1i1p1
#>         driving_model_id: CCCma-CanESM2
#>         driving_experiment_name: historical-r1
#>         driving_model_ensemble_member: r1i1p1
#>         realization: 1
#>         initialization_method: 1
#>         physics_version: 1
#>         forcing: GHG,Oz,SA,BC,OC,LU,Vl (GHG includes CO2,CH4,N2O,CFC11,effective CFC12)
#>         project_id: CanSISE
#>         model_id: CCCma-CanRCM4
#>         CORDEX_domain: NAM-44
#>         rcm_version_id: r2
#>         frequency: mon
#>         product: output
#>         CCCma_runid: nam44_v001_eia-001
#>         Conventions: CF-1.4
#>         creation_date: 2016-10-04-T17:30:32Z
#>         contact: cccma_info@ec.gc.ca
#>         references: http://www.cccma.ec.gc.ca/models
#>         history: created: 2016-10-04 17:31:06 by rcm2nc

解决方案 #

这份文件采用旋转极投影（rotated_latitude_longitude），极点坐标为（83°E，42.5°N），直接偏移范围的方法无效，需要通过投影转换或利用内置经纬度变量来实现坐标转换，以下是两种可行方法：

方法一：利用旋转投影CRS进行转换（规则网格适用）

先根据文件中的旋转极信息定义投影坐标系，再用terra::project()转换到标准WGS84经纬度：

library(terra)

# 读取数据
data <- rast("snw_NAM-44_CCCma-CanESM2_historical-r1_r1i1p1_CCCma-CanRCM4_r2_mon_195001-195012.nc")

# 定义旋转极投影的PROJ字符串
rot_crs <- "+proj=ob_tran +o_proj=latlon +lon_0=83 +o_lat_p=42.5 +ellps=WGS84 +no_defs"
# 为栅格设置正确的旋转投影CRS
crs(data) <- rot_crs

# 转换到标准WGS84经纬度（EPSG:4326）
data_wgs84 <- project(data, "EPSG:4326")

# 查看转换结果
data_wgs84

方法二：利用内置lon/lat变量构建栅格（不规则网格适用）

如果旋转投影的网格不规则，直接使用文件中存储的每个网格点的经纬度坐标，通过插值构建标准经纬度栅格：

library(terra)
library(ncdf4)

# 打开NetCDF文件
nc <- nc_open("snw_NAM-44_CCCma-CanESM2_historical-r1_r1i1p1_CCCma-CanRCM4_r2_mon_195001-195012.nc")

# 提取经纬度和雪量数据
lon <- ncvar_get(nc, "lon")
lat <- ncvar_get(nc, "lat")
snw <- ncvar_get(nc, "snw")
nc_close(nc)

# 调整雪量数据维度（适配terra的[行,列,时间层]格式）
snw_arr <- aperm(snw, c(2, 1, 3))

# 创建包含经纬度的点矢量
pts <- vect(cbind(as.vector(lon), as.vector(lat)), crs = "EPSG:4326")
# 为点矢量添加雪量属性
values(pts) <- as.data.frame(snw_arr)

# 创建目标栅格（覆盖北美区域，分辨率与原数据一致）
target_rast <- rast(ext(-180, -50, 15, 85), res = 0.44, crs = "EPSG:4326")

# 双线性插值到目标栅格
data_wgs84 <- rasterize(pts, target_rast, values(pts), method = "bilinear")

方法选择 #

• 方法一速度快，适合规则网格的旋转投影数据，精度可靠；
• 方法二更灵活，适用于不规则网格，或旋转投影CRS定义存在问题的场景，结果更贴合原始数据的经纬度信息。

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