23

I am migrating code from sp package to the newer sf package. My previous code I had a polygon SpatialDataFrame (censimentoMap) and a SpatialPointDataFrame (indirizzi.sp) and I got the polygon cell id ("Cell110") for each point laying within with the instruction below:

points.data <- over(indirizzi.sp, censimentoMap[,"Cell110"])

Actually I created two sf objects:

shape_sf <- st_read(dsn = shape_dsn) shape_sf <- st_transform(x=shape_sf, crs=crs_string) and indirizzi_sf = st_as_sf(df, coords = c("lng", "lat"), crs = crs_string)

And I am looking for the sf equivalent of the above instruction... Migth it be:

ids<-sapply(st_intersects(x=indirizzi_sf,y=shshape_sfpeCrif), function(z) if (length(z)==0) NA_integer_ else z[1]) cell_ids <- shape_sf[ids,"Cell110"]

2 Answers 2

24

You can get the same result by using st_join: First create a demo polygon and some points with sf.

library(sf)
library(magrittr)

poly <- st_as_sfc(c("POLYGON((0 0 , 0 1 , 1 1 , 1 0, 0 0))")) %>% 
  st_sf(ID = "poly1")    

pts <- st_as_sfc(c("POINT(0.5 0.5)",
                   "POINT(0.6 0.6)",
                   "POINT(3 3)")) %>%
  st_sf(ID = paste0("point", 1:3))

now see the result using over on sp objects

over(as(pts, "Spatial"), as(polys, "Spatial"))
>#      ID
># 1 poly1
># 2 poly1
># 3  <NA>

now equivalent with sf st_join

st_join(pts, poly, join = st_intersects)
># Simple feature collection with 3 features and 2 fields
># geometry type:  POINT
># dimension:      XY
># bbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
># epsg (SRID):    NA
># proj4string:    NA
>#     ID.x  ID.y               .
># 1 point1 poly1 POINT (0.5 0.5)
># 2 point2 poly1 POINT (0.6 0.6)
># 3 point3  <NA>     POINT (3 3)

or for the exact same result

as.data.frame(st_join(pts, poly, join = st_intersects))[2] %>% setNames("ID")

>#    ID
># 1 poly1
># 2 poly1
># 3  <NA>
6

Instead of using

st_join(pts, poly, join = st_intersects)

you can just use

st_intersection(pts, poly)

It is faster, and gives only points that are inside the poly.

Your Answer

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

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