如何在R中仅对荷兰国境内的栅格单元完成逐像元值相乘运算

依赖包 #

所需工具包为terra（SpatRaster对象处理核心包），可选geodata用于快速获取荷兰行政区划边界，若已有本地边界文件可忽略geodata。

# 安装依赖（首次运行执行）
install.packages(c("terra", "geodata"))
# 加载包
library(terra)
library(geodata)

实现步骤 #

• 第一步：加载并对齐荷兰国界矢量数据

# 方式1：在线获取荷兰国家边界（level=0对应国界，path为临时存储路径）
nld_border <- gadm(country = "NLD", level = 0, path = tempdir())
# 方式2：读取本地shp边界文件（替换为你自己的文件路径即可）
# nld_border <- vect("your_local_nld_border.shp")

# 对齐边界与栅格的坐标系（你现有栅格为WGS84，无特殊情况不需要额外转换）
nld_border <- project(nld_border, crs(CMAWRF[[1]]))

# 可选：验证两个栅格空间属性完全匹配，返回TRUE代表可直接运算
compareGeom(CMAWRF[[1]], pop_rast)

注：上述代码中pop_rast为你自己的1层人口栅格对象名，可替换为实际命名

• 第二步：生成荷兰区域掩膜，仅保留境内栅格值

# 对人口栅格做掩膜，荷兰境外单元值设为NA
pop_masked <- mask(pop_rast, nld_border)

• 第三步：逐单元逐图层执行相乘运算
  terra的栅格算术运算会自动按位置对齐单元，1层的人口掩膜栅格会自动循环匹配183层的风速栅格，境外NA单元运算后仍为NA，相当于不处理：

# 输出结果为183层SpatRaster，仅荷兰境内单元为「风速*人口」的计算结果
wind_pop_result <- CMAWRF[[1]] * pop_masked

可选验证 #

你可以通过以下操作确认结果符合预期：

• 运行plot(wind_pop_result[[1]])绘制第一层结果，叠加国界查看计算范围是否正确：lines(nld_border, col="red")
• 若需要将境外NA替换为原始风速值/其他固定值，可使用ifel(is.na(wind_pop_result), CMAWRF[[1]], wind_pop_result)自定义调整

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