9

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

2 Answers 2

7

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.)

3

I am sympathetic, this is a bit of a pain. You have to explode the slots, of the gBuffer polygon object, into individual polygons.

require(sp)
require(rgeos)
pts <- SpatialPoints(cbind(c(1, 1, 3), c(1, 2, 3))) 

b <- gBuffer(pts, width=0.6)

over(pts, b)

###########################################################
# explodes slots into individual polygons
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 <- SpatialPolygonsDataFrame(b.spolys, data.frame(FID=as.numeric(row.ids)) ) 
###########################################################

over(pts, b)
1
  • Have just tested and it seems promising, but my holes are no longer holes :(. Can we adapt this to keep holes together with the polys they're associated with? (I was heading along similar lines as your suggestion, with n <- length(buff@polygons[[1]]@Polygons); SpatialPolygons(mapply(function(x, y) Polygons(list(x), y), buff@polygons[[1]]@Polygons, seq_len(n))) ... I should have mentioned the hole problem in the post, sorry.
    – jbaums
    Commented Jan 21, 2015 at 20:09

Your Answer

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

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