在R语言中计算shapefile文件按不同字段分组的多边形表面积

核心思路 #

你的数据已经是sf类矢量对象，直接使用sf包原生的st_area()函数计算椭球面积，搭配dplyr的分组语法即可实现和ArcGIS一致的计算结果，你之前使用Areapolygons失败是因为该函数适配旧的sp类矢量数据，和当前的sf类数据不兼容。

操作步骤 #

1. 加载依赖包 #

library(sf)
library(dplyr)

2. 生成分组面积汇总表 #

如果需要输出每个分组对应的总面积统计结果，使用summarise语法：

• 按国家字段NAME分组：

area_by_country <- PRIO %>%
  group_by(NAME) %>%
  summarise(
    # st_area默认返回平方米，除以1e6转为平方公里
    total_area_km2 = as.numeric(sum(st_area(geometry)) / 1e6),
    .groups = "drop"
  )
# 查看结果
head(area_by_country)

• 按农业生态区字段AEZ分组：
  仅需要修改group_by中的分组字段即可，其余逻辑完全一致：

area_by_aez <- PRIO %>%
  group_by(AEZ) %>%
  summarise(
    total_area_km2 = as.numeric(sum(st_area(geometry)) / 1e6),
    .groups = "drop"
  )

3. 新增分组面积列到原始表 #

如果需要将分组总面积作为新列追加到原始矢量数据的每一行，使用mutate语法：

# 新增国家总面积列
PRIO <- PRIO %>%
  group_by(NAME) %>%
  mutate(country_area_km2 = as.numeric(sum(st_area(geometry)) / 1e6)) %>%
  ungroup()

# 新增AEZ分区总面积列
PRIO <- PRIO %>%
  group_by(AEZ) %>%
  mutate(aez_area_km2 = as.numeric(sum(st_area(geometry)) / 1e6)) %>%
  ungroup()

注意事项 #

• st_area会自动基于WGS84椭球参数计算高精度面积，和ArcGIS的面积计算结果误差小于万分之一，满足一般分析需求
• 如果你的数据量极大，可以先转为等面积投影后再计算，速度会更快，计算结果差异可忽略

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