Point in overlapping polygons using R

Question

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.

Answer 1

Example data (something you should provide)

library(raster)
n <- 100
set.seed(0)
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))
points(xy)

Now do:

e <- extract(pols, xy)

To see (part of) the result:

head(e)
#  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')