R: Tracer une surface 3D à partir de x, y, z

Question

Imaginez que j'ai une matrice de 3 colonnes
x, y, z où z est une fonction de x et y.

Je sais comment tracer un « diagramme de dispersion » de ces points avec plot3d(x,y,z)

Mais si je veux une surface plutôt que je dois utiliser d'autres commandes telles que surface3d Le problème est qu'il n'accepte pas la mêmes entrées que plot3d il semble avoir besoin d'une matrice avec

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

Comment puis-je obtenir cette matrice? J'ai essayé avec la commande interp, comme je le fais quand j'ai besoin d'utiliser des courbes de niveau.

Comment puis-je tracer une surface directement à partir de x, y, z sans calculer cette matrice? Si j'avais trop de points, cette matrice serait trop grande.

acclamations

Answer 1

Vous pouvez utiliser la fonction outer() pour générer. Jetez un oeil à la démo pour la fonction persp(), qui est une fonction graphique de base pour dessiner des graphiques en perspective pour les surfaces.

Voici leur premier exemple:

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) 

de même pour surface3d():

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

Answer 2

Si votre x et y coords ne sont pas sur une grille alors vous devez interpoler votre x, y , z surface sur un. Vous pouvez le faire avec kriging en utilisant l'un des paquets de géostatistique (geoR, gstat, autres) ou des techniques plus simples telles que la pondération de distance inverse.

Je suppose que la fonction 'interp' dont vous parlez provient du paquet akima. Notez que la matrice de sortie est indépendante de la taille de vos points d'entrée. Vous pourriez avoir 10000 points dans votre entrée et l'interpoler sur une grille de 10x10 si vous le vouliez. Par défaut Akima :: interp il ne sur une grille 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) 

Ça va regarder hérissés et des déchets parce que ses données aléatoires. J'espère que vos données ne sont pas!

Answer 3

rgl est génial, mais il faut un peu d'expérimentation pour obtenir les bons axes.

Si vous avez beaucoup de points, pourquoi ne pas en prélever un échantillon aléatoire, puis tracer la surface résultante. Vous pouvez ajouter plusieurs surfaces basées sur des échantillons provenant des mêmes données pour voir si le processus d'échantillonnage affecte horriblement vos données. Donc, voici une fonction assez horrible mais elle fait ce que je pense que vous voulez faire (mais sans échantillonnage). Étant donné une matrice (x, y, z) où z représente les hauteurs, elle trace à la fois les points et aussi une surface. Les limites sont qu'il ne peut y avoir qu'un z pour chaque paire (x, y). Les avions qui se replient sur eux-mêmes causeront des problèmes.

Le plot_points = T tracera les points individuels à partir desquels la surface est faite - ceci est utile pour vérifier que la surface et les points se rencontrent effectivement. Le plot_contour = T tracera un tracé de contour 2D sous la visualisation 3D. Définir la couleur à rainbow pour donner de jolies couleurs, tout le reste le mettra en gris, mais vous pouvez ensuite modifier la fonction pour donner une palette personnalisée. Cela fait l'affaire pour moi de toute façon, mais je suis sûr qu'il peut être rangé et optimisé. Le verbose = T imprime beaucoup de sortie que j'utilise pour déboguer la fonction au fur et à mesure qu'elle se casse.

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

Le add_rgl_model fait le même travail sans les options, mais recouvre une surface sur la 3DPlot existante.

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

Alors mon approche serait, essayer de le faire avec toutes vos données (je tracer facilement les surfaces générées par ~ 15k points). Si cela ne fonctionne pas, prenez plusieurs petits échantillons et tracez-les tous en même temps en utilisant ces fonctions.

Answer 4

Vous pourriez envisager d'utiliser Lattice. Dans cet exemple, j'ai défini une grille sur laquelle je veux tracer z ~ x, y. Ça ressemble à quelque chose comme ça. Notez que la majeure partie du code consiste simplement à créer une forme 3D que je tracerai en utilisant la fonction filaire.

Les variables "b" et "s" peuvent être x ou 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

Peut-être est-il en retard maintenant mais après Spacedman, avez-vous essayé duplicate = "strip" ou toute autre 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)