如何在R中匹配合并两类地理Shapefile文件

数据准备 #

已通过以下代码导入两类Shapefile并统一转换为WGS84坐标系：

Bigger_Units  <- sf::st_read("C:/Users/ME/OneDrive/Documents/hr_shape/HR_000a18a_e.shp", options = "ENCODING=WINDOWS-1252") %>%  
  sf::st_transform('+proj=longlat +datum=WGS84')

Smaller_Units <- sf::st_read("C:/Users/ME/OneDrive/Documents/shape7/lfsa000b16a_e.shp", options = "ENCODING=WINDOWS-1252") %>%  
  sf::st_transform('+proj=longlat +datum=WGS84')

数据结构示例 #

小单元数据（Smaller_Units）

> head(Smaller_Units)
Simple feature collection with 6 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -60.67924 ymin: 45.72791 xmax: -59.93641 ymax: 46.26758
CRS:           +proj=longlat +datum=WGS84
  CFSAUID PRUID                        PRNAME                       geometry
1     B1E    12 Nova Scotia / Nouvelle-Écosse MULTIPOLYGON (((-60.00359 4...
2     B1G    12 Nova Scotia / Nouvelle-Écosse MULTIPOLYGON (((-60.00298 4...
3     B1H    12 Nova Scotia / Nouvelle-Écosse MULTIPOLYGON (((-60.04267 4...
4     B1J    12 Nova Scotia / Nouvelle-Écosse MULTIPOLYGON (((-60.2913 46...
5     B1K    12 Nova Scotia / Nouvelle-Écosse MULTIPOLYGON (((-59.99001 4...
6     B1L    12 Nova Scotia / Nouvelle-Écosse MULTIPOLYGON (((-60.1559 46...

大单元数据（Bigger_Units）

> head(Bigger_Units)
Simple feature collection with 6 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -82.80346 ymin: 42.05986 xmax: -76.16181 ymax: 45.43812
CRS:           +proj=longlat +datum=WGS84
  HR_UID                                                  ENGNAME                                                                 FRENAME
1   3537                             City of Hamilton Health Unit                       Circonscription sanitaire de la cité de Hamilton
2   3538          Hastings and Prince Edward Counties Health Unit      Circonscription sanitaire des comtés de Hastings et Prince Edward
3   3539                                 Huron County Health Unit                            Circonscription sanitaire du comté de Huron
4   3540                                 Chatham-Kent Health Unit                               Circonscription sanitaire de Chatham-Kent
5   3541 Kingston, Frontenac and Lennox and Addington Health Unit Circonscription sanitaire de Kingston, Frontenac et Lennox et Addington
6   3542                                      Lambton Health Unit                                    Circonscription sanitaire de Lambton
  SHAPE_AREA SHAPE_LEN                       geometry
1 1212618763  173963.7 MULTIPOLYGON (((-79.86045 4...
2 9131935097  569876.5 MULTIPOLYGON (((-77.69017 4...
3 3696168189  335661.7 MULTIPOLYGON (((-80.97369 4...
4 3078455510  269377.6 MULTIPOLYGON (((-81.8351 42...
5 8174790022  460742.3 MULTIPOLYGON (((-76.79163 4...
6 4172673622  398518.1 MULTIPOLYGON (((-81.76871 4...

地理单元匹配方法 #

利用sf包的空间连接功能，可快速实现小单元（CFSAUID）与大单元（ENGNAME）的匹配，根据空间关系选择对应规则：

1. 完全包含匹配（优先推荐） #

如果小单元完全被大单元覆盖，使用st_within作为空间判断规则：

library(dplyr)

# 左连接保留所有小单元，匹配对应的大单元
match_result <- Smaller_Units %>%
  sf::st_join(Bigger_Units, join = sf::st_within) %>%
  # 提取目标字段
  select(smaller_units = CFSAUID, corresponding_bigger_unit = ENGNAME) %>%
  # 移除几何列，转为普通数据框
  sf::st_drop_geometry()

2. 相交匹配（处理部分重叠场景） #

如果小单元与大单元存在部分重叠，改用st_intersects规则：

match_result <- Smaller_Units %>%
  sf::st_join(Bigger_Units, join = sf::st_intersects) %>%
  select(smaller_units = CFSAUID, corresponding_bigger_unit = ENGNAME) %>%
  sf::st_drop_geometry()

处理小单元多归属情况 #

当一个小单元同时属于多个大单元时，可根据需求选择以下处理方式：

1. 保留所有匹配结果 #

直接保留多行记录，完整展示小单元的多归属关系：

# 结果示例：
# smaller_units    corresponding_bigger_unit
# 1           BXX City of Hamilton Health Unit
# 2           BXX     Huron County Health Unit

2. 保留面积占比最高的匹配 #

为每个小单元仅保留关联度最高的大单元（按交集面积占比判断）：

# 计算小单元与大单元的交集面积，筛选占比最高的匹配
intersect_areas <- sf::st_intersection(Smaller_Units, Bigger_Units) %>%
  mutate(intersect_area = sf::st_area(.)) %>%
  group_by(CFSAUID) %>%
  filter(intersect_area == max(intersect_area)) %>%
  ungroup() %>%
  select(smaller_units = CFSAUID, corresponding_bigger_unit = ENGNAME) %>%
  sf::st_drop_geometry()

3. 合并多匹配结果为一行 #

将小单元对应的多个大单元用分隔符合并为字符串：

merged_result <- match_result %>%
  group_by(smaller_units) %>%
  summarise(corresponding_bigger_unit = paste(corresponding_bigger_unit, collapse = ", ")) %>%
  ungroup()

输出目标格式 #

经过上述处理后，即可得到符合需求的匹配结果：

> head(match_result)
 smaller_units    corresponding_bigger_unit
1           B1E City of Hamilton Health Unit
2           B1G City of Hamilton Health Unit
3           B1H City of Hamilton Health Unit
4           B1J     Huron County Health Unit
5           B1K     Huron County Health Unit
6           B1L     Huron County Health Unit

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