sf包下非完美嵌套多边形按面积多数规则空间连接方案咨询

我有两个shapefile文件，目前遇到两个问题：

• 二者应为嵌套多边形结构，但三级（level 3）小多边形无法完全匹配落入二级（level 2）大多边形内部。
  
• 即便二者嵌套完全匹配，也没有可用的关联键：小多边形文件缺失可对应大多边形文件中level2_code字段的同名字段。

我目前使用sf进行空间操作，想要将三级小多边形以1:1的关系匹配到二级大多边形上，采用多数占比规则：如果一个小多边形被两个大多边形分割为60/40的面积占比，则将其匹配到覆盖其60%面积的大多边形下，请问最优的空间连接方案是什么？

示例数据 #

示例shapefile获取代码：

devtools::source_gist("3bf2ed930868a2807a58b5c04f700a1b")

二级大多边形数据（example_bigger） #

Simple feature collection with 27 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 6.63908 ymin: 5.194552 xmax: 7.432286 ymax: 5.938931
Geodetic CRS:  WGS 84
# A tibble: 27 x 2
   level2_code                                                             geometry
         <dbl>                                                   <MULTIPOLYGON [°]>
 1           1 (((7.321909 5.40493, 7.321575 5.401201, 7.322123 5.399188, 7.322649…
 2           2 (((7.308938 5.617926, 7.308858 5.613426, 7.308925 5.611024, 7.30905…
 3           3 (((7.361156 5.714312, 7.358485 5.713518, 7.354733 5.711535, 7.35242…
 4           4 (((7.388033 5.510387, 7.387133 5.507236, 7.385656 5.503802, 7.38476…
 5           5 (((7.178212 5.808309, 7.178176 5.806321, 7.17929 5.803731, 7.18097 …
 6           6 (((7.202845 5.920962, 7.206038 5.924097, 7.20778 5.925807, 7.209518…
 7           7 (((7.015881 5.895269, 7.015242 5.890107, 7.012308 5.886351, 7.01131…
 8           8 (((7.381608 5.705609, 7.382155 5.703596, 7.382979 5.70072, 7.384001…
 9           9 (((7.209132 5.60738, 7.211122 5.605943, 7.213938 5.604008, 7.216109…
10          10 (((7.265981 5.744175, 7.26594 5.740849, 7.265327 5.737696, 7.265033…
# … with 17 more rows

三级小多边形数据（example_smaller） #

Simple feature collection with 348 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 6.63908 ymin: 5.194552 xmax: 7.432286 ymax: 5.938931
Geodetic CRS:  WGS 84
# A tibble: 348 x 2
   level3_code                                                             geometry
         <dbl>                                                   <MULTIPOLYGON [°]>
 1           1 (((7.153985 5.786385, 7.151873 5.777959, 7.155913 5.768072, 7.14011…
 2           2 (((7.353153 5.870928, 7.34829 5.856211, 7.342018 5.842341, 7.314459…
 3           3 (((7.235717 5.510277, 7.235989 5.5113, 7.236185 5.512454, 7.23656 5…
 4           4 (((6.947564 5.754671, 6.946301 5.755986, 6.945236 5.756682, 6.95680…
 5           5 (((7.001083 5.550417, 7.00185 5.55142, 7.002059 5.552684, 7.002277 …
 6           6 (((7.043065 5.77843, 7.04189 5.777916, 7.041544 5.777805, 7.041274 …
 7           7 (((7.002059 5.552684, 7.00185 5.55142, 7.001083 5.550417, 6.999884 …
 8           8 (((7.262075 5.526829, 7.262966 5.528058, 7.264392 5.528519, 7.26564…
 9           9 (((7.07199 5.735022, 7.070997 5.735426, 7.063343 5.736378, 7.055127…
10          10 (((7.265549 5.575613, 7.265963 5.574931, 7.266252 5.574365, 7.26567…
# … with 338 more rows

最优实现方案 #

通过面积占比加权的空间连接即可实现需求，完整可运行代码如下：

library(sf)
library(dplyr)

# 1. 计算小多边形与相交的二级多边形的重叠面积
overlap_stats <- st_intersection(example_smaller, example_bigger) %>%
  mutate(overlap_area = st_area(.)) %>%
  st_drop_geometry()

# 2. 按小多边形分组，取重叠面积最大的二级多边形编码
match_rule <- overlap_stats %>%
  group_by(level3_code) %>%
  arrange(desc(overlap_area)) %>%
  slice(1) %>%
  ungroup() %>%
  select(level3_code, level2_code)

# 3. 匹配结果关联回原始小多边形数据
small_poly_with_level2 <- example_smaller %>%
  left_join(match_rule, by = "level3_code")

方案说明 #

• 完全符合多数占比匹配规则，自动将每个小多边形关联到重叠面积占比最高的二级多边形下
• 运行效率较高，st_intersection会自动跳过无相交关系的多边形组合，不需要全量计算
• 异常值可追溯：若存在小多边形未与任何二级多边形相交的情况，对应level2_code会返回NA，可单独筛选处理这类边界异常数据

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