David, I think you could also do this:
library(raster) a <- disaggregate(aggregate(pols)) ids <- over(pts, a) Now you have the aggregated polygons and for each point you know which polygon it belongs to. Perhaps continue with something like this: tab <- table(ids, as.integer(names(ids))) tab <- t(t(tab) * as.integer(colnames(tab))) x <- apply(tab, 1, function(i) paste(i[i!=0], collapse='_')) p <- SpatialPolygonsDataFrame(a, data.frame(points=x)) head(p) Robert On Sun, Jan 11, 2015 at 11:57 AM, Roger Bivand <roger.biv...@nhh.no> wrote: > On Sun, 11 Jan 2015, David Depew wrote: > >> Dear list, >> >> given a set of polygons created from points with a specified buffer >> distance, >> >> library(sp) >> library(rgeos) >> box <- readWKT("POLYGON((0 0, 0 1000, 1000 1000, 1000 0, 0 0))") >> plot(box) >> set.seed(1) >> pts <- spsample(box, n=40, type="random") >> pols <- gBuffer(pts, byid=TRUE, width=50) >> plot(pols, add=TRUE) >> >> >> I'd like to merge polygons where the buffers overlap >> i.e. >> merg.poly=gUnaryUnion(pols) >> plot(merg.poly) >> >> but I'm having difficulty figuring out how to specify that the merged >> polygons retain some unique Id. Clearly in the package help, one would >> specify this by use of the "id" argument, however since the polygons >> started as points, overlapping points do not have the same ID. Suppose one >> approach is to define a common ID for overlapping polys before using >> gUnaryUnion, but I'm at a loss where to begin with that approach. Any >> advice would be greatly appreciated, > > > Thanks for a self-contained example! > > Right, this isn't obvious, as you need to go to a list/graph representation > of the intersects predicate, which takes us outside rgeos. We'll compare > this with the equivalent for distance neighbours in spdep: > > library(spdep) > dnb <- dnearneigh(pts, 0, 100) # 2*50 buffer > nc <- n.comp.nb(dnb) # find graph components > merg.poly <- gUnaryUnion(pols, nc$comp.id) > plot(merg.poly) > > which jumps over the buffered objects. To use the buffered objects, do: > > gI <- gIntersects(pols, pols, byid=TRUE, returnDense=FALSE) # list > for (i in seq(along=gI)) {gI[[i]] <- gI[[i]][-(which(gI[[i]] == i))]; > if (length(gI[[i]]) == 0) gI[[i]] <- 0L} # remove i<->i edges > class(gI) <- "nb" > attr(gI, "region.id") <- names(gI) > nc1 <- n.comp.nb(gI) > all.equal(nc, nc1) # same result > > So what this is doing is having gIntersects return a list. Since i > intersects i, we have to drop the i<->i edges, pretend that this is an nb > object, and use the n.comp.nb() function to find the graph components. These > two approaches give the same results. There is a new vignette in spdep > starting to look at representations of neighbours as graph objects: > > http://cran.r-project.org/web/packages/spdep/vignettes/nb_igraph.html > > Hope this helps, > > Roger > >> >> Thanks, >> >> [[alternative HTML version deleted]] >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > -- > Roger Bivand > Department of Economics, Norwegian School of Economics, > Helleveien 30, N-5045 Bergen, Norway. > voice: +47 55 95 93 55; fax +47 55 95 91 00 > e-mail: roger.biv...@nhh.no > > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo