基于R的terra包重投影EURO CORDEX数据及提取值NA问题排查

一、基于EPSG代码的重投影方法（替代废弃的PROJ4/rgdal） #

必备库 #

• terra：当前R中处理栅格数据的核心库，替代rgdal的重投影功能
• sf：处理矢量数据的CRS转换，与terra完全兼容

正确重投影流程 #

EURO CORDEX的EUR-11数据集采用旋转极坐标投影，无官方EPSG代码，但可以通过以下两种方式实现栅格与矢量的CRS统一：

1. 将栅格重投影到矢量的EPSG（快速提取场景）

你的矢量点图层是WGS84（EPSG:4326），直接将栅格转到此坐标系：

library(terra)
library(sf)

# 读取数据
test <- rast(".../tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_sem_200012-200511.nc", subds="tas")
analysis_loc <- st_read("... .shp")

# 重投影栅格到EPSG:4326
test_wgs84 <- project(test, "EPSG:4326")

# 提取点对应栅格值
ext <- terra::extract(test_wgs84, analysis_loc)

2. 将矢量重投影到栅格的原始CRS（保留栅格精度场景）

避免栅格重投影带来的精度损失，将矢量点转换为栅格自带的旋转投影：

# 获取栅格原始CRS
raster_crs <- crs(test)

# 转换矢量到栅格CRS
analysis_loc_proj <- st_transform(analysis_loc, raster_crs)

# 提取值
ext <- terra::extract(test, analysis_loc_proj)

二、terra::extract返回全NA的原因及修复 #

绘图时栅格与点看似对齐，但提取全NA的核心问题是坐标数值不匹配：
你手动指定cor_crs投影栅格，但矢量点仍保留WGS84的经纬度数值。plot函数会自动做临时可视化对齐，但实际点的坐标并不在重投影后栅格的数值范围内，导致extract无法匹配到对应栅格单元。

修复步骤 #

1. 按上述方法统一栅格与矢量的CRS（二选一即可）
2. 验证坐标范围，确保点落在栅格内：

   # 查看栅格范围
   print(ext(test))
   # 查看转换后点的坐标
   print(st_coordinates(analysis_loc_proj))
   

3. 重新运行extract即可得到非NA值。

关键注意事项 #

• 不要手动编写PROJ4字符串，直接用crs(test)获取栅格自带的CRS参数，避免参数错误导致的坐标偏差
• 若NetCDF数据未自带CRS，可通过Copernicus CDS的数据集文档获取正确的投影参数，手动赋值给栅格：crs(test) <- "正确的PROJ字符串"

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