R: Vẽ một bề mặt 3D từ x, y, z

Question

tưởng tượng tôi có ma trận 3 cột
x, y, z trong đó z là hàm của x và y.

tôi biết làm thế nào để vẽ một "âm mưu phân tán" của những điểm này với plot3d(x,y,z)

Nhưng nếu tôi muốn có một bề mặt thay vào đó tôi phải sử dụng các lệnh khác như surface3d Vấn đề là nó không chấp nhận các đầu vào tương tự như plot3d có vẻ như cần một ma trận với

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

Làm cách nào để có được ma trận này? Tôi đã thử với lệnh interp, như tôi làm khi tôi cần sử dụng các ô đường viền.

Làm cách nào để vẽ một bề mặt trực tiếp từ x, y, z mà không tính toán ma trận này? Nếu tôi có quá nhiều điểm ma trận này sẽ là quá lớn.

cổ vũ

Answer 1

Bạn có thể sử dụng chức năng outer() để tạo.

Hãy xem bản trình diễn cho hàm persp(), là một chức năng đồ họa cơ sở để vẽ các ô phối cảnh cho các bề mặt.

Dưới đây là ví dụ đầu tiên của họ:

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) 

cũng áp dụng như surface3d():

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

Answer 2

Nếu x và coords y là không phải trên một lưới thì bạn cần phải suy x của bạn, y , z bề mặt lên một. Bạn có thể làm điều này với kriging bằng cách sử dụng bất kỳ gói địa lý (geoR, gstat, những người khác) hoặc các kỹ thuật đơn giản hơn như trọng số khoảng cách nghịch đảo.

Tôi đoán hàm 'interp' mà bạn đề cập đến là từ gói akima. Lưu ý rằng ma trận đầu ra là độc lập với kích thước của các điểm đầu vào của bạn. Bạn có thể có 10000 điểm trong đầu vào của bạn và nội suy mà vào một lưới 10x10 nếu bạn muốn. Theo mặc định akima :: interp thực hiện trên lưới 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) 

Điều đó sẽ có vẻ đáng sợ và rác vì dữ liệu ngẫu nhiên của nó. Hy vọng rằng dữ liệu của bạn là không!

Answer 3

rgl thật tuyệt, nhưng phải mất một chút thời gian thử nghiệm để có được trục đúng.

Nếu bạn có nhiều điểm, tại sao không lấy mẫu ngẫu nhiên từ chúng, và sau đó vẽ bề mặt kết quả. Bạn có thể thêm một số bề mặt tất cả dựa trên mẫu từ cùng một dữ liệu để xem liệu quá trình lấy mẫu có ảnh hưởng khủng khiếp đến dữ liệu của bạn hay không.

Vì vậy, đây là một chức năng khá khủng khiếp nhưng nó làm những gì tôi nghĩ rằng bạn muốn nó làm (nhưng không có lấy mẫu). Cho một ma trận (x, y, z) trong đó z là chiều cao nó sẽ vẽ cả hai điểm và cũng là một bề mặt. Hạn chế là chỉ có thể có một z cho mỗi cặp (x, y). Vì vậy, máy bay lặp lại chính mình sẽ gây ra vấn đề.

plot_points = T sẽ vẽ các điểm riêng lẻ mà từ đó bề mặt được tạo ra - điều này rất hữu ích để kiểm tra xem bề mặt và các điểm có thực sự gặp nhau hay không. plot_contour = T sẽ vẽ đồ thị đường bao 2d bên dưới hình ảnh 3D. Đặt màu thành rainbow để cho màu sắc đẹp, bất kỳ thứ gì khác sẽ đặt màu xám, nhưng sau đó bạn có thể thay đổi hàm để cung cấp bảng màu tùy chỉnh. Tuy nhiên, điều này có thể giúp tôi nhưng tôi chắc chắn rằng nó có thể được dọn dẹp và tối ưu hóa. Các verbose = T in ra rất nhiều đầu ra mà tôi sử dụng để gỡ lỗi các chức năng như và khi nó phá vỡ.

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 thực hiện cùng một công việc mà không có tùy chọn, nhưng phủ lên bề mặt trên vùng hiện có 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") 
} 

Vì vậy, cách tiếp cận của tôi sẽ là, cố gắng làm điều đó với tất cả dữ liệu của bạn (Tôi dễ dàng vẽ các bề mặt được tạo từ ~ 15k điểm). Nếu điều đó không hiệu quả, hãy lấy một vài mẫu nhỏ hơn và vẽ tất cả chúng cùng một lúc bằng các hàm này.

Answer 4

Bạn có thể xem bằng cách sử dụng mạng. Trong ví dụ này, tôi đã xác định một lưới mà tôi muốn vẽ đồ thị z ~ x, y. Nó trông như thế này Lưu ý rằng hầu hết các mã chỉ là xây dựng một hình dạng 3D mà tôi vẽ bằng cách sử dụng chức năng wireframe.

Biến "b" và "s" có thể là x hoặc 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

Có thể là muộn nhưng sau Spacedman, bạn có thử sao chép = "strip" hoặc bất kỳ tùy chọn nào khác không?

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)