R中raster包提取NetCDF数据的投影匹配与坐标问题咨询

问题2解答：坐标不匹配的原因 #

你对投影信息的理解存在两个误区：

1. 混淆了大地基准面和坐标系：proj4string里的datum=WGS84仅代表该投影使用WGS84作为大地基准面，串开头明确标注+proj=lcc，说明栅格使用的是兰伯特等角圆锥投影坐标系，属于平面投影坐标系，单位是米（串内+units=m也佐证了这一点）。你看到的111500、571500这类数值是投影后的平面米制坐标，完全符合LCC投影的特征，不存在格式错误。
2. 你构造坐标点时存在两个低级错误：一是经纬度填反了——该数据集覆盖中欧区域，投影中心经度为参数标注的13.33°E，正常点位经度应该在10-16区间、纬度在46-49区间，你写的lon=47.4, lat=11.5对应的位置在索马里以东的印度洋上，根本不在数据集覆盖范围内；二是没有给坐标点设置对应WGS84经纬度的坐标系信息，相当于直接把经纬度数值当成米制平面坐标去匹配栅格，点位直接飘到栅格范围外，自然返回NA。

问题1解答：正确的坐标输入与提取方法 #

核心逻辑是先给输入的经纬度点绑定正确的WGS84地理坐标系信息，再将点投影转换为和栅格一致的LCC坐标系，再做提取即可，可直接运行的代码如下：

library(raster)
library(sp)

# 读入NC栅格数据
b <- brick("precipitationintensity14rcp45modelsensemblemedianannual.nc")

# 构造目标点位坐标，注意：第一个参数是经度（lon）、第二个是纬度（lat），不要写反
coords <- data.frame(lon = 11.5, lat = 47.4)
coordinates(coords) <- c("lon", "lat")
# 给输入点设置WGS84经纬度坐标系（EPSG:4326是通用WGS84编码）
proj4string(coords) <- CRS("+init=epsg:4326")

# 将点转换为栅格对应的LCC投影坐标系
coords_trans <- spTransform(coords, CRS(proj4string(b)))

# 执行值提取
extract(b$X1981.07.15, coords_trans)

如果后续需要处理大量NC文件，更推荐用替代raster的新版terra包，坐标转换和提取的逻辑完全一致，运算速度更快。

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