R语言从栅格Raster Brick中按经纬度裁剪指定区域的问题求助

核心问题原因 #

• 旧版raster包对sf格式的矢量对象兼容性不足，直接传入crop()函数会出现识别异常
• 手动给栅格赋值CRS容易出现投影匹配偏差，优先读取原始数据自带的投影信息校验

基于原有raster包的修正代码 #

url <- "http://thredds.cdip.ucsd.edu/thredds/fileServer/cdip/model/MOP_grids/CA_0.01_nowcast.nc"
options(timeout = 1000)
# download.file返回值为状态码，无需赋值给数据变量
download.file(url, "/Users/mycomp/Desktop/wave_data.nc")
data_set <- "/Users/mycomp/Desktop/wave_data.nc"
waves <- brick(data_set, sub = "waveHs")
# 先校验原始栅格自带的CRS，确认和ROI的4326一致再后续操作
print(crs(waves))
ROI <- st_polygon(list(rbind(c(-121.0062, 33.10625), c(-121.0062, 34.90625), 
                  c(-118.7438, 34.90625), c(-118.7438, 33.10625), c(-121.0062, 33.10625))))
ROI <- st_sfc(ROI, crs = 4326)
# 把sf矢量转成raster包可识别的Extent对象
roi_extent <- extent(ROI)
# 提取第一时间层
t_1 <- subset(waves, 1)
# 先按范围裁剪
t_1_crop <- crop(t_1, roi_extent)
# 再按多边形掩膜，把ROI外的数值设为NA
t_1_mask <- mask(t_1_crop, ROI)

更推荐的terra包实现方案（raster包已停止维护更新） #

library(terra)
library(sf)
data_set <- "/Users/mycomp/Desktop/wave_data.nc"
# 用terra的rast接口读取nc文件，兼容性和性能都优于raster包
waves <- rast(data_set, subds = "waveHs")
ROI <- st_polygon(list(rbind(c(-121.0062, 33.10625), c(-121.0062, 34.90625), 
                  c(-118.7438, 34.90625), c(-118.7438, 33.10625), c(-121.0062, 33.10625))))
ROI <- st_sfc(ROI, crs = 4326)
# terra的crop原生支持sf对象，设置mask=TRUE可一步完成范围裁剪+区域掩膜
t_1_crop <- crop(waves[[1]], ROI, mask = TRUE)

额外注意事项 #

• 如果裁剪后还是空值，检查ROI的坐标顺序：sf多边形要求坐标为经度在前，纬度在后，且闭合方向符合矢量规范
• 可先调用plot(waves[[1]])查看原始栅格的坐标范围，确认ROI落在原始栅格的覆盖区间内

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