How can I do a spatial join with the sf package using st_join()

I'm also working my way around the features of the sf package, so apologies if this is not correct or there are better ways. I think one problem here is that if building the geometries like in your example you are not obtaining what you think:

> ptsSimple feature collection with 1 feature and 0 fieldsgeometry type:  MULTIPOINTdimension:      XYbbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3epsg (SRID):    NAproj4string:    NA                     st_sfc.mpt.1 MULTIPOINT(0.5 0.5, 0.6 0.6...> polysSimple feature collection with 1 feature and 0 fieldsgeometry type:  MULTIPOLYGONdimension:      XYbbox:           xmin: 0 ymin: 0 xmax: 2 ymax: 2epsg (SRID):    NAproj4string:    NA                    st_sfc.mpol.1 MULTIPOLYGON(((0 0, 1 0, 1 ...

You can see that you have only one "feature" both in pts and in polys. This means that you are building one "multipolygon" feature (that is, a polygon constituted by 3 parts), instead thatn three different polygons. The same goes for the points.

After digging a bit, I found this different (and in my opinion easier) way to build the geometries, using WKT notation:

polys <- st_as_sfc(c("POLYGON((0 0 , 0 1 , 1 1 , 1 0, 0 0))",                     "POLYGON((0 0 , 0 2 , 2 2 , 2 0, 0 0 ))",                      "POLYGON((0 0 , 0 -1 , -1 -1 , -1 0, 0 0))")) %>%   st_sf(ID = paste0("poly", 1:3))    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))> polysSimple feature collection with 3 features and 1 fieldgeometry type:  POLYGONdimension:      XYbbox:           xmin: -1 ymin: -1 xmax: 2 ymax: 2epsg (SRID):    NAproj4string:    NA     ID                              .1 poly1 POLYGON((0 0, 0 1, 1 1, 1 0...2 poly2 POLYGON((0 0, 0 2, 2 2, 2 0...3 poly3 POLYGON((0 0, 0 -1, -1 -1, ...> ptsSimple feature collection with 3 features and 1 fieldgeometry type:  POINTdimension:      XYbbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3epsg (SRID):    NAproj4string:    NA      ID              .1 point1 POINT(0.5 0.5)2 point2 POINT(0.6 0.6)3 point3     POINT(3 3)

you can see that now both polys and pts have three features.

We can now find the "intersection matrix" using:

# Determine which points fall inside which polygonspi <- st_contains(polys,pts, sparse = F) %>%   as.data.frame() %>%   mutate(polys = polys$ID) %>%   select(dim(pi)[2],1:dim(pi)[1])colnames(pi)[2:dim(pi)[2]] = levels(pts$ID)> pi  polys point1 point2 point31 poly1   TRUE   TRUE  FALSE2 poly2   TRUE   TRUE  FALSE3 poly3  FALSE  FALSE  FALSE

meaning (as pointed out @symbolixau in the comments) that polygons 1 and 2 contain points 1 and 2, while polygon 3 doesn't contain any points. Point 3 is instead not contained in any polygon.

HTH.

I see a different output:

> # Determine which points fall inside which polygons> st_join(pts, polys, join = st_contains)Simple feature collection with 1 feature and 0 fieldsgeometry type:  MULTIPOINTdimension:      XYbbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3epsg (SRID):    NAproj4string:    NA                        geometry1 MULTIPOINT(0.5 0.5, 0.6 0.6...

was this with the most recent CRAN version of sf?

Note, the original set of multipoint and multipolygon can be 'cast' to point and polygon, without creating new objects:

st_contains(polys %>% st_cast("POLYGON"), pts %>% st_cast("POINT"), sparse = F)#      [,1]  [,2]  [,3]#[1,]  TRUE  TRUE FALSE#[2,]  TRUE  TRUE FALSE#[3,] FALSE FALSE FALSE