I'm working with a dataframe containing longitude and latitude for each point. I have a shapefile containing mutually exclusive polygons. I would like to find the index of the polygon it where each point is contained. Is there a specific function that helps me achieve this? I've been trying with the sf package, but I'm open to doing it with another one. Any help is greatly appreciated.
Asked
Active
Viewed 56 times
-1
-
1It's easier to help you if you include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input and desired output that can be used to test and verify possible solutions. – MrFlick Jun 29 '20 at 22:59
-
1Your third sentence, which seems crucial, is very unclear. – Edward Jun 29 '20 at 23:02
-
1You might want to check the function sf::st_contains(). – Henry Cyranka Jun 29 '20 at 23:41
-
Another option: `sp::point.in.polygon` – r2evans Jun 30 '20 at 02:08
-
Try looking at this :https://stackoverflow.com/questions/3647744/intersecting-points-and-polygons-in-r – Shubham Pujan Jun 30 '20 at 07:08
1 Answers
1
I believe you may be looking for function sf::st_intersects() - in combination with sparse = TRUE setting it returns a list, which can be in this use case (points & a set of non-overlapping polygons) converted to a vector easily.
Consider this example, built on the North Carolina shapefile shipped with {sf}
library(sf)
# as shapefile included with sf package
shape <- st_read(system.file("shape/nc.shp", package="sf")) %>%
st_transform(4326) # WGS84 is a good default
# three semi random cities
cities <- data.frame(name = c("Raleigh", "Greensboro", "Wilmington"),
x = c(-78.633333, -79.819444, -77.912222),
y = c(35.766667, 36.08, 34.223333)) %>%
st_as_sf(coords = c("x", "y"), crs = 4326) # again, WGS84
# plot cities on full map
plot(st_geometry(shape))
plot(cities, add = T, pch = 4)
# this is your index
index_of_intersection <- st_intersects(cities, shape, sparse = T) %>%
as.numeric()
# plot on subsetted map to doublecheck
plot(st_geometry(shape[index_of_intersection, ]))
plot(cities, add = T, pch = 4)
Jindra Lacko
- 2,670
- 2
- 12
- 26