如何使用“ks”和“rgl”从 R 中构建的 3D 内核密度图中提取值答案

【问题标题】：How to extract values from a 3D kernel density plot built in R using 'ks' and 'rgl'如何使用“ks”和“rgl”从 R 中构建的 3D 内核密度图中提取值
【发布时间】：2014-08-04 21:18:26
【问题描述】：

我一直在使用“ks”包和“rgl”包来生成 3D 内核密度估计和这些的 3D 图。第一部分效果很好（下面的简短示例）。我无法弄清楚的是，是否可以首先提取用于构建内核的给定 xyz 位置的内核值。换句话说，提取 3D 图中点的值，类似于用于 'raster' 包中的 2D 表面的提取命令。有没有人有做这样的事情的经验可以为我指明正确的方向？非常感谢。 -DJ

library("rgl")
library("ks")

# call the plug-in bandwidth estimator
H.pi <- Hpi(b,binned=TRUE) ## b is a matrix of x,y,z points

# calculate the kernel densities
fhat2 <- kde(b, H=H.pi)

#plot the 50% and 95% kernels in gray and blue
plot(fhat2,cont=c(50,95),colors=c("gray","blue"),drawpoints=TRUE
    ,xlab="", ylab="", zlab="",size=2, ptcol="white", add=FALSE, box=TRUE, axes=TRUE) 




#Structure of fhat2. Original df consists of ~6000 points.  

List of 9
 $ x          : num [1:6173, 1:3] -497654 -497654 -497929 -498205 -498205 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:6173] "50" "57" "70" "73" ...
  .. ..$ : chr [1:3] "x" "max_dep" "y"

$ eval.points:List of 3
  ..$ : num [1:51] -550880 -546806 -542733 -538659 -534586 ...
  ..$ : num [1:51] -7.9 -4.91 -1.93 1.06 4.05 ...
  ..$ : num [1:51] -376920 -374221 -371522 -368823 -366124 ...

$ estimate   : num [1:51, 1:51, 1:51] 0 0 0 0 0 ...

$ H          : num [1:3, 1:3] 3.93e+07 -2.97e+03 8.95e+06 -2.97e+03 2.63e+01 ...

$ gridtype   : chr [1:3] "linear" "linear" "linear"

$ gridded    : logi TRUE

$ binned     : logi FALSE

$ names      : chr [1:3] "x" "max_dep" "y"

$ w          : num [1:6173] 1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "class")= chr "kde"

【问题讨论】：

将 fhat2 转换为光栅，可能那里有一个“image() xyz”列表。什么是 str(fhat2)？
应@mdsummer 的要求添加了 str(fhat2)。但是转换为光栅或使用 image() 不会导致我丢失 3D 信息吗？在我的例子中，z 是一个深度值，内核密度值是在 Lat/Lon/Depth 空间中产生的。所以也许这更好地描述为4D？似乎 image() 会使这些信息变平，这不是我需要的。

标签： r 3d extract rgl kernel-density

【解决方案1】：

试试这个

## from ?plot.kde
library(ks)
library(MASS)
 data(iris)

 ## trivariate example
 fhat <- kde(x=iris[,1:3])

## this indicates the orientation
image(fhat$eval.points[[1]], fhat$eval.points[[2]], apply(fhat$estimate, 1:2, sum))
points(fhat$x[,1:2])

library(raster)

## convert to RasterBrick from raw array
## with correct orientation relative to that of ?base::image
b <- brick(fhat$estimate[,ncol(fhat$estimate):1,], 
    xmn = min(fhat$eval.points[[1]]), xmx = max(fhat$eval.points[[1]]), ymn = min(fhat$eval.points[[2]]), ymx = max(fhat$eval.points[[2]]), 
    transpose = TRUE)

## check orientation
plot(calc(b, sum))
points(fhat$x[,1:2])

现在我们很高兴，因为光栅功率很好。

plot(b)

## note this is a matrix with nrows = nrow(fhat$x), ncols = nlayers(b)
extract(b, fhat$x[,1:2])

【讨论】：

谢谢@mdsumner。因此，如果我正确理解这一点，那么您正在做的是构建一系列垂直于 z 轴的光栅表面 - 在本例中为花瓣长度。每个栅格表面都由从数组 fhat$estimate 中获取的网格值组成（看起来 kde 构建了其中的 51 个）。然后提取给定 xy（此处为萼片长度和萼片宽度）坐标的每个光栅表面的值，而不管 z 值如何。因此，对于 51 个切片中的每一个切片，我都会得到一个给定 xy 的内核密度值。听起来对吗？
是的，你必须自己做 sub-z，它不是内置的，但如果你需要更多帮助，它很容易。

【解决方案2】：

答案也可能在 eval.points 中。进行更多研究，您似乎可以在此处输入自己的积分，因此您可以输入用于构建 kde 的积分或一组全新的积分。

【讨论】：