如何高效将经纬度观测点匹配到对应多边形并生成归属标记

问题背景 #

我有一个包含经纬度观测坐标的dataframe，以及三个覆盖观测区域的多边形。需要在dataframe中新增一列，用简短标识标记每个观测点所属的多边形。当前通过逐个调用sf::st_within判断归属，再用循环生成loc列，代码繁琐，希望找到更高效的一次性实现方法。

示例数据与多边形定义 #

观测数据 #

testdata <- structure(list(lon = c(-79.30315, -79.29561, -79.29572, -79.29833, 
-79.29603, -79.29659, -79.29097, -79.29347, -79.29347, -79.294, 
-79.29644, -79.29653, -79.27027, -79.27656, -79.26938, -79.28683, 
-79.29256, -79.28848, -79.29256, -79.29097, -79.292277, -79.29311, 
-79.29259, -79.292277, -79.29654, -79.29123, -79.3037, -79.28743, 
-79.29641, -79.29661, -79.29295, -79.29101, -79.2951, -79.29652, 
-79.29035, -79.29183, -79.29661, -79.30317, -79.28734, -79.28734, 
-79.296509, -79.24673, -79.24699, -79.24689, -79.24693, -79.26633, 
-79.2585, -79.25765, -79.27642, -79.2754, -79.28565, -79.29443, 
-79.28264, -79.24703, -79.24706, -79.24706, -79.247, -79.24863, 
-79.30019, -79.3023, -79.30154, -79.24672, -79.24672, -79.24672, 
-79.2969, -79.25789, -79.26389, -79.26274, -79.26976, -79.29689, 
-79.26999, -79.29676, -79.29634, -79.29671, -79.29676, -79.29739, 
-79.30154, -79.28073, -79.28101, -79.28164, -79.28447, -79.28244, 
-79.28674, -79.28123, -79.28515, -79.29752, -79.29449, -79.30283, 
-79.26671, -79.29492, -79.24765, -79.23946, -79.24728, -79.24826, 
-79.2477), lat = c(25.69511, 25.69913, 25.69909, 25.69808, 25.69791, 
25.6987, 25.6976, 25.69887, 25.69887, 25.69812, 25.698905, 25.69876, 
25.69545, 25.69664, 25.69499, 25.69707, 25.69877, 25.69738, 25.69877, 
25.69786, 25.698363, 25.69876, 25.69884, 25.698363, 25.699164, 
25.69808, 25.69149, 25.69682, 25.69919, 25.69897, 25.6987, 25.69785, 
25.69937, 25.69917, 25.69759, 25.69846, 25.69897, 25.69231, 25.69592, 
25.69592, 25.699178, 25.73256, 25.73151, 25.73248, 25.7323, 25.73618, 
25.75886, 25.76018, 25.69645, 25.69439, 25.69597, 25.69989, 25.69769, 
25.73171, 25.73161, 25.73156, 25.23161, 25.73258, 25.69319, 25.6926, 
25.69706, 25.73245, 25.73245, 25.73245, 25.70004, 25.77368, 25.7471, 
25.72693, 25.72724, 25.6997, 25.72725, 25.69891, 25.69924, 25.69885, 
25.69891, 25.69776, 25.69712, 25.72679, 25.72655, 25.73001, 25.72484, 
25.72546, 25.72225, 25.72625, 25.72549, 25.69781, 25.69989, 25.69358, 
25.69701, 25.69955, 25.73175, 25.73411, 25.73237, 25.73462, 25.73189
)), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 
12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 40L, 41L, 42L, 43L, 
44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 
57L, 58L, 59L, 60L, 65L, 66L, 67L, 68L, 69L, 70L, 71L, 72L, 73L, 
74L, 75L, 76L, 77L, 78L, 79L, 80L, 81L, 82L, 83L, 84L, 85L, 86L, 
87L, 88L, 89L, 90L, 91L, 92L, 93L, 94L, 95L, 96L, 97L, 98L, 99L, 
100L, 135L, 136L, 137L, 138L, 139L, 140L, 141L, 142L, 143L, 144L, 
145L, 146L, 147L, 148L, 149L, 150L, 151L, 152L), class = "data.frame")

多边形定义 #

# POL1
coords1 <- matrix(c(-79.0877, 25.7213, 
                     -79.32764, 25.69931, 
                     -79.32665, 25.20054, 
                     -79.04101, 25.21296),
                   ncol = 2, byrow = T)
coords1 <- rbind(coords1, coords1[1, ]) # 闭合多边形
pol_1 <- sf::st_sfc(sf::st_polygon(list(coords1)), crs = 4326)

# POL2
coords2 <- matrix(c(-79.25322, 25.76284,
                      -79.21888, 25.76624, 
                      -79.21459, 25.72063,
                      -79.25613, 25.72349),
                    ncol = 2, byrow = T)
coords2 <- rbind(coords2, coords2[1, ]) # 闭合多边形
pol_2 <- sf::st_sfc(sf::st_polygon(list(coords2)), crs = 4326)

# POL3
coords3 <- matrix(c(-79.25446, 25.7759, 
                      -79.26378, 25.70252, 
                      -79.31184,25.71772, 
                      -79.26071, 25.77646), 
                    ncol = 2, byrow = T)
coords3 <- rbind(coords3, coords3[1, ]) # 闭合多边形
pol_3 <- sf::st_sfc(sf::st_polygon(list(coords3)), crs = 4326)

高效解决方案 #

核心思路是将三个多边形合并为一个带标识列的sf对象，然后通过空间连接一次性完成点与多边形的归属匹配，避免多次调用st_within和循环判断。

步骤1：合并多边形并添加标识 #

将三个多边形整合成一个sf dataframe，新增loc列标记每个多边形的标识（"One"、"Two"、"Three"）：

library(sf)

# 将单个多边形转为sf dataframe并添加标识
pol_df <- rbind(
  st_sf(loc = "One", geometry = pol_1),
  st_sf(loc = "Two", geometry = pol_2),
  st_sf(loc = "Three", geometry = pol_3)
)

步骤2：将观测数据转为sf对象 #

test_sf <- st_as_sf(testdata, coords = c("lon", "lat"), crs = 4326)

步骤3：空间连接匹配归属 #

使用st_join结合st_within谓词，一次性将每个点匹配到所属的多边形标识：

# 空间连接，保留所有观测点，匹配所属多边形
result <- st_join(test_sf, pol_df, join = st_within, left = TRUE)

# 移除几何列，转为普通dataframe（可选）
result_df <- st_drop_geometry(result)

说明 #

• st_join的join = st_within指定用"点在多边形内"的规则匹配
• left = TRUE确保所有观测点都被保留，不在任何多边形内的点loc列会是NA
• 若存在点同时属于多个多边形，st_join会返回多行匹配结果，可根据需求用dplyr::distinct或其他逻辑去重

对比原方法的优势 #

• 代码更简洁，无需多次调用st_within和手动循环判断
• 空间连接是sf库的优化操作，处理大规模数据时效率更高
• 扩展性更强，新增多边形只需在合并时添加对应行即可

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