0

I've been using the 'ks' package along with the 'rgl' package to produce 3D kernel density estimates and 3D plots of these. This first part has worked out fine (brief example below). What I can't figure out is if it's possible to extract the values of the kernels for the given xyz locations used to build the kernels in the first place. In other words, extract the values for points in a 3D plot, akin to the extract command used for 2D surfaces in the 'raster' package. Does anyone have experience doing something like this that can point me in the right direction? Thanks much. -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"
dxj
  • 513
  • 4
  • 6
  • Convert the fhat2 to raster, probably there is an "image() xyz" list in there. What is str(fhat2)? – mdsumner Aug 04 '14 at 22:19
  • Added str(fhat2) as requested by @mdsummer. But wouldn't converting to a raster or using image() cause me to lose the 3D information? In my case the z is a depth value and the kernel density values are produced in Lat/Lon/Depth space. So maybe this is better described as 4D? It seems like image() would flatten this information which is not what I need. – dxj Aug 05 '14 at 00:23

2 Answers2

1

Try this

## 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])

Now we are happy because raster power is good. 

plot(b)

## note this is a matrix with nrows = nrow(fhat$x), ncols = nlayers(b)
extract(b, fhat$x[,1:2])
mdsumner
  • 29,099
  • 6
  • 83
  • 91
  • Thanks @mdsumner. So if I'm understanding this correctly, what you're doing is building a series of raster surfaces perpendicular to the z axis - in this case,petal.length. Each raster surface consists of gridded values taken from the array fhat$estimate (looks like kde builds 51 of these). You then extract the value of each raster surface for a given xy (here, sepal length and sepal width) coordinate, regardless of the z value. So then I'm left with a kernel density value for a given xy across each of the 51 slices. Does that sound right? – dxj Aug 06 '14 at 00:52
  • Yes you have to do the sub-z yourself, it's not built in but easy enough if you need more help. – mdsumner Aug 06 '14 at 05:54
0

The answer may also lie in eval.points. Researching more it looks like you can enter your own points here, so you can potentially enter the points used to build the kde or an entirely new set of points.

dxj
  • 513
  • 4
  • 6