1

I am using the ks package from R to estimate 2d space utilization using distance and depth information. What I would like to do is to use the 95% contour output to get the maximum vertical and horizontal distance. So essentially, I want to be able to get the dimensions or measurements of the resulting 95% contour.

Here is a piece of code with as an example,

require(ks)

dist<-c(1650,1300,3713,3718)
depth<-c(22,19.5,20.5,8.60)

dd<-data.frame(cbind(dist,depth))

## auto bandwidth selection
H.pi2<-Hpi(dd,binned=TRUE)*1
ddhat<-kde(dd,H=H.pi2)

plot(ddhat,cont=c(95),lwd=1.5,display="filled.contour2",col=c(NA,"palegreen"),
     xlab="",ylab="",las=1,ann=F,bty="l",xaxs="i",yaxs="i",
     xlim=c(0,max(dd[,1]+dd[,1]*0.4)),ylim=c(60,-3))

Any information about how to do this will be very helpful. Thanks in advance,

user1626688
  • 1,583
  • 4
  • 18
  • 27

2 Answers2

1

To create a 95% contour polygon from your 'kde' object:

library(raster)
im.kde <- image2Grid (list(x = ddhat$eval.points[[1]], y =  ddhat$eval.points[[2]], z = ddhat$estimate))
kr <- raster(im.kde)

It is likely that one will want to resample this raster to a higher resolution before constructing polygons, and include the following two lines, before creation of the polygon object:

new.rast <- raster(extent(im.kde),res = c(50,50))
kr <- resample(kr, new.rast)

bin.kr <- kr
bin.kr[bin.kr < contourLevels(k, prob = 0.05)]<-NA
bin.kr[bin.kr > 0]<-1
k.poly<-rasterToPolygons(bin.kr,dissolve=T)

Note that the results are similar, but not identical, to Hawthorne Beier's GME function 'kde'. He does use the kde function from ks, but must do something slightly different for the output polygon.

Dave Gregovich
  • 11
  • 3
0

At the moment I'm going for the "any information" prize rather than attempting a final answer. The ks:::plot.kde function dispatches to ks:::plotkde.2d in this case. It works its magic through side effects and I cannot get these functions to return values that can be inspected in code. You would need to hack the plotkde.2d function to return the values used to plot the contour lines. You can visualize what is in ddhat$estimate with:

persp(ddhat$estimate)

It appears that contourLevels examines the estimate-matrix and finds the value at which greater than the specified % of the total density will reside.

> contourLevels(ddhat, 0.95)
         95% 
1.891981e-05

And then draws the contout based on which values exceed that level. (I just haven't found the code that does that yet.)

IRTFM
  • 258,963
  • 21
  • 364
  • 487