使用terra::intersect计算相交多边形分区域面积的问题咨询

参照terra::intersect示例操作，代码如下：

library(terra)
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
e <- ext(5.6, 6, 49.55, 49.7)
x <- intersect(v, e)
p <- vect(c("POLYGON ((5.8 49.8, 6 49.9, 6.15 49.8, 6 49.6, 5.8 49.8))",
            "POLYGON ((6.3 49.9, 6.2 49.7, 6.3 49.6, 6.5 49.8, 6.3 49.9))"), crs=crs(v))
values(p) <- data.frame(pid=1:2, area=expanse(p))
y <- intersect(v, p)

需求是统计每个p多边形分别位于Diekirch、Grevenmacher、Luxembourg区域的面积（单位：公顷）。尝试以下代码后，结果和expanse(v)一致，无法得到目标数据：

lapply(y, FUN=expanse, unit="ha")
expanse(v)

yy <- union(v,p)
lapply(yy, FUN=expanse, unit="ha")

你之前的操作没有利用intersect返回结果的属性信息，也未对交集后的单个多边形面积做分组汇总。y <- intersect(v,p)生成的是v和p的重叠区域多边形，每个多边形同时携带v的区域属性（比如NAME_1，对应Diekirch等区域）和p的pid属性，直接基于这些属性分组计算面积即可：

方法1：使用Base R #

# 给每个交集多边形添加面积字段（公顷）
y$area_ha <- expanse(y, unit = "ha")

# 按p的pid和v的区域名称分组汇总总面积
area_summary <- aggregate(area_ha ~ pid + NAME_1, data = y, sum)
print(area_summary)

方法2：使用dplyr（需先安装包） #

library(dplyr)

# 转换为数据框后分组计算
area_summary <- as.data.frame(y) %>%
  mutate(area_ha = expanse(y, unit = "ha")) %>%
  group_by(pid, NAME_1) %>%
  summarise(total_area_ha = sum(area_ha), .groups = "drop")

print(area_summary)

运行后会得到类似如下的结果，清晰展示每个p多边形在三个区域的分布面积：

pid    NAME_1 total_area_ha
1   1   Diekirch      XXXXXXXX
2   1 Luxembourg      XXXXXXXX
3   2 Grevenmacher      XXXXXXXX

内容的提问来源于stack exchange，提问作者we need a Mat. Stat.