Calculating the area of a convex hull

A reader recently posted a comment on my plotting convex hull post asking how to calculate the area of a convex hull. While this never occurred to me before, I decided to make a quick post on how to do it. To make this snippet of script work, you’ll need to install the package “sp” using the install.pacakages() function.

First I’ll create some dummy data, in this case, it will be a simple box (so we can make sure that the area is correct):

box.coords <- matrix(c(1, 1,
                       4, 1,
                       1, 4,
                       4, 4), nrow = 4, ncol = 2, byrow = T)
plot(box.coords[,1], box.coords[,2], pch = 20)

We have a square with sides of length 3. So, we know beforehand that we should get an area of 9. Next, we will use the chull() to extract the coordinates that make up the convex hull (in this case, all of the points, since we only have 4 points!). Also, we’ll need to close the convex hull by repeating the first row in the last row. The reason for this will be apparent in minute. Finally, we create a new matrix with the convex hull points, that includes the repeated coordinates to close off the convex hull:

box.hpts <- chull(x = box.coords[,1], y = box.coords[,2])
box.hpts <- c(box.hpts, box.hpts[1])
box.chull.coords <- box.coords[box.hpts,]

Finally, we load the “sp” package so we can use the function Polygon(). If you read the help files for Polygon(), you’ll see the first argument is the a two row matrix of coordinates (x, y), with the coordinates in the first repeated again in the last row (aha!) This will close off the polygon and allow Polygon() to calculate an area. To extract the area of the polygon, you need to access the attributes of a SpatialPolygon object, using the @area symbol:

chull.poly <- Polygon(box.chull.coords, hole=F)
chull.area <- chull.poly@area
> chull.area
[1] 9

Not too hard! You can do this with any polygon, as long as it is a valid SpatialPolygon object. Here is a quick examples using more realistic data:

library(sp)

x1 <- rnorm(100, 0.8, 0.3)
y1 <- rnorm(100, 0.8, 0.3)
hpts <- chull(x = x1, y = y1)
hpts <- c(hpts, hpts[1])
xy.coords <- cbind(x1, y1)
chull.coords <- xy.coords[hpts,]
chull.poly <- Polygon(chull.coords, hole=F)
chull.area <- chull.poly@area

Just remember, that calculating an area of a convex hull only makes sense if the units of the X and Y axes are the same. Hope this answers the question!

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s