如何在R中计算轮廓内的面积？

Question

我想知道在R中是否可以计算轮廓内的面积。

例如，从以下轮廓得到的面积：

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)

抱歉，我知道这段代码可能很复杂。如果它太难以阅读，请提供任何示例，您可以向我展示如何计算简单生成轮廓的面积。

感谢任何帮助。

- Burton Guster

---

你是指轮廓内部的区域吗？ - Phonon

轮廓线是如何定义的（以一个例子为例）...你想要什么区域？如果有边界的坐标，那就相对容易了。 - John

是的。轮廓内的区域。@John：我将把一个例子编辑到我的原始问题中。 - Burton Guster

好的，既然你提供了一个例子，我想知道你是想让该等高线下的体积为0，还是想要由该等高线定义的表面积。如果是表面积，你是指三维总面积还是仅指二维面积。 - John

@John：只需要计算二维区域。所以，如果我的轮廓看起来像一个圆形，我想要计算那个圆形对象的面积。 - Burton Guster

Answer 1

我假设您正在使用由contourLines返回的对象（每个级别都有x和y组件的未命名列表）。我原本希望在易于访问的位置找到它，但实际上是在一个PDF文件中找到了一个算法，我模糊地记得在http://finzi.psych.upenn.edu/R/library/PBSmapping/doc/PBSmapping-UG.pdf上看到过（请参见标为“-11-”的pdf第19页）（补充说明：维基百科关于“多边形”的文章引用了这篇关于测量公式的讨论：http://www.maa.org/pubs/Calc_articles/ma063.pdf，这证明了我使用abs()的理由。）

建立一个示例：

 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

在等级==130的区域（选择此区域是因为没有两个等级为130的区域，也不符合任何情节边界）如下：

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

Answer 2

感谢@DWin提供可重现的例子，以及（我最喜欢的R包！）和的作者们...

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

得到与@DWin相同的答案，这很令人安心。(可能是相同的算法，但在包中用Fortran例程实现...)

- Ben Bolker

---

我并不值得太多的赞誉，因为我只是去了help(contourlInes)页面并从那里获取了示例。@BurtonGuster应该随意给Ben打勾。我应该向Ben捐赠几千个积分，以表彰他在这里和R-help邮件列表上发布的所有优秀工作。同样也要感谢Hadley和Gabor Grothendeick。 - IRTFM