我需要在 R 中找到从表面到特定轮廓的体积。以 R 帮助文件中的轮廓示例:
x <- 10*1:nrow(volcano)
y <- 10*1:ncol(volcano)
contour(x,y,volcano)
给定结果图,我如何找到从特定轮廓线到表面的体积。
在实践中,我将使用 bkde2D 来获取散点图的密度图。由此我可以制作等高线图,但我想确定结果图中由各种密度截止值定义的体积。
函数contour
只绘制轮廓线但不返回任何值。您需要使用的是功能contourLines
。
cL <- contourLines(x,y,volcano)
从那里,您可以通过以下方式计算每条等高线的面积:
area <- rep(0,length(cL))
for(i in 1:length(cL)){
d <- data.frame(cL[[i]]$x,cL[[i]]$y)
sa <- sb <- 0
for(j in 1:(nrow(d)-1)){
sa <- sa+d[j,1]*d[j+1,2]
sb <- sb+d[j,2]*d[j+1,1]
}
area[i] <- abs((sa-sb)/2)
}
area
[1] 1.413924e+05 3.109685e+04 2.431528e+04 2.049473e+04 6.705976e+04 3.202145e+05 1.720469e+03
[8] 2.926802e+05 2.335421e+05 1.834791e+05 1.326162e+05 4.672784e+02 9.419792e+04 5.121851e+03
[15] 5.126860e+04 3.660862e-01 1.216750e+03 2.051307e+04 4.670745e+02 4.146927e+03
现在,如果你想要两条等高线之间的音量(比如 120 和 130 级之间):
level1 <- 120
level2 <- 130
levels <- unlist(lapply(cL,function(x)x$level))
base <- (1:length(cL))[level==level1]
top <- (1:length(cL))[level==level2]
vol <- (level[top]-level[base])*(area[base]+area[top])/2
vol
[1] 2631111
这就是我所能做到的,因为如果下一条等高线被分成几个扇区,我不知道如何进行。