R：从x，y，z绘制三维曲面
想象我有一个3列matrix
x，y，z其中z是x和y的函数。
我知道如何用plot3d(x,y,z)绘制这些点的“散点图”
但是，如果我想要一个表面，而不是我必须使用其他命令如surface3d问题是，它不接受相同的input作为plot3d它似乎需要一个matrix与
(nº elements of z) = (n of elements of x) * (n of elements of x)
我怎样才能得到这个matrix？ 我已经尝试了interp命令，就像我需要使用轮廓图一样。
如何直接从x，y，z绘制曲面而不计算这个matrix？ 如果我的分数太多，这个matrix就太大了。
干杯
如果你的x和y坐标不在网格上，那么你需要将你的x，y，z曲面插值到一个上。 您可以使用任何地理统计软件包（geoR，gstat，其他）或简单的技术（如反距离加权）进行克里格处理。
我猜你所提到的“interp”function是来自akima软件包。 请注意，输出matrix与input点的大小无关。 你可以在你的input中有10000个点，并且如果你想要的话把它内插到一个10×10的网格上。 默认情况下，akima :: interp会将其转换为40×40网格：
> 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)
这将看起来尖刻和垃圾，因为它的随机数据。 希望你的数据不是！
你可以看看使用莱迪思。 在这个例子中，我已经定义了一个我想要绘制z〜x，y的网格。 看起来像这样 请注意，大部分代码只是build立一个3Dgraphics，我使用线框function进行绘图。
variables“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))
你可以使用函数outer()来生成它。
看一看函数persp()的演示，这是一个绘制曲面透视图的基本graphics函数。
这是他们的第一个例子：
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)
这同样适用于surface3d() ：
require(rgl) surface3d(x, y, z)
rgl很棒，但需要一些实验才能正确地获取轴。
如果你有很多的分数，为什么不从他们那里随机抽样，然后绘制结果表面。 您可以根据来自相同数据的样本添加多个曲面，以查看采样过程是否会严重影响您的数据。
所以，这是一个非常可怕的function，但它做我认为你想要做的（但没有采样）。 给定一个matrix（x，y，z），其中z是高度，它将绘制点和曲面。 限制是每个（x，y）对只能有一个z。 所以飞回自己的飞机会造成问题。
plot_points = T将绘制曲面的各个点 – 这对检查曲面和点实际是否plot_points = T很有用。 plot_contour = T将绘制三维可视化下面的二维等高线图。 设置颜色为rainbow给漂亮的颜色，其他任何东西都会把它设置为灰色，但是你可以改变function给一个自定义的调色板。 无论如何，这对我来说是个诀窍，但我相信它可以被整理和优化。 verbose = T打印出大量的输出，我用它来debugging函数，当它打破。
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) } }
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") }
所以我的方法是，试着用你所有的数据来做（我可以轻松地绘制从〜15k点生成的曲面）。 如果这样做不起作用，则使用这些函数将一些较小的样本一次全部绘制出来。
也许现在已经晚了，但是在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)
