[R-sig-Geo] merging polygons and IDs

[Robert J. Hijmans]

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.Bivand at 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 at 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.Bivand at nhh.no
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo