I would like to perform a spatial overlay to identify the polygon within which a set of points fall. However I first want to buffer and dissolve the polygons, such that points falling anywhere within merged polygons (but not within any holes) will be labelled similarly by the overlay procedure.
Unfortunately, the buffering and/or dissolving process I'm using is reducing the SpatialPolygons
object to a multi-polygon. Using gBuffer
creates a multi-polygon when byid=FALSE
, but fails to union overlapping polygons when byid=TRUE
. In the latter case, subsequently dissolving polygons with gUnaryUnion
again creates a multi-polygon. Overlaying SpatialPoints
with these multi-polygon results in all contained points being reported as falling within polygon 1.
Here's a related toy example with buffered points, where I'd like to do an overlay of the points with the buffered polygons:
library(sp)
library(rgeos)
pts <- SpatialPoints(cbind(c(1, 1, 2, 3), c(1, 2, 1.5, 2.5)))
plot(gBuffer(pts, width=0.6), lwd=2)
points(pts, pch=20, cex=2)
text(coordinates(pts), labels=seq_len(length(pts)), pos=4, font=2)
And the results of some overlays...
With
byid=FALSE
:b <- gBuffer(pts, width=0.6) over(pts, b) # [1] 1 1 1 1
With
byid=TRUE
:b2 <- gBuffer(pts, width=0.6, byid=TRUE) over(pts, b2) # [1] 1 2 3 4
With
byid=TRUE
and subsequentgUnaryUnion
:b3 <- gUnaryUnion(b2) over(pts, b3) # [1] 1 1 1 1
I am looking for the correct method to achieve the outcome 1 1 1 2
, i.e. points 1, 2 and 3 fall within polygon 1 (previously polygons 1, 2, and 3, which have been merged), and point 4 falls within polygon 2 (originally polygon 4).
EDIT
My true data include holes, and ideally I'd like to be able to conserve them. I've updated the example data above to include a hole after buffering.
Applying @JeffreyEvans's approach to these polys:
polys <- b@polygons[[1]]@Polygons
pl <- vector("list", length(polys))
for (i in 1:length(polys)) { pl[i] <- Polygons(list(polys[[i]]), i) }
b.spolys <- SpatialPolygons(pl)
row.ids <- sapply(slot(b.spolys, "polygons"), function(i) slot(i, "ID"))
b.exploded <- SpatialPolygonsDataFrame(b.spolys, data.frame(FID=as.numeric(row.ids)))
results in:
over(pts, b.exploded)
# FID
# 1 2
# 2 2
# 3 2
# 4 1
which is the expected outcome (I'm not too fussed about the poly order - 2,2,2,1
vs. 1,1,1,2
), but b.exploded
has lost the information about the hole, instead representing it as a non-hole polygon. It obviously doesn't affect the outcome for the toy example, but an overlay involving points located in the hole will be misleading, e.g.:
in_hole <- SpatialPoints(cbind(1.375, 1.5))
over(in_hole, b.exploded)
# FID
# 1 2