[R-sig-Geo] calculate the area of an alpha shape

[Ashton Shortridge]

Dear André,

You're right! I assumed something that wasn't true - that $x in the ashape.obj 
were the points on the hull. That's what you get for free help, I guess.... 

Fortunately this seems to work. Notice how I index $x by the alpha.extremes 
field to extract just the points I want for bds:

library(sp)
bds <- achull.obj$ashape.obj$x[achull.obj$ashape.obj$alpha.extremes,] 
bds <- rbind(bds, bds[1,])          # close the ring
ashape <- Polygon(bds)            # convert to sp Polygon  
print(ashape at area)

# The following plots the points in ashape superimposed on the hulls
plot(achull.obj, col = "blue")
plot(achull.obj$ashape.obj, add = TRUE, col = "red")
points(coordinates(ashape), pch=16, cex=2, col='gold')


On Monday, November 04, 2013 04:32:42 PM Proosdij, Andre van wrote:
> Many thanks Ashton!
> 
> It does help a lot, but I'm afraid, it's not completely correct yet. One
> would expect that the area of the alpha shape is larger than the area of
> the alpha-convex hull (straight line segments instead of arcs). However,
> applying your suggestion, the area is smaller! If I plot the points and the
> alpha shape, the plot is correct. The area of the alpha shape seems to be
> calculated a polygon that includes the point (1.5, 2.0). See the updated
> example below. I don't know what happens, but apparently, the list of
> points that defines the polygon is incorrect. Any thoughts of how to adjust
> this?
> 
> Best,
> 
> André
> 
> 
> x <- c(1,1,1.5,2.8,3,3) # original set of points
> y <- c(1,3,2,3,4,1) # original set of points
> #x <- c(1,1,2.8,3,3) # points without (1.5, 2.0)
> #y <- c(1,3,3,4,1) # points without (1.5, 2.0)
> z <- cbind(x,y)
> plot(z)
> 
> ### Calculate the alpha convex hull
> library(alphahull)
> alpha <- 3 # Define alpha > 0
> achull.obj <- ahull(z, alpha = alpha)
> 
> ### Plot the sampled points and the alpha-convex hull object
> plot(achull.obj, add = TRUE, col = "blue")
> 
> ### Get the length of the alpha-convex hull
> achulllength <- achull.obj$length
> achulllength
> 
> ### Get the area of the alpha-convex hull
> achullarea <- areaahull(achull.obj)
> achullarea
> 
> ### Plot the alpha shape
> plot(achull.obj$ashape.obj, add = TRUE, col = "red")
> 
> ### Get the length of the alpha shape
> ashapelength <- achull.obj$ashape.obj$length
> ashapelength
> 
> ### Get the area of the alpha shape
> library(sp)
> bds <- achull.obj$ashape.obj$x  # matrix of coordinates in ashape
> bds <- rbind(bds, bds[1,]) # close the ring
> ashape <- Polygon(bds) # convert to sp Polygon
> ashapearea <- ashape at area
> ashapearea
> 
> 
> 
> 
> Ir. A.S.J. van Proosdij
> PhD Student
> t   +31 317 48 1198 
> e  andre.vanproosdij at wur.nl
> 
> Naturalis Biodiversity Center, section NHN
> Biosystematics Group, Wageningen University
> Gen. Foulkesweg 37, 6703 BL Wageningen, the Netherlands
> 
> www.bis.wur.nl
> 
> -----Original Message-----
> From: Ashton Shortridge [mailto:ashton at msu.edu]
> Sent: maandag 4 november 2013 14:35
> To: r-sig-geo at r-project.org
> Cc: Proosdij, Andre van
> Subject: Re: [R-sig-Geo] calculate the area of an alpha shape
> 
> Dear André,
> 
> There may be a more elegant solution, but this works fine:
> library(sp)
> bds <- achull.obj$ashape.obj$x  # matrix of coordinates in ashape
> bds <- rbind(bds, bds[1,])          # close the ring
> ashape <- Polygon(bds)            # convert to sp Polygon
> print(ashape at area)
> 
> On Monday, November 04, 2013 11:58:26 AM Proosdij, Andre van wrote:
> > Dear colleagues,
> > 
> > I wish to calculate the area covered by the alpha shape based on a
> > sample of points in a 2D plane. For the alpha-convex hull this works
> > fine using the function ahull in the alphahull package. However, I
> > want to get the area of the slightly different alpha shape. How can I
> > do this? Below a short example.
> > 
> > Many thanks in advance for your answers and ideas!
> > 
> > 
> > André

-----
Ashton Shortridge
Associate Professor			ashton at msu.edu
Dept of Geography			http://www.msu.edu/~ashton
235 Geography Building		ph (517) 432-3561
Michigan State University		fx (517) 432-1671