R语言计算加拿大多边形唯一地理质心的问题求助

问题原因 #

你遇到的问题是因为下载的StatCan shapefile中包含多部件多边形（MULTIPOLYGON）——也就是一个ADA区域由多个不连续的多边形部件组成。默认情况下，st_centroid或st_point_on_surface会为每个部件生成一个点，导致单个ADA区域显示多个质心。

解决步骤 #

1. 检查几何类型：先确认数据中是否存在多部件多边形

table(st_geometry_type(ada))

如果输出包含MULTIPOLYGON，说明确实存在多部件对象。

2. 按特征计算质心：使用st_centroid的by_feature = TRUE参数（默认即为该值，确保对每个ADA特征计算质心而非单个部件），同时建议先转换到平面坐标系（比如你用的32617）计算质心，结果会更准确，之后再转回WGS84。

3. 清理冗余代码：你原代码中重复定义了shapefile_1和sf_cent，这会导致逻辑混乱，需要移除重复部分。

修正后的完整代码 #

library(dplyr)
library(sf)
library(ggplot2)

# 下载并解压shapefile
url_1 <- "https://www12.statcan.gc.ca/census-recensement/alternative_alternatif.cfm?l=eng&dispext=zip&teng=lada000b21a_e.zip&k=%20%20%20151162&loc=//www12.statcan.gc.ca/census-recensement/2021/geo/sip-pis/boundary-limites/files-fichiers/lada000b21a_e.zip"
download.file(url_1, destfile = "lada000b21a_e.zip")
unzip("lada000b21a_e.zip")

# 读取数据
ada <- st_read("lada000b21a_e.shp")

# 转换到平面坐标系计算质心（确保质心计算准确），再转回WGS84
sf_cent <- ada %>%
  st_transform(32617) %>%
  st_centroid(by_feature = TRUE) %>%  # 每个ADA特征对应一个质心
  st_transform(4326)

# 可视化：多边形用平面坐标系，质心用WGS84（ggplot会自动处理CRS匹配）
ggplot() + 
  geom_sf(data = ada %>% st_transform(32617), fill = 'white') +
  geom_sf(data = sf_cent, color = 'red', size = 1)

额外说明 #

如果某些ADA区域的多部件多边形确实需要合并为单个多边形（比如后续分析需要），可以用以下代码先合并：

ada_single <- ada %>%
  group_by(ADAUID) %>%
  summarise(geometry = st_union(geometry), .groups = "drop")

之后再对ada_single计算质心即可。

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