Определите, в каком многоугольнике лежит точка - PullRequest
Новая
       89

Определите, в каком многоугольнике лежит точка

0 голосов
/ 06 мая 2020

У меня есть два фрейма данных: один с точками (широта / долгота) и один с многоугольниками. Я хочу определить, в каком многоугольнике лежит точка. Есть ли способ сделать это, кроме написания функции points.in.polygon () из пакета sp?

Ответы [ 2 ]

2 голосов
/ 06 мая 2020
Пакет

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)
1 голос
/ 06 мая 2020

Я думаю, что более в библиотеке sp может быть тем, что вы ищете. Вот упрощенный пример, адаптированный со страницы документации, указанной выше. 

# 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')

Over Results

...