R: Trazar una superficie 3D a partir de x, y, z

Question

imagine que tengo una matriz de 3 columnas
x, y, z donde z es una función de xey.

sé cómo trazar una "nube de puntos" de estos puntos con plot3d(x,y,z)

Pero si quiero una superficie en lugar de eso hay que utilizar otros comandos como Surface3D El problema es que no acepta la mismas entradas que dibujar3d que parece necesitar una matriz con

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

¿Cómo puedo obtener esta matriz? Lo he intentado con el comando interp, como hago cuando necesito usar diagramas de contorno.

¿Cómo puedo trazar una superficie directamente desde x, y, z sin calcular esta matriz? Si tuviera demasiados puntos, esta matriz sería demasiado grande.

aplausos

Answer 1

Puede utilizar la función outer() para generarlo.

Eche un vistazo a la demostración de la función persp(), que es una función gráfica básica para dibujar gráficos de perspectiva para superficies.

Aquí es su primer ejemplo:

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) 

Lo mismo se aplica a surface3d():

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

Answer 2

Si su xy coordenadas Y no sean en una cuadrícula, entonces usted necesita para interpolar su x, y , z superficie en uno. Puede hacer esto con kriging usando cualquiera de los paquetes de geoestadística (geoR, gstat, otros) o técnicas más simples como la ponderación de distancia inversa.

Supongo que la función 'interp' que mencionas proviene del paquete akima. Tenga en cuenta que la matriz de salida es independiente del tamaño de sus puntos de entrada. Podría tener 10000 puntos en su entrada e interpolarlos en una cuadrícula de 10x10 si quisiera. Por defecto Akima :: interp lo hace en una rejilla de 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) 

Eso va a mirar de punta y la basura debido a que sus datos aleatorios. ¡Espero que tus datos no!

Answer 3

rgl es genial, pero requiere un poco de experimentación para obtener los ejes correctos.

Si tiene muchos puntos, ¿por qué no toma una muestra al azar de ellos y luego grafica la superficie resultante? Puede agregar varias superficies, todas basadas en muestras de los mismos datos, para ver si el proceso de muestreo está afectando horriblemente sus datos.

Por lo tanto, esta es una función bastante horrible pero hace lo que creo que quiere que haga (pero sin el muestreo). Dada una matriz (x, y, z) donde z es la altura, trazará los puntos y también una superficie. Las limitaciones son que solo puede haber una z para cada par (x, y). Entonces los aviones que vuelven sobre sí mismos causarán problemas.

El plot_points = T trazará los puntos individuales desde los que se fabrica la superficie; esto es útil para verificar que la superficie y los puntos se encuentren realmente. El plot_contour = T trazará un gráfico de contorno 2D debajo de la visualización 3D. Establezca el color en rainbow para dar colores bonitos, cualquier otra cosa lo configurará en gris, pero luego puede modificar la función para obtener una paleta personalizada. Esto me funciona de todos modos, pero estoy seguro de que puede ser arreglado y optimizado. El verbose = T imprime una gran cantidad de salida que utilizo para depurar la función cuando se rompe.

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) 
    } 
} 

El add_rgl_model hace el mismo trabajo sin las opciones, pero se superpone a una superficie sobre la 3DPlot existente.

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") 
} 

Así que mi enfoque sería, tratar de hacerlo con todos sus datos (I trazo fácilmente superficies generadas a partir de ~ 15k puntos). Si eso no funciona, tome varias muestras más pequeñas y grábelas todas a la vez usando estas funciones.

Answer 4

Puede consultar el uso de Lattice. En este ejemplo he definido una grilla sobre la cual quiero trazar z ~ x, y. Se ve algo como esto. Tenga en cuenta que la mayor parte del código es simplemente la construcción de una forma 3D que trazo utilizando la función de estructura de alambre.

Las variables "b" y "s" podrían ser x o 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)) 

Answer 5

Quizás es tarde ahora pero siguiendo Spacedman, ¿has probado duplicate = "strip" o alguna otra opción?

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)