[R] Tutorial on rgl Graphics

[Duncan Murdoch]

On 7/25/2008 11:01 AM, Tom La Bone wrote:
> After looking around a bit more I found the example I was looking for --
> plotlm3d, which I found on the R wiki
> 
>      http://wiki.r-project.org/rwiki/doku.php?id=graph_gallery:new-graphics
> 
> The original author was John Fox, and it was modified by Jose Claudio Faria
> and Duncan Murdoch. Below is a simplified version of the function and two
> examples. One example was presented as a test of the function and it works
> fine. The second example is the plot I want to make and I can't seem to get
> the scale on the x and y axes correct. Being unfamiliar with rgl, can
> someone provide a hint on how to get the scales right? Thanks for the help.
> 
> Tom
> 
> 
> 
> plotlm3d <- function (x, y, z,
>                       surface        = T,
>                       model          = 'z ~ x + y',
>                       simple.axes    = T,
>                       box            = F,
>                       xlab           = deparse(substitute(x)),
>                       ylab           = deparse(substitute(y)),
>                       zlab           = deparse(substitute(z)),
>                       surface.col    = c('blue', 'orange', 'red', 'green',
>                                          'magenta', 'cyan', 'yellow',
> 'gray', 'brown'),
>                       point.col      = 'yellow',
>                       grid.col       = material3d("color"),
>                       grid           = T,
>                       grid.lines     = 26,
>                       sphere.factor  = 1,
>                       threshold      = 0.01)
> {
>   require(rgl)
>   require(mgcv)
>   xlab; ylab; zlab
>   size <- max(c(x,y,z))/100 * sphere.factor
>   if (size > threshold)
>     spheres3d(x, y, z, color = point.col, radius = size)
>   else
>     points3d(x, y, z, color = point.col)
>   aspect3d(c(1, 1, 1))
>   if (surface) {
>     xvals <- seq(min(x), max(x), length = grid.lines)
>     yvals <- seq(min(y), max(y), length = grid.lines)
>     dat  <- expand.grid(x = xvals, y = yvals)
>     for (i in 1:length(model)) {
>       mod <- lm(formula(model[i]))
>       zhat <- matrix(predict(mod, newdata = dat), grid.lines, grid.lines)
>       surface3d(xvals, yvals, zhat, color = surface.col[i], alpha = 0.5, lit
> = F)
>       if (grid)
>         surface3d(xvals, yvals, zhat, color = grid.col, alpha = 0.5,
>         lit = F, front = 'lines', back = 'lines') }}
>   if(simple.axes) {
>     axes3d(c('x', 'y', 'z'))
>     title3d(xlab = xlab, ylab = ylab, zlab = zlab)
>   }
>   else
>     decorate3d(xlab = xlab, ylab = ylab, zlab = zlab, box = box)
> }
> 
> #This is an example of a 3D scatterplot that works fine
> x <- c( 274,  180,  375,  205,   86,  265,   98,  330,  195,   53,
>        430,  372,  236,  157,  370)
> y <- c(2450, 3254, 3802, 2838, 2347, 3782, 3008, 2450, 2137, 2560,
>       4020, 4427, 2660, 2088, 2605)
> z <- c( 162,  120,  223,  131,   67,  169,   81,  192,  116,   55,
>        252,  232,  144,  103,  212)
> open3d()
> plotlm3d(x, y, z,
>          surface = T,
>          model   = 'z ~ x + y',
>          xlab    = 'x',
>          ylab    = 'y',
>          zlab    = 'z')
> 
> #This is the plot I am trying to make - scales on x and y axes are wrong
> x <- c(0.3405,0.1220,0.1028,0.08451,0.05668,0.0345,0.003788,0.002121)
> y <- c(0.3460,0.1227,0.1097,0.09666,0.07677,0.06278,0.02168,0.01303)
> z <- c(2720,1150,1010,790,482,358,78,35)
> open3d()
> plotlm3d(x, y, z,
>          surface = T,
>          model   = 'z ~ x + y - 1',
>          xlab    = 'x',
>          ylab    = 'y',
>          zlab    = 'z')


I think you're seeing the effect of the inaccurate bounding box 
calculation described in ?spheres3d.  The problem is that you're asking 
for a sphere, but your axes are on wildly different scales.  The 
bounding box calculation for that situation fails.

You can work around this by telling rgl to ignore the extent of the 
spheres (via par3d(ignoreExtent=TRUE)), and plot some points in the 
corners of the bounding box you really want.  Set their alpha to 0 and 
they'll be invisible.

Some day I'll probably fix this, but it's likely to be a while.

Duncan Murdoch