【问题标题】：Determine in which polygon a point lies确定点位于哪个多边形中
【发布时间】：2020-05-06 04:54:06
【问题描述】：

我有两个数据框：一个带有点（纬度/经度），一个带有多边形。我想确定一个点位于哪个多边形中。除了从 sp-package 编写 points.in.polygon() 的函数之外，还有其他方法吗？

【问题讨论】：

【解决方案1】：

sf 包为您提供了大部分空间操作。 sf 包替换 sp 作为 R 中此类操作的黄金标准。

在这里，您正在寻找具有 within 匹配的空间连接。

数据

library(sf)
nc <- st_transform(st_read(system.file("shape/nc.shp", package="sf")), 2264)                
gr = st_sf(
  label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "),
  geom = st_make_grid(nc))
gr <- gr %>% st_centroid() %>% st_as_sf()

gr 是点对象

gr
Simple feature collection with 100 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 538593.4 ymin: 98164 xmax: 2920556 ymax: 994368.4
epsg (SRID):    2264
proj4string:    +proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
First 10 features:
   label                   geom
1   J  1 POINT (538593.4 98164)
2   J  2 POINT (803255.9 98164)
3   J  3  POINT (1067918 98164)
4   J  4  POINT (1332581 98164)
5   J  5  POINT (1597243 98164)
6   J  6  POINT (1861906 98164)
7   J  7  POINT (2126568 98164)
8   J  8  POINT (2391231 98164)
9   J  9  POINT (2655893 98164)
10  J 10  POINT (2920556 98164)

多边形中的点

简单

points_in_poly <- gr %>% st_join(nc, join = st_within, left = FALSE)

产生：

points_in_poly
Simple feature collection with 55 features and 15 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 538593.4 ymin: 98164 xmax: 2920556 ymax: 994368.4
epsg (SRID):    2264
proj4string:    +proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
First 10 features:
   label  AREA PERIMETER CNTY_ CNTY_ID      NAME  FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79
7   J  7 0.212     2.024  2241    2241 Brunswick 37019  37019       10  2181     5     659  2655     6     841
17  I  7 0.240     2.365  2232    2232  Columbus 37047  37047       24  3350    15    1431  4144    17    1832
27  H  7 0.225     2.107  2162    2162    Bladen 37017  37017        9  1782     8     818  2052     5    1023
28  H  8 0.214     2.152  2185    2185    Pender 37141  37141       71  1228     4     580  1602     3     763
35  G  5 0.163     1.716  2095    2095     Union 37179  37179       90  3915     4    1034  5273     9    1348
36  G  6 0.080     1.188  2123    2123  Scotland 37165  37165       83  2255     8    1206  2617    16    1436
37  G  7 0.225     2.107  2162    2162    Bladen 37017  37017        9  1782     8     818  2052     5    1023
38  G  8 0.204     1.871  2100    2100    Duplin 37061  37061       31  2483     4    1061  2777     7    1227
39  G  9 0.125     2.868  2156    2156  Carteret 37031  37031       16  2414     5     341  3339     4     487
41  F  1 0.051     1.096  2109    2109      Clay 37043  37043       22   284     0       1   419     0       5
                        geom
7      POINT (2126568 98164)
17  POINT (2126568 197742.3)
27  POINT (2126568 297320.5)
28  POINT (2391231 297320.5)
35  POINT (1597243 396898.8)
36  POINT (1861906 396898.8)
37  POINT (2126568 396898.8)
38  POINT (2391231 396898.8)
39  POINT (2655893 396898.8)
41 POINT (538593.4 496477.1)

【讨论】：

【解决方案2】：

我认为sp 库中的over 可能是您正在寻找的。这是改编自上面链接的文档页面的简化示例。

# example polygon dataframe
r1 <- cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409,180162, 180114), 
            c(332349, 332057, 332342, 333250, 333558, 333676, 332618, 332413, 332349))
r2 <- cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 179524, 179979, 180042),
            c(332373, 332026, 331426, 330889, 330683, 331133, 331623, 332152, 332357, 332373))
r3 <- cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
            c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 329783, 329665, 329720, 329933, 330478, 331062, 331086))
r4 <- cbind(c(180304, 180403,179632,179420,180304),
            c(332791, 333204, 333635, 333058, 332791))
sr1 <- Polygons(list(Polygon(r1)),"r1")
sr2 <- Polygons(list(Polygon(r2)),"r2")
sr3 <- Polygons(list(Polygon(r3)),"r3")
sr4 <- Polygons(list(Polygon(r4)),"r4")
sr <- SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf <- SpatialPolygonsDataFrame(sr, data.frame(poly=c("r1","r2","r3","r4"),
                                                row.names=c("r1","r2","r3","r4")))

# load some points to test
data(meuse)
points <- meuse[,c('x','y')]
coordinates(points) = ~x+y

# test which polygon each point is over
result <- over(points, srdf)

# check the results
plot(points, col=result$poly)
legend("topleft", legend=levels(result$poly), pch=16, 
       col=c('black','red','green','blue'))
polygon(r1, border='black')
polygon(r2, border='red')
polygon(r3, border='green')
polygon(r4, border='blue')

【讨论】：