2010-10-20 6 views
28

Immagina di avere una matrice a 3 colonne
x, y, z dove z è una funzione di x e y.R: Tracciare una superficie 3D da x, y, z

so come tracciare un "grafico a dispersione" di questi punti con plot3d(x,y,z)

Ma se voglio una superficie invece devo usare altri comandi come surface3d Il problema è che non accetta il stessi ingressi come tracciare3d sembra avere bisogno di una matrice con

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

Come posso ottenere questa matrice? Ho provato con il comando interp, come faccio quando devo usare i grafici del profilo.

Come è possibile tracciare una superficie direttamente da x, y, z senza calcolare questa matrice? Se avessi troppi punti questa matrice sarebbe troppo grande.

applausi

+1

da dove proviene plot3d ??? qual è il pacchetto che stai usando? –

+0

rgl, ma potrei usare qualsiasi suggerimento – skan

risposta

6

È possibile utilizzare la funzione outer() per generarlo.

Dai un'occhiata alla demo per la funzione persp(), che è una funzione grafica di base per disegnare grafici prospettici per le superfici.

Ecco il loro primo esempio:

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 stesso vale per surface3d():

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

Ciao. Non ho una funzione esplicita. I miei dati vengono esportati da un altro programma, provengono da una complessa ottimizzazione, potremmo dire che provengono da una scatola nera. Ho x, y e z ma non ho alcun mezzo per calcolare altri valori di z per le altre combinazioni x e y. – skan

26

Se il x e coordinate Y non sono su una griglia allora avete bisogno di interpolare il tuo x, y , z superficie su uno. Puoi farlo con krig usando uno dei pacchetti di geostatistica (geoR, gstat, altri) o tecniche più semplici come la pesatura a distanza inversa.

Sto indovinando la funzione "interp" che si menziona dal pacchetto akima. Nota che la matrice di output è indipendente dalla dimensione dei tuoi punti di input. Potresti avere 10000 punti nel tuo input e interpolarlo su una griglia 10x10 se lo desideri. Di Akima impostazione predefinita :: interp lo fa su una griglia 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) 

Che ti aspetto spinoso e spazzatura perché i suoi dati casuali. Spero che i tuoi dati non siano!

+1

Ciao. Questo è quello che ho fatto a volte ma a volte non funziona. Esempio di Fow nel tuo esempio Ricevo molte NA quando uso interp. – skan

+2

Il problema è che interp non funziona se i miei punti sono collineari – skan

+1

Come potrei farlo con misc3d parametric3d o altra funzione? – skan

5

rgl è ottimo, ma richiede un po 'di sperimentazione per ottenere gli assi giusti.

Se si dispone di molti punti, perché non prelevarne un campione casuale, quindi tracciare la superficie risultante. Puoi aggiungere diverse superfici tutte basate su campioni dagli stessi dati per vedere se il processo di campionamento ha un impatto terribile sui tuoi dati.

Quindi, ecco una funzione piuttosto orribile ma fa quello che penso tu voglia che faccia (ma senza il campionamento). Data una matrice (x, y, z) dove z è l'altezza, traccerà sia i punti che una superficie. Le limitazioni sono che ci può essere solo una z per ciascuna coppia (x, y). Quindi gli aerei che si ricollocano su se stessi causeranno problemi.

Il plot_points = T traccia i singoli punti da cui viene creata la superficie, utile per verificare che la superficie e i punti effettivamente si incontrino. Lo plot_contour = T traccerà un grafico di contorno 2d sotto la visualizzazione 3d. Imposta il colore su rainbow per ottenere dei bei colori, qualsiasi altra cosa lo imposterà in grigio, ma puoi modificare la funzione per dare una tavolozza personalizzata. Questo è il trucco per me comunque, ma sono sicuro che può essere riordinato e ottimizzato. Il verbose = T stampa un sacco di output che uso per eseguire il debug della funzione come e quando si interrompe.

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

Il add_rgl_model fa lo stesso lavoro senza le opzioni, ma sovrappone una superficie sulla 3dplot esistente.

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

Quindi il mio approccio sarebbe quello di, provare a farlo con tutti i tuoi dati (ho facilmente tracciare superfici generate da ~ 15k punti). Se ciò non funziona, prendi diversi campioni più piccoli e tracciali tutti insieme usando queste funzioni.

5

Si potrebbe usare Lattice. In questo esempio ho definito una griglia sulla quale voglio tracciare z ~ x, y. Sembra qualcosa del genere. Si noti che la maggior parte del codice è solo la creazione di una forma 3D che viene tracciata utilizzando la funzione wireframe.

Le variabili "b" e "s" potrebbero essere 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)) 
3

Forse è tardi ma seguendo Spacedman, hai provato duplicato = "strip" o qualsiasi altra opzione?

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)