6

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.

joran
  • 169,992
  • 32
  • 429
  • 468
Burton Guster
  • 2,213
  • 8
  • 31
  • 29
  • 1
    Do you mean area *within* a contour? – Phonon Dec 06 '11 at 20:41
  • How is the contour defined (an example)... what area do you want? If you've got coordinates of a boundary it's relatively easy. – John Dec 06 '11 at 20:48
  • Yes. The area within a contour. @John: I'll edit an example into my original question. – Burton Guster Dec 06 '11 at 20:58
  • OK, so now that you've provided an example I'm wondering if you want the volume under that contour to 0 or if you want the area of the surface defined by that contour. And, if the area of the surface, do you mean the total area 3-dimensionally or just the 2d area. – John Dec 06 '11 at 21:20
  • @John: Just the 2-dimensional area. So if my contour looks like a circle, I'm looking for the area of that circle like object. – Burton Guster Dec 06 '11 at 21:31

2 Answers2

6

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
IRTFM
  • 258,963
  • 21
  • 364
  • 487
5

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 ...)

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • 2
    I don't deserve much credit since I just went to the `help(contourlInes)` page and grabbed the example there. @BurtonGuster should feel free to give Ben the check mark. I should probably donate a few thousand points to Ben for all the good work he has posted here and on the R-help mailing list. Likewise to Hadley and Gabor Grothendeick – IRTFM Dec 07 '11 at 02:37