CODEWITHSUNDEEP
rareacontour
Burton Guster
Burton Guster

Reputation: 2223

How can I calculate the area within a contour in R?

I'm wondering if it is possible to caclulate the area within a contour in R.

For example, the area of the contour that results from:

sw<-loess(m~l+d)
mypredict<-predict(sw, fitdata) # Where fitdata is a data.frame of an x and y matrix

contour(x=seq(from=-2, to=2, length=30), y=seq(from=0, to=5, length=30), z=mypredict)

Sorry, I know this code might be convoluted. If it's too tough to read. Any example where you can show me how to calculate the area of a simply generated contour would be helpful.

Thanks for any help.

Upvotes: 6

Views: 3218

Answers (2)

Ben Bolker
Ben Bolker

Reputation: 226771

Thanks to @DWin for reproducible example, and to the authors of sos (my favourite R package!) and splancs ...

library(sos)
findFn("area polygon compute")
library(splancs)
with(clines[[9]],areapl(cbind(x,y)))

Gets the same answer as @DWin, which is comforting. (Presumably it's the same algorithm, but implemented within a Fortran routine in the splancs package ...)

Upvotes: 5

IRTFM
IRTFM

Reputation: 263461

I'm going to assume you are working with an object returned by contourLines. (An unnamed list with x and y components at each level.) I was expecting to find this in an easy to access location but instead found a pdf file that provided an algorithm which I vaguely remember seeing http://finzi.psych.upenn.edu/R/library/PBSmapping/doc/PBSmapping-UG.pdf (See pdf page 19, labeled "-11-") (Added note: The Wikipedia article on "polygon" cites this discussion of the Surveyors' Formula: http://www.maa.org/pubs/Calc_articles/ma063.pdf , which justifies my use of abs().)

Building an example:

 x <- 10*1:nrow(volcano)
 y <- 10*1:ncol(volcano)
contour(x, y, volcano); 
clines <- contourLines(x, y, volcano)
x <- clines[[9]][["x"]]
 y <- clines[[9]][["y"]]
 level <- clines[[9]][["level"]]
 level
#[1] 130

The area at level == 130 (chosen because there are not two 130 levels and it doesn't meet any of the plot boundaries) is then:

A = 0.5* abs( sum( x[1:(length(x)-1)]*y[2:length(x)] - y[1:(length(x)-1)]*x[2:length(x)] ) )
A
#[1] 233542.1

Upvotes: 6

Related Questions