If you use the
sf package, then
st_intersection can get you all the intersections, with none of the non-intersections...
Sample data - some points and two polygons which overlap. The polygons have a "Region" attribute which is A or B:
Then if I do:
inters = st_intersection(points,polys)
Name Region geom
4 D A POINT (-15.73582 13.59678)
5 R A POINT (-15.68296 13.56941)
6 H A POINT (-15.63011 13.57446)
7 KK A POINT (-15.58938 13.56983)
8 IIJ A POINT (-15.53783 13.5972)
10 Y A POINT (-15.46808 13.56646)
10.1 Y B POINT (-15.46808 13.56646)
11 ksjdkj B POINT (-15.46981 13.64015)
12 kdjihj B POINT (-15.39616 13.61362)
14 ksjicc B POINT (-15.27443 13.62752)
And the columns are the attributes from each overlapping point-polygon combination -
geometry comes from
Region comes from
The points that aren't in polygons do not appear in this result.
If speed is an issue, your problem is trivially parallel over polygons, so you can divide your polygons into N sets and run on N processors for a potential N-fold speed increase.
For example, suppose you have 1200 polygons in
polys. Suppose computing
st_intersection(points, polys) takes 10 seconds. Computing the intersection with half the polygons with
st_intersection(points, polys[1:600,]) should take five seconds. If you can do that at the same time as
st_intersection(points, polys[601:1200,]) then you can do all 1200 polygon tests in five seconds - and then a short time to join the two results.
If you have an 8-core processor then you could divide your polygons into 8 groups, and get an 8-fold speedup, using R's various parallel packages.
If you have access to some cloud resources, like Amazon EC2 or Google Cloud or Azure, you can "spin up" 600 virtual machines, have each one test 2 polygons, and get a potential massive speedup.