How to create 3D - MATLAB style - surface plots in R

Question

I find it challenging to create aesthetically pleasing 3D surfaces in R. I am familiar with the solutions (persp, image, wireframe, lattice, rgl and several other solutions in other questions in SO), but the results are not nice.

Is it possible to create 3D surface plots in R like in MATLAB?

Here is the MATLAB code
    % Create a grid of x and y points
    points = linspace(-2, 0, 20);
    [X, Y] = meshgrid(points, -points);

    % Define the function Z = f(X,Y)
    Z = 2./exp((X-.5).^2+Y.^2)-2./exp((X+.5).^2+Y.^2);

    % "phong" lighting is good for curved, interpolated surfaces. "gouraud"
    % is also good for curved surfaces
    surf(X, Y, Z); view(30, 30);
    shading interp;
    light;
    lighting phong;
    title('lighting phong', 'FontName', 'Courier', 'FontSize', 14);

The plot is modern, colorful, aesthetically pleasing, the code syntax is very readable.

Is this possible in base R?

Answer 1

jet.colors is the R-answer to one of hte Matlab color palettes:

points = seq(-2, 0, length=20)
#create a grid
XY = expand.grid(X=points,Y=-points)
# A z-function 
Zf <- function(X,Y){
     2./exp((X-.5)^2+Y^2)-2./exp((X+.5)^2+Y^2);
     }
# populate a surface
Z <- Zf(XY$X, XY$Y)
zlim <- range(Z)
zlen <- zlim[2] - zlim[1] + 1

jet.colors <-   # function from grDevices package
      colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                      "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
colorzjet <- jet.colors(100)  # 100 separate color 
require(rgl)
open3d()
rgl.surface(x=points, y=matrix(Z,20), 
            coords=c(1,3,2),z=-points, 
            color=colorzjet[ findInterval(Z, seq(min(Z), max(Z), length=100))] )
axes3d()
rgl.snapshot("copyMatlabstyle.png")

I will admit that getting the colors to line up with the "Z-axis" (which is actually the rgl y-axis) seemed very unintuitive. If you want the shiny, specular effect that Matlab delivers you can play with the angle of illumination.

You can also add or remove lighting:

clear3d(type = "lights")
light3d(theta=0, phi=0)
light3d(theta=0, phi=0)  # twice as much light.

After:

 grid3d("x")
 grid3d("y")
 grid3d("z")

 rgl.snapshot("copyMatlabstyle3.png")

You could have put the y-grid "behind" the surface with:

grid3d("y+")

Similar tweaks to the axes3d or axis3d calls could move the location of the scales.

For further examples, look at http://rgm3.lab.nig.ac.jp/RGM/R_image_list and search for 'plot3d' which brings up examples of the R2BayesX::plot3d function, Look at Karline Soetaert's plot3D package vignette, "50 ways to plot a volcano"