R中基于WorldClim栅格统计秘鲁196省平均温极值的高效方法

目前你已经完成了秘鲁范围的栅格裁剪，接下来执行**矢量分区栅格值统计（Zonal Statistics）**操作，即可快速得到每个省份对应的气温均值，最终生成符合要求的196行data.frame。以下是两种经过性能优化的实现路径：

优先选择：Terra包实现（计算效率最高） #

Terra是R语言中替代旧版raster包的新一代栅格处理工具，核心计算逻辑用C语言实现，同等数据量下比旧包快3-5倍，尤其适合1km分辨率的大尺度栅格处理，实现代码如下：

# 加载依赖包
library(terra)
library(sf)

# 读入矢量和栅格数据
my_shape <- vect("multipolygon shapfile with 196 provinces inside Peru")
r_max <- rast("最高气温栅格文件路径")
r_min <- rast("最低气温栅格文件路径")

# *关键校验步骤*：对齐矢量与栅格的坐标参考系统，避免统计偏差
my_shape <- project(my_shape, crs(r_max))

# 将最高温、最低温栅格堆叠为单个多波段栅格，实现一次计算出两个变量
r_temp_stack <- c(r_max, r_min)

# 裁剪栅格到秘鲁边界范围，同时掩膜掉省界外的无效值（比单纯crop精度更高）
r_crop <- crop(r_temp_stack, my_shape, mask = TRUE)

# 核心统计步骤：按196个省的矢量边界，计算每个区域覆盖栅格的均值
province_temp_stats <- zonal(
  x = r_crop,
  z = my_shape,
  fun = "mean",
  na.rm = TRUE
)

# 将统计结果与省份属性（省名、省ID等）拼接，转为196行的data.frame
result_df <- cbind(as.data.frame(my_shape), province_temp_stats)

兼容旧代码：raster包实现 #

如果你已经沿用旧的raster包写了部分流程，可以通过extract函数实现同等功能，不需要重写全部已有代码：

# 加载依赖包
library(raster)
library(sf)

# 读入已有数据流程
my_shape <- st_read("multipolygon shapfile with 196 provinces inside Peru")
r <- raster("raster file of temp with 1km2 grid")
r2 <- crop(r, extent(my_shape))

# 对齐投影
my_shape <- st_transform(my_shape, crs = st_crs(r))

# 掩膜掉省界外的栅格值
r2 <- mask(r2, my_shape)

# 按矢量边界计算每个省的栅格均值
province_temp_stats <- extract(
  x = r2,
  y = my_shape,
  fun = mean,
  na.rm = TRUE,
  df = TRUE
)

# 拼接属性列得到最终data.frame
result_df <- cbind(st_drop_geometry(my_shape), province_temp_stats[, -1])

效率优化注意点 #

• 不要用逐省循环的方式提取栅格值，内置的zonal/extract函数是底层C实现的计算逻辑，比手写R循环效率高一个数量级
• 同时计算最高温、最低温时，尽量将两个栅格堆叠为多波段栅格一次统计，不要分开两次运行计算流程，可节省一半运行时间
• 统计时务必开启na.rm = TRUE参数，排除省界边缘覆盖的空值、无效值对均值的干扰
• 相比单纯做crop裁剪，增加mask步骤去掉省界外的无效栅格，可减少后续统计阶段的无效计算量

内容的提问来源于stack exchange，提问作者Alfonso R.