如何使用R语言terra包按类别和多边形计算分区统计？

我想用R语言的terra包，同时按高程类别和多边形两个维度对高程数据做分区统计。目前已实现单多边形层面的均值计算，现有代码及输出如下：

现有实现代码 #

library(terra)

# 读取栅格文件
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
plot(r)

# 读取多边形文件
pathshp <- system.file("ex/lux.shp", package = "terra")
v <- vect(pathshp)

plot(v, add = T)

# 重分类高程栅格
m_r <- c(-Inf, 200, 1,  
         200, 300, 2,
         300, 400, 3,
         400, 500, 4,
         500, Inf, 5)

rclmat <- matrix(m_r, ncol=3, byrow=TRUE)
rc <- terra::classify(r, rclmat, include.lowest=TRUE)
plot(rc)

# 多边形层面均值计算
mean_elev <- terra::extract(r, v, mean, na.rm=TRUE, bind = T, 
                               exact = T)
as.data.frame(mean_elev)

现有输出 #

ID_1       NAME_1 ID_2           NAME_2 AREA    POP elevation
1     1     Diekirch    1         Clervaux  312  18081  467.3792
2     1     Diekirch    2         Diekirch  218  32543  334.6856
3     1     Diekirch    3          Redange  259  18664  377.2070
4     1     Diekirch    4          Vianden   76   5163  372.2499
5     1     Diekirch    5            Wiltz  263  16735  418.7867
6     2 Grevenmacher    6       Echternach  188  18899  314.7698
7     2 Grevenmacher    7           Remich  129  22366  240.2105
8     2 Grevenmacher   12     Grevenmacher  210  29828  283.2307
9     3   Luxembourg    8         Capellen  185  48187  329.8955
10    3   Luxembourg    9 Esch-sur-Alzette  251 176820  310.3833
11    3   Luxembourg   10       Luxembourg  237 182607  314.0103
12    3   Luxembourg   11           Mersch  233  32112  313.5930

期望输出 #

需要同时按多边形和高程类别统计，输出格式如下（包含多边形属性、高程类别、该类别在多边形内的占比、类别均值）：

解决方案代码 #

结合terra提取数据和dplyr分组统计，实现双维度统计：

library(terra)
library(dplyr)

# 读取栅格与矢量数据
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
pathshp <- system.file("ex/lux.shp", package = "terra")
v <- vect(pathshp)

# 重分类高程栅格
m_r <- c(-Inf, 200, 1,  
         200, 300, 2,
         300, 400, 3,
         400, 500, 4,
         500, Inf, 5)
rclmat <- matrix(m_r, ncol=3, byrow=TRUE)
rc <- classify(r, rclmat, include.lowest=TRUE)

# 合并原始高程与分类栅格，方便一次性提取
r_stack <- c(r, rc)
names(r_stack) <- c("elevation", "elev_class")

# 提取每个多边形内的所有栅格值（不绑定矢量，便于后续分组）
extracted_data <- extract(r_stack, v, exact = TRUE, bind = FALSE)

# 按多边形ID和高程类别分组统计
final_stats <- extracted_data %>%
  # 分组统计每个类别在多边形内的数量和均值
  group_by(ID, elev_class) %>%
  summarise(
    class_count = n(),
    mean_elev = mean(elevation, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  # 计算每个多边形的总栅格数，用于计算占比
  left_join(
    extracted_data %>% group_by(ID) %>% summarise(total_count = n()),
    by = "ID"
  ) %>%
  mutate(
    percentage = round((class_count / total_count) * 100, 1)
  ) %>%
  # 关联多边形的属性信息
  left_join(
    as.data.frame(v),
    by = c("ID" = "ID_2")
  ) %>%
  # 整理列顺序，匹配期望格式
  select(ID_1, NAME_1, ID_2 = ID, NAME_2, elev_class, percentage, mean_elev) %>%
  # 按行政区和类别排序
  arrange(ID_1, ID_2, elev_class)

# 查看结果
print(final_stats)

关键说明 #

1. 栅格合并：将原始高程栅格和分类栅格合并为一个栅格栈，一次提取两个变量，避免重复操作
2. 数据提取：extract设置bind=FALSE，返回数据框格式，更适合用dplyr做分组统计
3. 分组计算：先按多边形ID和高程类别分组，统计类别数量和均值；再计算每个类别在多边形内的占比
4. 属性关联：将统计结果与多边形的属性信息（如行政区名称）关联，输出完整的多维度统计数据

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