R: построение трехмерной поверхности по x, y, z

Question

представьте, у меня есть матрица из 3 столбцов
х, у, г где z является функцией от x и y.

Я знаю, как построить «график рассеяния» этих точек с plot3d(x,y,z)

Но если я хочу вместо поверхности, я должен использовать другие команды, такие как surface3d Проблема в том, что он не принимает те же входные данные, что plot3d похоже, нужна матрица с

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

Как я могу получить эту матрицу? Я пробовал с помощью команды interp, как я делаю, когда мне нужно использовать контурные графики.

Как я могу построить поверхность непосредственно из x, y, z без вычисления этой матрицы? Если бы у меня было слишком много очков, эта матрица была бы слишком большой.

ура

Answer 1

Если ваши координаты x и y не находятся на сетке, вам нужно интерполировать поверхность x, y, z на одну.Вы можете сделать это с помощью кригинга, используя любой из пакетов геостатистики (geoR, gstat и т. Д.) Или более простые методы, такие как взвешивание по обратному расстоянию.

Я предполагаю, что упоминаемая вами функция 'interp' из пакета akima,Обратите внимание, что выходная матрица не зависит от размера ваших входных точек.Вы можете иметь 10000 точек на входе и интерполировать их на сетку 10х10, если хотите.По умолчанию akima :: interp делает это на сетке 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)

Это будет выглядеть колючим и мусором из-за случайных данных.Надеюсь, ваши данные не!

Answer 2

Вы можете посмотреть на использование решетки. В этом примере я определил сетку, по которой я хочу построить z ~ x, y. Это выглядит примерно так. Обратите внимание, что большая часть кода просто создает трехмерную фигуру, которую я строю, используя функцию каркаса.

Переменные "b" и "s" могут быть x или 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 3

Вы можете использовать функцию outer() для ее генерации.

Взгляните на демонстрацию для функции persp(), которая является базовой графической функцией для рисования перспективных графиков для поверхностей.

Вот их первый пример:

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)

То же самое относится к surface3d():

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

Answer 4

rgl отлично, но нужно немного поэкспериментировать, чтобы получить правильные оси.

Если у вас много точек, почему бы не взять случайную выборку из них, а затем построить результирующую поверхность.Вы можете добавить несколько поверхностей на основе выборок из одних и тех же данных, чтобы увидеть, как процесс выборки ужасно влияет на ваши данные.

Итак, вот довольно ужасная функция, но она делает то, что, я думаю, вы хотите, чтобы она делала (но без выборки).Учитывая матрицу (x, y, z), где z - высоты, она построит как точки, так и поверхность.Ограничения заключаются в том, что для каждой пары (x, y) может быть только один z.Таким образом, плоскости, которые зацикливаются на самих себе, вызовут проблемы.

plot_points = T построит отдельные точки, из которых сделана поверхность - это полезно для проверки того, что поверхность и точки действительно совпадают.plot_contour = T отобразит двухмерный контур ниже трехмерной визуализации.Установите цвет на rainbow, чтобы получить красивые цвета, все остальное будет серым, но затем вы можете изменить функцию, чтобы получить собственную палитру.В любом случае, это помогает мне, но я уверен, что это можно привести в порядок и оптимизировать.verbose = T выводит много выходных данных, которые я использую для отладки функции по мере ее выхода из строя.

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 выполняет ту же работу без параметров, но накладывает поверхность на существующий 3d-график.

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

Итак, мой подход заключается в том, чтобы попытаться сделать это со всеми вашими данными (я легко рисую поверхности, сгенерированные из ~ 15k точек).Если это не сработает, возьмите несколько небольших выборок и нанесите их все сразу, используя эти функции.

Answer 5

Может быть, уже поздно, но, следуя Spacedman, вы пытались использовать duplicate = "strip" или какой-либо другой вариант?

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)