[r] R: Plotting a 3D surface from x, y, z

imagine I have a 3 columns matrix
x, y, z where z is a function of x and y.

I know how to plot a "scatter plot" of these points with plot3d(x,y,z)

But if I want a surface instead I must use other commands such as surface3d The problem is that it doesn't accept the same inputs as plot3d it seems to need a matrix with

(nº elements of z) = (n of elements of x) * (n of elements of x)

How can I get this matrix? I've tried with the command interp, as I do when I need to use contour plots.

How can I plot a surface directly from x,y,z without calculating this matrix? If I had too many points this matrix would be too big.

cheers

If your x and y coords are not on a grid then you need to interpolate your x,y,z surface onto one. You can do this with kriging using any of the geostatistics packages (geoR, gstat, others) or simpler techniques such as inverse distance weighting.

I'm guessing the 'interp' function you mention is from the akima package. Note that the output matrix is independent of the size of your input points. You could have 10000 points in your input and interpolate that onto a 10x10 grid if you wanted. By default akima::interp does it onto a 40x40 grid:

require(akima)
require(rgl)

x = runif(1000)
y = runif(1000)
z = rnorm(1000)
s = interp(x,y,z)
> dim(s$z)
[1] 40 40
surface3d(s$x,s$y,s$z)

That'll look spiky and rubbish because its random data. Hopefully your data isnt!

rgl is great, but takes a bit of experimentation to get the axes right.

If you have a lot of points, why not take a random sample from them, and then plot the resulting surface. You can add several surfaces all based on samples from the same data to see if the process of sampling is horribly affecting your data.

So, here is a pretty horrible function but it does what I think you want it to do (but without the sampling). Given a matrix (x, y, z) where z is the heights it will plot both the points and also a surface. Limitations are that there can only be one z for each (x,y) pair. So planes which loop back over themselves will cause problems.

The plot_points = T will plot the individual points from which the surface is made - this is useful to check that the surface and the points actually meet up. The plot_contour = T will plot a 2d contour plot below the 3d visualization. Set colour to rainbow to give pretty colours, anything else will set it to grey, but then you can alter the function to give a custom palette. This does the trick for me anyway, but I'm sure that it can be tidied up and optimized. The verbose = T prints out a lot of output which I use to debug the function as and when it breaks.

plot_rgl_model_a <- function(fdata, plot_contour = T, plot_points = T, 
                             verbose = F, colour = "rainbow", smoother = F){
  ## takes a model in long form, in the format
  ## 1st column x
  ## 2nd is y,
  ## 3rd is z (height)
  ## and draws an rgl model

  ## includes a contour plot below and plots the points in blue
  ## if these are set to TRUE

  # note that x has to be ascending, followed by y
  if (verbose) print(head(fdata))

  fdata <- fdata[order(fdata[, 1], fdata[, 2]), ]
  if (verbose) print(head(fdata))
  ##
  require(reshape2)
  require(rgl)
  orig_names <- colnames(fdata)
  colnames(fdata) <- c("x", "y", "z")
  fdata <- as.data.frame(fdata)

  ## work out the min and max of x,y,z
  xlimits <- c(min(fdata$x, na.rm = T), max(fdata$x, na.rm = T))
  ylimits <- c(min(fdata$y, na.rm = T), max(fdata$y, na.rm = T))
  zlimits <- c(min(fdata$z, na.rm = T), max(fdata$z, na.rm = T))
  l <- list (x = xlimits, y = ylimits, z = zlimits)
  xyz <- do.call(expand.grid, l)
  if (verbose) print(xyz)
  x_boundaries <- xyz$x
  if (verbose) print(class(xyz$x))
  y_boundaries <- xyz$y
  if (verbose) print(class(xyz$y))
  z_boundaries <- xyz$z
  if (verbose) print(class(xyz$z))
  if (verbose) print(paste(x_boundaries, y_boundaries, z_boundaries, sep = ";"))

  # now turn fdata into a wide format for use with the rgl.surface
  fdata[, 2] <- as.character(fdata[, 2])
  fdata[, 3] <- as.character(fdata[, 3])
  #if (verbose) print(class(fdata[, 2]))
  wide_form <- dcast(fdata, y ~ x, value_var = "z")
  if (verbose) print(head(wide_form))
  wide_form_values <- as.matrix(wide_form[, 2:ncol(wide_form)])  
  if (verbose) print(wide_form_values)
  x_values <- as.numeric(colnames(wide_form[2:ncol(wide_form)]))
  y_values <- as.numeric(wide_form[, 1])
  if (verbose) print(x_values)
  if (verbose) print(y_values)
  wide_form_values <- wide_form_values[order(y_values), order(x_values)]
  wide_form_values <- as.numeric(wide_form_values)
  x_values <- x_values[order(x_values)]
  y_values <- y_values[order(y_values)]
  if (verbose) print(x_values)
  if (verbose) print(y_values)

  if (verbose) print(dim(wide_form_values))
  if (verbose) print(length(x_values))
  if (verbose) print(length(y_values))

  zlim <- range(wide_form_values)
  if (verbose) print(zlim)
  zlen <- zlim[2] - zlim[1] + 1
  if (verbose) print(zlen)

  if (colour == "rainbow"){
    colourut <- rainbow(zlen, alpha = 0)
    if (verbose) print(colourut)
    col <- colourut[ wide_form_values - zlim[1] + 1]
    # if (verbose) print(col)
  } else {
    col <- "grey"
    if (verbose) print(table(col2))
  }


  open3d()
  plot3d(x_boundaries, y_boundaries, z_boundaries, 
         box = T, col = "black",  xlab = orig_names[1], 
         ylab = orig_names[2], zlab = orig_names[3])

  rgl.surface(z = x_values,  ## these are all different because
              x = y_values,  ## of the confusing way that 
              y = wide_form_values,  ## rgl.surface works! - y is the height!
              coords = c(2,3,1),
              color = col,
              alpha = 1.0,
              lit = F,
              smooth = smoother)

  if (plot_points){
    # plot points in red just to be on the safe side!
    points3d(fdata, col = "blue")
  }

  if (plot_contour){
    # plot the plane underneath
    flat_matrix <- wide_form_values
    if (verbose) print(flat_matrix)
    y_intercept <- (zlim[2] - zlim[1]) * (-2/3) # put the flat matrix 1/2 the distance below the lower height 
    flat_matrix[which(flat_matrix != y_intercept)] <- y_intercept
    if (verbose) print(flat_matrix)

    rgl.surface(z = x_values,  ## these are all different because
                x = y_values,  ## of the confusing way that 
                y = flat_matrix,  ## rgl.surface works! - y is the height!
                coords = c(2,3,1),
                color = col,
                alpha = 1.0,
                smooth = smoother)
  }
}

The add_rgl_model does the same job without the options, but overlays a surface onto the existing 3dplot.

add_rgl_model <- function(fdata){

  ## takes a model in long form, in the format
  ## 1st column x
  ## 2nd is y,
  ## 3rd is z (height)
  ## and draws an rgl model

  ##
  # note that x has to be ascending, followed by y
  print(head(fdata))

  fdata <- fdata[order(fdata[, 1], fdata[, 2]), ]

  print(head(fdata))
  ##
  require(reshape2)
  require(rgl)
  orig_names <- colnames(fdata)

  #print(head(fdata))
  colnames(fdata) <- c("x", "y", "z")
  fdata <- as.data.frame(fdata)

  ## work out the min and max of x,y,z
  xlimits <- c(min(fdata$x, na.rm = T), max(fdata$x, na.rm = T))
  ylimits <- c(min(fdata$y, na.rm = T), max(fdata$y, na.rm = T))
  zlimits <- c(min(fdata$z, na.rm = T), max(fdata$z, na.rm = T))
  l <- list (x = xlimits, y = ylimits, z = zlimits)
  xyz <- do.call(expand.grid, l)
  #print(xyz)
  x_boundaries <- xyz$x
  #print(class(xyz$x))
  y_boundaries <- xyz$y
  #print(class(xyz$y))
  z_boundaries <- xyz$z
  #print(class(xyz$z))

  # now turn fdata into a wide format for use with the rgl.surface
  fdata[, 2] <- as.character(fdata[, 2])
  fdata[, 3] <- as.character(fdata[, 3])
  #print(class(fdata[, 2]))
  wide_form <- dcast(fdata, y ~ x, value_var = "z")
  print(head(wide_form))
  wide_form_values <- as.matrix(wide_form[, 2:ncol(wide_form)])  
  x_values <- as.numeric(colnames(wide_form[2:ncol(wide_form)]))
  y_values <- as.numeric(wide_form[, 1])
  print(x_values)
  print(y_values)
  wide_form_values <- wide_form_values[order(y_values), order(x_values)]
  x_values <- x_values[order(x_values)]
  y_values <- y_values[order(y_values)]
  print(x_values)
  print(y_values)

  print(dim(wide_form_values))
  print(length(x_values))
  print(length(y_values))

  rgl.surface(z = x_values,  ## these are all different because
              x = y_values,  ## of the confusing way that 
              y = wide_form_values,  ## rgl.surface works!
              coords = c(2,3,1),
              alpha = .8)
  # plot points in red just to be on the safe side!
  points3d(fdata, col = "red")
}

So my approach would be to, try to do it with all your data (I easily plot surfaces generated from ~15k points). If that doesn't work, take several smaller samples and plot them all at once using these functions.

You could look at using Lattice. In this example I have defined a grid over which I want to plot z~x,y. It looks something like this. Note that most of the code is just building a 3D shape that I plot using the wireframe function.

The variables "b" and "s" could be x or y.

require(lattice)

# begin generating my 3D shape
b <- seq(from=0, to=20,by=0.5)
s <- seq(from=0, to=20,by=0.5)
payoff <- expand.grid(b=b,s=s)
payoff$payoff <- payoff$b - payoff$s
payoff$payoff[payoff$payoff < -1] <- -1
# end generating my 3D shape


wireframe(payoff ~ s * b, payoff, shade = TRUE, aspect = c(1, 1),
    light.source = c(10,10,10), main = "Study 1",
    scales = list(z.ticks=5,arrows=FALSE, col="black", font=10, tck=0.5),
    screen = list(z = 40, x = -75, y = 0))

You can use the function outer() to generate it.

Have a look at the demo for the function persp(), which is a base graphics function to draw perspective plots for surfaces.

Here is their first example:

x <- seq(-10, 10, length.out = 50)  
y <- x  
rotsinc <- function(x,y) {
    sinc <- function(x) { y <- sin(x)/x ; y[is.na(y)] <- 1; y }  
    10 * sinc( sqrt(x^2+y^2) )  
}

z <- outer(x, y, rotsinc)  
persp(x, y, z)

The same applies to surface3d():

require(rgl)  
surface3d(x, y, z)

Maybe is late now but following Spacedman, did you try duplicate="strip" or any other option?

x=runif(1000)
y=runif(1000)
z=rnorm(1000)
s=interp(x,y,z,duplicate="strip")
surface3d(s$x,s$y,s$z,color="blue")
points3d(s)