How Can I Do a Spatial Join with the Sf Package Using St_Join()
Here's a Toy Example I've Been Wrestling with # Make Points Point1 # Determine Which Points Fall Inside Which Polygons > St_Join(Pts, Polys, Join =...
Here's a toy example I've been wrestling with
# Make points
point1 <- c(.5, .5)
point2 <- c(.6, .6)
point3 <- c(3, 3)
mpt <- st_multipoint(rbind(point1, point2, point3)) # create multipoint
# Make polygons
square1 <- rbind(c(0, 0), c(1, 0), c(1,1), c(0, 1), c(0, 0))
square2 <- rbind(c(0, 0), c(2, 0), c(2,2), c(0, 2), c(0, 0))
square3 <- rbind(c(0, 0), c(-1, 0), c(-1,-1), c(0, -1), c(0, 0))
mpol <- st_multipolygon(list(list(square1), list(square2), list(square2))) # create multipolygon
# Convert to class 'sf'
pts <- st_sf(st_sfc(mpt))
polys <- st_sf(st_sfc(mpol))
# Determine which points fall inside which polygons
st_join(pts, polys, join = st_contains)
The last line produces
Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) :
cannot coerce class "c("sfc_MULTIPOINT", "sfc")" to a data.frame
How can I do a spatial join to determine which points fall inside which polygons?
3 Answers
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:
> pts
Simple feature collection with 1 feature and 0 fields
geometry type: MULTIPOINT
dimension: XY
bbox: xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID): NA
proj4string: NA
st_sfc.mpt.
1 MULTIPOINT(0.5 0.5, 0.6 0.6...
> polys
Simple feature collection with 1 feature and 0 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 0 ymin: 0 xmax: 2 ymax: 2
epsg (SRID): NA
proj4string: 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))
> polys
Simple feature collection with 3 features and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: -1 ymin: -1 xmax: 2 ymax: 2
epsg (SRID): NA
proj4string: 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, ...
> pts
Simple feature collection with 3 features and 1 field
geometry type: POINT
dimension: XY
bbox: xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID): NA
proj4string: 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 polygons
pi <- 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 point3
1 poly1 TRUE TRUE FALSE
2 poly2 TRUE TRUE FALSE
3 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 fields
geometry type: MULTIPOINT
dimension: XY
bbox: xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID): NA
proj4string: NA
geometry
1 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
Marcus Vance is a cybersecurity auditor and technology writer dedicated to educating the public about online safety, data privacy regulations, enterprise security, and emerging cyber threats.
