[R] Draw level lines on the surface of a bivariate function

[Duncan Murdoch]

On 3/11/2006 12:17 PM, Cuvelier Etienne wrote:
> Hello,
> Is it possible to draw level lines  on the surface of a bivariate function?
> 
> In the following example, to draw surface and levels lines for a 
> multivariate normal law,
> I use persp, trans3d, contourLines and lines,
> but if  the lines are correctly drawn, some parts of them are, of 
> course, visible
> even if they are drawn on a non visible "face".
> 
> Any suggestion to avoid this problem ?

You could try to work out which parts of the contour should be hidden, 
based on their direction---but that's probably hard.

You could try the rgl package, but you'll need some fiddling, because 
when you draw lines at the same place you've drawn a surface, it's 
mostly rounding error that determines which one shows up.  rgl also 
doesn't do the automatic rescaling that persp does.  However, this code 
does something pretty close to what you want:

library(rgl)
open3d()
bg3d("white")
surface3d(x,y,20*z,col='red')
for (i in 1:length(lev))
   lines3d(lev[[i]]$x, lev[[i]]$y, 20*lev[[i]]$level+0.01, col='black')
bbox3d(back='lines',xat=pretty(range(x)), yat=pretty(range(y)), 
zat=pretty(range(z))*20, zlab=pretty(range(z)))

Duncan Murdoch

> Thank you
> 
> Etienne
> 
> Example :
> 
> trans3d <- function(x,y,z, pmat)
> {
>    tr <- cbind(x,y,z,1) %*% pmat
>    list(x = tr[,1]/tr[,4], y= tr[,2]/tr[,4])
> }
> 
> mean1 = c(2,4)
> cov1 = matrix(c(0.7, 0.2,0.2,0.7), ncol=2)
> 
> x = seq(0, 5, by = 0.1)
> y = seq(0,10, by= 0.1)
> z = matrix(nrow = length(x), ncol=length(y))
> 
> for(i in 1:length(x))
>   for(j in 1:length(y))
>   {
>     z[i,j]=dmvnorm(c(x[i],y[j]),mean1,cov1)
>   }
> 
> pmat1 = persp(x, y, z, col= "red",shade = 0.25, border=NA)
> 
> lev =contourLines(x,y,z, nlevel = 10)
> for(i in 1:length(lev))
>   lines(trans3d(lev[[i]]$x, lev[[i]]$y, lev[[i]]$level, pmat1), col="black")
> 
> 
> --
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html