Dear R-helpers,
thanks for yesterday's speeding-up tip. Here is my next query:
I have lots of polygons (not necessarily convex ones, and they never
have holes) given by x,y coordinates.
I want to get the polygon that is the union of these polygons. This is
my current method, but I am hoping there is a faster method (up to
thousands of polygons, each with ca. 40 xy points).
Example:
library(gpclib)
# A polygon
leaf <- structure(c(0, 1, 12.9, 16.5, 18.8, 17, 16.8, 15.5, 12.1, 8.2,
6.3, 5, 2, 0, -1.5, -4.3, -6.6, -10.3, -14.8, -19.4, -22.2, -23.5,
-22.2, -17.6, -7.8, 0, 0, -2.4, 2.8, 8.9, 19.9, 33.9, 34.8, 40.4,
49.7, 69.2, 77.4, 83.4, 91.4, 99, 92.8, 87.3, 81.2, 71.1, 57.6,
45.4, 39.2, 26, 15.6, 5.3, 0.6, 0), .Dim = c(26L, 2L), .Dimnames = list(
NULL, c("X", "Y")))
# Lots of polygons:
releaf <- function(leaf)cbind(leaf[,1]+rnorm(1,0,50),leaf[,2]+rnorm(1,0,50))
leaves <- replicate(500, releaf(leaf), simplify=FALSE)
# Make into gpc.poly class:
leaves <- lapply(leaves, as, "gpc.poly")
# Make union .....
system.time({
leavesoutline <- union(leaves[[1]], leaves[[2]])
for(i in 3:length(leaves))leavesoutline <- union(leavesoutline, leaves[[i]])
})
# about 1sec here.
# Check it:
plot(leavesoutline)
thanks!
Remko
-------------------------------------------------
Remko Duursma
Research Lecturer
Centre for Plants and the Environment
University of Western Sydney
Hawkesbury Campus
Richmond NSW 2753
Dept of Biological Science
Macquarie University
North Ryde NSW 2109
Australia
Mobile: +61 (0)422 096908
www.remkoduursma.com
Hi, I think you could use a concave hull from the alphahull package, http://yihui.name/en/2010/04/alphahull-an-r-package-for-alpha-convex-hull/ It may be difficult to find the right parameters if the polygons differ widely in edge lengths, though. HTH, baptiste On 2 June 2010 03:53, Remko Duursma <remkoduursma at gmail.com> wrote:> Dear R-helpers, > > thanks for yesterday's speeding-up tip. Here is my next query: > > I have lots of polygons (not necessarily convex ones, and they never > have holes) given by x,y coordinates. > > I want to get the polygon that is the union of these polygons. This is > my current method, but I am hoping there is a faster method (up to > thousands of polygons, each with ca. 40 xy points). > > Example: > > library(gpclib) > > # A polygon > leaf <- structure(c(0, 1, 12.9, 16.5, 18.8, 17, 16.8, 15.5, 12.1, 8.2, > 6.3, 5, 2, 0, -1.5, -4.3, -6.6, -10.3, -14.8, -19.4, -22.2, -23.5, > -22.2, -17.6, -7.8, 0, 0, -2.4, 2.8, 8.9, 19.9, 33.9, 34.8, 40.4, > 49.7, 69.2, 77.4, 83.4, 91.4, 99, 92.8, 87.3, 81.2, 71.1, 57.6, > 45.4, 39.2, 26, 15.6, 5.3, 0.6, 0), .Dim = c(26L, 2L), .Dimnames = list( > ? ?NULL, c("X", "Y"))) > > # Lots of polygons: > releaf <- function(leaf)cbind(leaf[,1]+rnorm(1,0,50),leaf[,2]+rnorm(1,0,50)) > leaves <- replicate(500, releaf(leaf), simplify=FALSE) > > # Make into gpc.poly class: > leaves <- lapply(leaves, as, "gpc.poly") > > # Make union ..... > system.time({ > leavesoutline <- union(leaves[[1]], leaves[[2]]) > for(i in 3:length(leaves))leavesoutline <- union(leavesoutline, leaves[[i]]) > }) > # about 1sec here. > > # Check it: > plot(leavesoutline) > > > > thanks! > > Remko > > > ------------------------------------------------- > Remko Duursma > Research Lecturer > > Centre for Plants and the Environment > University of Western Sydney > Hawkesbury Campus > Richmond NSW 2753 > > Dept of Biological Science > Macquarie University > North Ryde NSW 2109 > Australia > > Mobile: +61 (0)422 096908 > www.remkoduursma.com > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >
Reduce might work. Not sure about the speed advantages though. It does
simplify code.
Unionall <- function(x) Reduce('union', x)
leaveout <- Unionall(leaves)
On Tue, Jun 1, 2010 at 9:53 PM, Remko Duursma
<remkoduursma@gmail.com>wrote:
> Dear R-helpers,
>
> thanks for yesterday's speeding-up tip. Here is my next query:
>
> I have lots of polygons (not necessarily convex ones, and they never
> have holes) given by x,y coordinates.
>
> I want to get the polygon that is the union of these polygons. This is
> my current method, but I am hoping there is a faster method (up to
> thousands of polygons, each with ca. 40 xy points).
>
> Example:
>
> library(gpclib)
>
> # A polygon
> leaf <- structure(c(0, 1, 12.9, 16.5, 18.8, 17, 16.8, 15.5, 12.1, 8.2,
> 6.3, 5, 2, 0, -1.5, -4.3, -6.6, -10.3, -14.8, -19.4, -22.2, -23.5,
> -22.2, -17.6, -7.8, 0, 0, -2.4, 2.8, 8.9, 19.9, 33.9, 34.8, 40.4,
> 49.7, 69.2, 77.4, 83.4, 91.4, 99, 92.8, 87.3, 81.2, 71.1, 57.6,
> 45.4, 39.2, 26, 15.6, 5.3, 0.6, 0), .Dim = c(26L, 2L), .Dimnames = list(
> NULL, c("X", "Y")))
>
> # Lots of polygons:
> releaf <-
> function(leaf)cbind(leaf[,1]+rnorm(1,0,50),leaf[,2]+rnorm(1,0,50))
> leaves <- replicate(500, releaf(leaf), simplify=FALSE)
>
> # Make into gpc.poly class:
> leaves <- lapply(leaves, as, "gpc.poly")
>
> # Make union .....
> system.time({
> leavesoutline <- union(leaves[[1]], leaves[[2]])
> for(i in 3:length(leaves))leavesoutline <- union(leavesoutline,
> leaves[[i]])
> })
> # about 1sec here.
>
> # Check it:
> plot(leavesoutline)
>
>
>
> thanks!
>
> Remko
>
>
> -------------------------------------------------
> Remko Duursma
> Research Lecturer
>
> Centre for Plants and the Environment
> University of Western Sydney
> Hawkesbury Campus
> Richmond NSW 2753
>
> Dept of Biological Science
> Macquarie University
> North Ryde NSW 2109
> Australia
>
> Mobile: +61 (0)422 096908
> www.remkoduursma.com
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]