绘制网格的3D地形图

Question

我想使用R语言来制作类似以下图片的图形，可能需要使用ggplot（尽管我怀疑这不可能，因为据我所知它没有三维能力）。我通常从raster包中获取栅格文件数据，但我可以将其转换为最适合的格式。 这张图取自于"Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change"，第1.14图。我不知道是哪个软件生成了这张图。

我猜想唯一的方法是使用lattice::cloud (wireframe) 或类似的东西？我似乎找不到任何办法来强制wireframe具有柱状图而不是表面图的属性；此外，基于海拔高度进行着色，在保持网格水平的情况下，可能是不可能的。

- AF7

---

也许你可以在这个答案的基础上进行构建：https://dev59.com/FYLba4cB1Zd3GeqPYgK_#25042436。它很灵活，允许你使用单独的颜色映射。 - koekenbakker

@koekenbakker 谢谢。我已经找到了一个可能的解决方案（请参见下面的我的答案），但对于非常大的图表来说，它可能会很慢（我确实是指非常大的图表）。我会尝试抽出一些时间来看看那两个答案。 - AF7

感谢。最终我找到了一个解决方案（感谢Ben Bolker），请见下面的答案。 - AF7

你为什么想要制作这样的图形？颜色和高度似乎编码了相同的数据，因此其中一个是多余的，并且3D性质意味着一些正方形被遮挡。只需进行平面彩色光栅绘图，并暂时放弃眼花缭乱的效果。 - Spacedman

@Spacedman 我也做了一个平面光栅图。我有不同分辨率的数据集，直接比较3D图形可以更容易地理解差异。只需查看《气候变化2013：物理科学基础——第五次评估报告工作组I对政府间气候变化专门委员会的贡献》，图1.14，你就会明白我的意思。此外，我想做可旋转的3D图形，这样可以放大并观察所选区域的差异。然而，最终目的主要是为了学习一些新的东西，以备将来有用。 - AF7

Answer 1

我找到了一些时间来研究这个问题。最终我采取了以下方法（链接的R代码）： http://pastebin.com/dA2nNNS0 感谢评论中指导我的人。

但我仍然有保存文件的问题，这是另一个问题的材料。

来自pastebin的代码：

library(raster) 
library(rgl)

r <- raster("topo12.nc", varname="topo") #Read file
cols <- terrain.colors(3100) #Define colors

binplot.3d <- function(x,y,z,alpha=1,topcol="#ff0000",sidecol="#aaaaaa"){ #Binplotting function
  save <- par3d(skipRedraw=TRUE)
  on.exit(par3d(save))

  x1<-c(rep(c(x[1],x[2],x[2],x[1]),3),rep(x[1],4),rep(x[2],4))
  z1<-c(rep(0,4),rep(c(0,0,z,z),4))
  y1<-c(y[1],y[1],y[2],y[2],rep(y[1],4),rep(y[2],4),rep(c(y[1],y[2],y[2],y[1]),2))
  x2<-c(rep(c(x[1],x[1],x[2],x[2]),2),rep(c(x[1],x[2],rep(x[1],3),rep(x[2],3)),2))
  z2<-c(rep(c(0,z),4),rep(0,8),rep(z,8) )
  y2<-c(rep(y[1],4),rep(y[2],4),rep(c(rep(y[1],3),rep(y[2],3),y[1],y[2]),2) )
  rgl.quads(x1,z1,y1,col=rep(sidecol,each=4),alpha=alpha)
  rgl.quads(c(x[1],x[2],x[2],x[1]),rep(z,4),c(y[1],y[1],y[2],y[2]),
              col=rep(topcol,each=4),alpha=1) 
  rgl.lines(x2,z2,y2,col="#000000")
}

cat("Row ( of", dim(r)[1],"):")
for (row in 1:dim(r)[1]) { #Plotting loop
    for (col in 1:dim(r)[2]) {
    if (round(r[row, col]) < 1) {
    binplot.3d(c(col-1,col), c(row-1,row), r[row, col]/500, alpha=1, topcol="cadetblue3")
    } else {
    binplot.3d(c(col-1,col), c(row-1,row), r[row, col]/500, alpha=1, topcol=cols[round(r[row, col])])
#   cat(round(r[row, col]), "\t")
    }
    }
    cat(row, "")
}

- AF7