R:用x, y, z绘制三维曲面。

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

假设我有一个3列矩阵x, y, z, z是x和y的函数。

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

我知道如何用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

但是，如果我想要一个表面，我必须使用其他的命令，比如surface3d，问题是它不接受与plot3d相同的输入，它似乎需要一个矩阵。

(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.

我怎么能在不计算这个矩阵的情况下直接画出x,y,z的曲面?如果我有太多的点这个矩阵会太大。

cheers

干杯

5 个解决方案

#1

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.

如果你的x和y坐标不在网格上，那么你需要把你的x,y,z曲面插到一个上面。您可以使用kriging使用任何地理统计包(geoR, gstat，其他)或更简单的技术，比如反向距离加权。

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:

我猜你提到的“interp”功能是来自akima软件包。注意，输出矩阵与输入点的大小无关。你可以在你的输入中有10000个点，如果你想的话，可以在10x10网格上进行插值。默认情况下，akima::interp在40x40网格上:

> 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!

因为它的随机数据，看起来会很乱，很垃圾。希望你的数据不是!

#2

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.

看一下函数persp()的演示，它是一个基本的图形函数，用于绘制曲面的透视图。

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():

同样适用于surface3d():

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

#3

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

rgl很棒，但是需要一些实验才能得到正确的坐标轴。

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.

这是一个很糟糕的函数但它做的是我认为你想要它做的(但是没有抽样)。给定一个矩阵(x, y, z) z是高度，它会同时画出点和曲面。限制是每个(x,y)对只能有一个z。因此，对自己进行循环的飞机会造成问题。

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_points = T将绘制出表面的各个点——这对于检查表面和点实际相遇是有用的。plot_contour = T将在3d可视化下绘制2d等高线图。将颜色设置为彩虹以提供漂亮的颜色，任何其他的颜色都会将其设置为灰色，但是您可以修改该函数以提供自定义调色板。无论如何，这对我来说都是一个技巧，但我确信它可以被整理和优化。verbose = T输出了很多输出，我用它来调试功能，当它中断时。

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在没有选项的情况下执行相同的任务，但在现有的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.

所以我的方法是，试着用你所有的数据(我很容易地绘制出从~15k点生成的曲面)。如果这种方法不起作用，可以使用几个较小的样本，并在使用这些函数时将它们全部绘制出来。

#4

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.

你可以用晶格。在这个例子中，我定义了一个网格，我想画出z~x,y。它看起来像这样。注意，大部分代码只是构建一个三维形状，我用线框函数绘图。

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

变量b和s可以是x或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))

#5

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

也许现在已经很晚了，但是跟随Spacedman，你试过重复的="strip"还是其他的选择?

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)