如何解决NOAA EPSG:9001栅格重投影至EPSG:4326的报错问题？

问题背景 #

处理NOAA全球温度栅格数据时，需将EPSG:9001投影（经度范围0-360）转换为EPSG:4326投影（适配世界国家矢量数据，经度范围-180-180），执行projectRaster时出现报错：

Error in if (yn == yx) { : missing value where TRUE/FALSE needed
In addition: Warning message:
In `dim<-`(`*tmp*`, value = c(nr, nc)) :
  NAs introduced by coercion to integer range

原始栅格范围：

> extent(raster.nc)
class      : Extent 
xmin       : 0 
xmax       : 360 
ymin       : -90 
ymax       : 90 

矢量数据范围：

extent(world.shape)
class      : Extent 
xmin       : -180 
xmax       : 180 
ymin       : -89 
ymax       : 83.6236 

解决方案 #

1. 先调整经度范围至-180-180 #

报错核心是0-360与-180-180的经度范围冲突，先通过rotate()函数转换栅格：

# 转换经度范围并调整栅格数据排列
raster.nc_180 <- rotate(raster.nc)

2. 执行投影转换 #

调整范围后，直接指定EPSG:4326的CRS进行投影：

# 明确目标投影
target_crs <- CRS(SRS_string = "EPSG:4326")
# 或用完整字符串："+proj=longlat +datum=WGS84 +no_defs"

# 执行双线性插值投影
raster_proj <- projectRaster(raster.nc_180, crs = target_crs, method = "bilinear")

3. 验证结果 #

检查转换后栅格的范围是否匹配：

extent(raster_proj)

正常会输出经度-180至180、纬度-90至90的范围，此时可与EPSG:4326的矢量数据正常匹配。

额外排查点 #

• 若原始栅格未正确定义CRS，先手动指定EPSG:9001的完整参数：

crs(raster.nc) <- "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

• 避免通过projection(world.shape)获取CRS，直接指定EPSG代码更稳定，防止矢量数据的CRS字符串格式异常。

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