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

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

我想知道是否有可能在R中的轮廓内caclulate区域。

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.

谢谢你的帮助。

2 个解决方案

#1

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

感谢@DWin的可重复性示例,以及sos(我最喜欢的R包!)和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 ...)

得到与@DWin相同的答案,这很令人欣慰。 (据推测它是相同的算法,但在splancs包中的Fortran例程中实现...)

#2

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

我将假设您正在使用contourLines返回的对象。 (一个在每个级别都有x和y组件的未命名列表。)我希望在一个易于访问的位置找到它,但是找到了一个提供算法的pdf文件,我依旧记得看到http://finzi.psych.upenn .edu / R / library / PBSmapping / doc / PBSmapping-UG.pdf(参见pdf第19页,标记为“-11-”)(补充说明:*关于“多边形”的文章引用了对Surveyors公式的讨论:http ://www.maa.org/pubs/Calc_articles/ma063.pdf,证明我使用了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:

级别== 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

#1