para - poligono en r

R: trazando una superficie 3D de x, y, z

Question

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

Sé cómo trazar un "diagrama de dispersión" de estos puntos con plot3d(x,y,z)

Pero si quiero una superficie en su lugar, debo usar otros comandos, como surface3d. El problema es que no acepta las mismas entradas que plot3d, parece que necesita 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.

aclamaciones

Answer 1

Podrías mirar usando 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 2

Puede usar 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í está 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 3

Si sus coordenadas xey no están en una cuadrícula, entonces necesita interpolar su superficie x, y, z en una. 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 grilla 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 se verá puntiagudo y basura porque sus datos son aleatorios. ¡Espero que tus datos no!

Answer 4

Tal vez es tarde ahora, pero después de 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)

Answer 5

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

Si tiene muchos puntos, ¿por qué no tomar una muestra aleatoria de ellos y luego graficar 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.

Entonces, esta es una función bastante horrible, pero hace lo que creo que quieres 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 los puntos individuales a partir de los cuales se plot_points = T la superficie; esto es útil para verificar que la superficie y los puntos se encuentren realmente. El plot_contour = T un gráfico de contornos 2d debajo de la visualización en 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 como y 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 realiza el mismo trabajo sin las opciones, pero superpone una superficie en el 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 intentar, tratar de hacerlo con todos sus datos (trazar fácilmente las 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.