19

I am working with shapefiles in R, one is point.shp the other is a polygon.shp. Now, I would like to intersect the points with the polygon, meaning that all the values from the polygon should be attached to the table of the point.shp.

I tried overlay() and spRbind in package sp, but nothing did what I expected them to do.

Could anyone give me a hint?

3 Answers 3

25

With the new sf package this is now fast and easy:

library(sf)
out <- st_intersection(points, poly)

Additional options

If you do not want all fields from the polygon added to the point feature, just call dplyr::select() on the polygon feature before:

library(magrittr)
library(dplyr)
library(sf)

poly %>% 
  select(column-name1, column-name2, etc.) -> poly

out <- st_intersection(points, poly)

If you encounter issues, make sure that your polygon is valid:

st_is_valid(poly)

If you see some FALSE outputs here, try to make it valid:

poly <- st_make_valid(poly) 

Note that these 'valid' functions depend on a sf installation compiled with liblwgeom.

3
  • 1
    Unfortunately this didn't work for my polgon of type SpatialPolygonsDataFrame and SpatialPointsDataFrame. I get the error message: no applicable method for 'st_intersection' applied to an object of class "c('SpatialPointsDataFrame', 'SpatialPoints', 'Spatial', 'SpatialVector')
    – yPennylane
    Commented Jun 28, 2018 at 14:55
  • 2
    @yPennylane You need to do this with sf objects. The classes you mention are from package sp. Read in your data using the sf package, i.e. using the st_read() function.
    – pat-s
    Commented Jun 29, 2018 at 9:46
  • @Jens Maybe it's time to switch the accepted answer? The old one seems deprecated.
    – pat-s
    Commented Jul 25, 2022 at 17:25
15

If you do overlay(pts, polys) where pts is a SpatialPointsDataFrame object and polys is a SpatialPolygonsDataFrame object then you get back a vector the same length as the points giving the row of the polygons data frame. So all you then need to do to combine the polygon data onto the points data frame is:

 o = overlay(pts, polys)
 pts@data = cbind(pts@data, polys[o,])

HOWEVER! If any of your points fall outside all your polygons, then overlay returns an NA, which will cause polys[o,] to fail, so either make sure all your points are inside polygons or you'll have to think of another way to assign values for points outside the polygon...

2
  • Thank you ! Somehow I overlooked the second line! Works perfectly well. but, I will post another question in a second: what to do if one wants to attach a simple data.frame to a spatialpolygondataframe?
    – Jens
    Commented Sep 6, 2010 at 10:09
  • 8
    @spacedman this answer is deprecated! The overlay function is now over. Its not clear what the polys object in your second line is referring to.
    – colin
    Commented Mar 4, 2016 at 21:02
4

You do this in one line with point.in.poly fom spatialEco package.

library(spatialEco)

new_shape <- point.in.poly(pts, polys)

from the documentation: point.in.poly "intersects point and polygon feature classes and adds polygon attributes to points".

1
  • 2
    Note that spatialEco::point.in.poly is just a wrapper of sp::over. Commented Sep 1, 2017 at 13:07

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.