I have shapefile with about 15,000 polygons. Some polygons overlap with other polygons. Most of them don't overlap with any other polygon. I also have a points data set with about 100 Million points (lat/lng).

I need to assign polygon(s) to each of the points. If I use over function I assign only one polygon to any give point. How do I assign two or more polygons to a point if that point is in multiple polygons.

I cannot run a loop polygon-by-polygon or point-by-point because there are too many polygons and too many points. However, processing with bunch of points/polygons in every loop, like a million points in one loop and/or 1000 polygons in one loop, is OK.


Example data (something you should provide)

n <- 100
x <- runif(n) * 360 - 180
y <- runif(n) * 180 - 90
xy <- cbind(x, y)
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
pols <- spPolygons(p1, p2, p3)

plot(pols, col=rainbow(3, alpha=.5))

Now do:

e <- extract(pols, xy)

To see (part of) the result:

#  point.ID poly.ID
#1        1       2
#2        2       1
#3        3      NA
#4        4      NA
#5        5      NA

To see points that overlap with multiple polygons:

tb <- table(e[,1])
i <- as.integer(names(tb[tb>1]))
points(xy[i,], pch=20, col='red')
| improve this answer | |
  • Thanks for the response. This is half of the solution I was looking for. This tells me which points belong to more than one polygon. But this does not give me what those polygons are. – Madhava Mahishi Jul 18 '16 at 4:40
  • It does give you that. See head(e) (I have added that) – Robert Hijmans Jul 19 '16 at 21:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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