[GIS] buffer and union in R without reducing to a multi-polygon

rsp

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)

enter image description here

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 subsequent gUnaryUnion:

    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

Best Answer

It turns out that sp::disaggregate can be used to pull apart the singlepart polygons.

pts <- SpatialPoints(cbind(c(1, 1, 2, 3), c(1, 2, 1.5, 2.5))) 
b <- gBuffer(pts, width=0.6) 
over(pts, disaggregate(b))

# [1] 1 1 1 2

in_hole <- SpatialPoints(cbind(1.375, 1.5))
over(in_hole, disaggregate(b))

# [1] NA

(raster::aggregate can be used to combine them again.)

Related Question