Points in multiple polygons using R

Question

Currently I have two data.frames, one of polygons (poly.x, poly.y, enum) and one of points (pt.x, pt.y) where enum is the id of the polygon. I am trying to determine which points belong to which polygons so I get a data.frame of (pt.x, pt.y, enum).

My first attempt uses point.in.polygon from the sp package and lapply functions to find which polygon(s) the point belongs to. While my code works, it takes a long time on large data sets.

My second attempt uses over also from the sp package, cobbled together from questions on gis stackexchange. While it is much faster, I cannot seem to get the correct output from over as it is a dataframe of 1s and NAs.

Below I've included a simplified working example (npoly can be changed to test the speed of different methods) as well as my working attempt using sp::point.in.polygon and nonsensical output from my sp::over attempt. I'm not fussed which method I end up using as long as it's fast.

Any help would be much appreciated!

#-------------------------------------------
# Libraries
library(ggplot2)  # sample plots
library(dplyr)    # bind_rows(), etc
library(sp)       # spatial data
# Sample data
npoly = 100
#  polygons
localpolydf <- data.frame(
  x = rep(c(0, 1, 1, 0), npoly) + rep(0:(npoly-1), each = 4),
  y = rep(c(0, 0, 1, 1), npoly),
  enum = rep(1:npoly, each = 4))
#  points
offsetdf <- data.frame(
  x = seq(min(localpolydf$x) - 0.5, max(localpolydf$x) + 0.5, by = 0.5),
  y = runif(npoly*2 + 3, 0, 1))
# Sample plot
ggplot() + 
  geom_polygon(aes(x, y, group = enum), 
               localpolydf, fill = NA, colour = "black") +
  geom_point(aes(x, y), offsetdf)

#-------------------------------------------
# Dplyr and lapply solution for point.in.polygon
ptm <- proc.time() # Start timer
#  create lists
offsetlist <- split(offsetdf, rownames(offsetdf))
polygonlist <- split(localpolydf, localpolydf$enum)
# lapply over each pt in offsetlist
pts <- lapply(offsetlist, function(pt) {
  # lapply over each polygon in polygonlist
  ptpoly <- lapply(polygonlist, function(poly) {
    data.frame(
      enum = poly$enum[1],
      ptin = point.in.polygon(pt[1,1], pt[1,2], poly$x, poly$y))
  })
  ptpoly <- bind_rows(ptpoly) %>% filter(ptin != 0)
  if (nrow(ptpoly) == 0) return(data.frame(x = pt$x, y = pt$y, enum = NA, ptin = NA))
  ptpoly$x = pt$x
  ptpoly$y = pt$y
  return(ptpoly[c("x", "y", "enum", "ptin")])
})
pts_apply <- bind_rows(pts)
proc.time() - ptm # end timer

#-------------------------------------------
# Attempted sp solution for over
ptm <- proc.time() # Start timer
# Split the dataframe into a list based on enum and then remove enum from df in the list
polygonlist <- split(localpolydf, localpolydf$enum)
polygonlist <- lapply(polygonlist, function(x) x[,c("x", "y")])
# Convert the list to Polygon, then create a Polygons object
polygonsp <- sapply(polygonlist, Polygon)
polygonsp <- Polygons(polygonsp, ID = 1)
polygonsp <- SpatialPolygons(list(polygonsp))
plot(polygonsp)
# Convert points to coordinates
offsetps <- offsetdf
coordinates(offsetps) <- ~x+y
points(offsetps$x, offsetps$y)
# Determine polygons points are in
pts_sp <- over(offsetps, polygonsp)
proc.time() - ptm # end timer

#===========================================
# Output
#  Apply: point.in.polygon
> head(pts_apply)
   x   y         enum ptin
1 -0.5 0.2218138   NA   NA
2  4.0 0.9785541    4    2
3  4.0 0.9785541    5    2
4 49.0 0.3971479   49    2
5 49.0 0.3971479   50    2
6 49.5 0.1177206   50    1
user  system elapsed 
4.434   0.002   4.435 
# SP: over
> head(pts_sp)
1  2  3  4  5  6 
NA  1  1 NA  1 NA 
user  system elapsed 
0.048   0.000   0.047

Chris · Accepted Answer · 2018-05-28 00:08:30Z

An alternative to using over is to use sf::intersection as the sf package is becoming more and more popular.

Getting the data into sf objects took me a little bit of work but if you are working with external data you can just read in with st_read and it will already be in the correct form.

Here is how to approach:

library(tidyverse)
library(sf)

# convert into st_polygon friendly format (all polygons must be closed)
# must be a nicer way to do this!
localpoly <- localpolydf %>% split(localpolydf$enum) %>% 
lapply(function(x) rbind(x,x[1,])) %>%
lapply(function(x) x[,1:2]) %>%
lapply(function(x) list(as.matrix(x))) %>%
lapply(function(x) st_polygon(x))

# convert points into sf object
points <- st_as_sf(offsetdf,coords=c('x','y'),remove = F)

#convert polygons to sf object and add id column
polys <- localpoly %>% st_sfc() %>% st_sf(geom=.) %>% 
mutate(id=factor(1:100)) 

#find intersection
joined <- polys  %>% st_intersection(points) 


# Sample plot
ggplot() + geom_sf(data=polys) +
geom_sf(data=joined %>% filter(id %in% c(1:10)),aes(col=id)) +
lims(x=c(0,10))

Note that to use geom_sf at the time of writing you will need to install the development version of ggplot.

plot output:

Roman Luštrik · Accepted Answer · 2018-05-12 07:33:16Z

1

over returns an index of points inside a geometry. Perhaps something like this:

xy <- offsetps[names(na.omit(pts_sp == 1)), ]

plot(polygonsp, axes = 1, xlim = c(0, 10))
points(offsetps)
points(xy, col = "red")

answered May 12, 2018 at 7:33

Roman Luštrik

70.9k25 gold badges160 silver badges203 bronze badges

2 Comments

Alwin Over a year ago

Thanks for the heads up! However, I was mainly after a column name of polygon IDs (enum) next to xy, ie for any (x, y) which polygon (in this case a square) is it in?

Roman Luštrik Over a year ago

@Alwin next to xy? Can you make a smaller reproducible example and show (perhaps color by hand?) what polygons you're after?

Alwin · Accepted Answer · 2018-05-12 12:56:25Z

After having another look, I realised Roman did pts_sp == 1 because I only had 1 ID for all of my squares, i.e. when I did ID = 1.

Once I fixed that, I was able to a column with ID = enum. To handle points in multiple polygons I can use returnList = TRUE and add additional lines to convert the list to a data.frame but it isn't necessar here.

# Attempted sp solution
ptm <- proc.time() # Start timer
# Split the dataframe into a list based on enum and then remove enum from df in the list
polygonlist <- split(localpolydf, localpolydf$enum)
# Convert the list to Polygon, then create a Polygons object
polygonsp <- sapply(polygonlist, function(poly){
  Polygons(list(Polygon(poly[, c("x", "y")])), ID = poly[1, "enum"])
})
# polygonsp <- Polygons(polygonsp, ID = 1)
polygonsp <- SpatialPolygons(polygonsp)
plot(polygonsp)
# Convert points to coordinates
offsetps <- offsetdf
coordinates(offsetps) <- ~x+y
points(offsetps$x, offsetps$y)
# Determine polygons points are in
pts_sp <- over(offsetps, polygonsp)
pts_sp <- data.frame(
  x = offsetps$x, y = offsetps$y,
  enum = unique(localpolydf$enum)[pts_sp])
proc.time() - ptm # end timer

Collectives™ on Stack Overflow

Points in multiple polygons using R

3 Answers 3

Comments

2 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

Comments

2 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related