Determine in which polygon a point lies

Question

I have two dataframes: one with points (lat/lon) and one with polygons. I want to determine in which polygon a point lies. Is there a way do to this, besides writing a function of the points.in.polygon() from the sp-package?

Answer 1

sf package gives you most of the spatial operations. sf package replaced sp as a gold standard for such operations in R.

Here, you are looking for a spatial join with a within matching.

Data

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 is point object

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 polygon

Simply

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

which yields:

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)

Answer 2

I think over in the sp library might be what you are looking for. Here's a simplified example adapted from the documentation page linked above.

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