如何用单层Shapefile掩膜多层NetCDF？输出无穷值问题排查

1. 验证空间重叠性 #

先确认栅格与Shapefile存在空间交集，否则掩膜后必然全为NA（被writeRaster转成无穷大）：

library(rgeos)
# 转换栅格范围为多边形
raster_ext <- as(extent(prcp_2000), "SpatialPolygons")
# 检查是否重叠
print(gIntersects(raster_ext, shapefile))

若返回FALSE，说明坐标系转换错误或数据范围完全不匹配，需重新核对CRS设置或数据来源。

2. 确保CRS完全一致 #

即使执行了spTransform，仍可能存在投影参数细微差异（如椭球体、单位），需强制匹配：

# 先确认栅格CRS是否正确识别
print(crs(prcp_2000))
# 若识别错误，手动指定（示例为WGS84，替换为实际CRS）
crs(prcp_2000) <- "+proj=longlat +datum=WGS84 +no_defs"
# 重新转换Shapefile
shapefile <- spTransform(shapefile, crs(prcp_2000))
# 验证一致性
print(identical(crs(shapefile), crs(prcp_2000)))

3. 调整掩膜与裁剪顺序 #

先裁剪再掩膜，既减少计算量，也能提前发现无重叠问题：

prcp_2000 <- crop(prcp_2000, shapefile)
prcp_2000 <- mask(prcp_2000, shapefile)

4. 修复Shapefile无效几何 #

无效多边形会导致掩膜操作无法正确识别重叠区域，需先修复：

shapefile <- gBuffer(shapefile, byid = TRUE, width = 0)

5. 控制NA编码格式 #

writeRaster输出CDF时默认将NA转为无穷大，需手动指定合法的缺省值：

writeRaster(prcp_2000, outfile, overwrite=TRUE, format="CDF", 
            varname="prcp", varunit="mm/day", 
            longname="mm of precipitation per day", 
            xname="lon", yname="lat", zname="day", 
            zunit="days since 1900-01-01",
            NAflag = -9999) # 用常用的-9999替代无穷大

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