0

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 

3 Answers 3

1

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: plot output

0
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")

enter image description here

2
  • 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?
    – Alwin
    May 12, 2018 at 10:35
  • @Alwin next to xy? Can you make a smaller reproducible example and show (perhaps color by hand?) what polygons you're after? May 13, 2018 at 10:09
1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.