如何在 R 中创建带有拟合曲线的 3D 条形图

How to create a 3D barplot with fitted curve in R

我正在尝试在 R 中制作类似于 this 的图形。基本上 X 和 Y 代表坐标(例如 lon/lat),Z 代表海拔。我想绘制一个以高程作为条高的 3D 条形图,然后显示穿过条的拟合 Z 值的平滑曲线。我能够使用 barplot3d 获得 3d 条形图,但我不确定是否可以在其上添加拟合曲线。有谁知道如何做到这一点?我在下面有一些示例代码展示了我到目前为止所做的尝试。

library(rgl)
library(barplot3d)
library(tidyverse)

x_mat<- matrix(rep(-1:1, each=3),nrow = 3) #x coordinates
y_mat<- matrix(rev(rep(-1:1, each=3)),nrow = 3, byrow = TRUE) #y coordinates

df<- data.frame(x = as.vector(x_mat), y = as.vector(y_mat)) #dataframe
set.seed(5)
df<- df %>% mutate(z= x^2+ y^2 + rnorm(n = 9, mean = 0, sd = 0.1)) #add elevation values

m<- lm(z ~ I(x^2)+I(y^2)+I(x*y)+x+y, data = df) #fitted curve

rgl.open()
rgl::plot3d(m)
barplot3d(rows=3,cols=3, z=df$z,scalexy=1, gap=0, alpha=0.4,theta=30,phi=50,
          topcolors = "gray", gridlines = TRUE)

更新:更复杂的例子

曲线似乎没有与非对称曲线的当前解决方案正确相交(在最新的答案中已修复)。

library(rgl)
library(barplot3d)
library(tidyverse)

x_mat<- matrix(rep(-1:1, each=3),nrow = 3) #x coordinates
y_mat<- matrix(rev(rep(-1:1, each=3)),nrow = 3, byrow = TRUE) #y coordinates
A<- 0.2
B<- 0.2
C<- 0.4
D<- 0.4
E<- 0

df<- data.frame(x = as.vector(x_mat), y = as.vector(y_mat)) #dataframe
df<- df %>% mutate(z= A*x^2 + B*y^2 + C*x*y + D*x + E*y) #add elevation values
z_mat<- matrix(data = df$z, nrow=3)

m<- lm(z ~ I(x^2)+I(y^2)+I(x*y)+x+y, data = df) #fitted curve
df$zpred<- predict(m, data.frame(x=df$x, y=df$y))
round(df$z-df$zpred,10) #Predictions should fit observations almost exactly (i.e. intersect exactly with bars)
#0 0 0 0 0 0 0 0 0 

n<- 10
xvals <- seq(-1, 1, len = n)
xmat <- replicate(n, seq(1.5, 3.5, len = n))
ymat <- t(xmat)
pred <- expand.grid(x = xvals, y = xvals)
zmat <- matrix(predict(m, pred), nrow = n, ncol = n)

barplot3d(rows=3,cols=3, z=df$z, gap=0, alpha=0.4, phi = 45,
          topcolors = "gray", sidecolors = "cyan", linecolors= "blue", gridlines = FALSE, zlabels = FALSE)
surface3d(x = xmat, y = zmat,  z = ymat-5, color = "purple", alpha = 0.7)
axes3d()

正如您可能知道的那样,您的 3d 表面相对于其应有的位置旋转了 90 度。这不是你的错;这只是 barplot3d 的绘制方式与其他 rgl 形状的区别。您还需要进行一些移动和重新缩放以使其适合。

barplot3d(rows=3,cols=3, z=df$z, gap=0, alpha=0.4, phi = 45,
          topcolors = "gray", gridlines = TRUE)

xvals <- seq(-1.5, 1.5, len = 10)
xmat <- replicate(10, seq(1, 4, len = 10))
ymat <- t(xmat)
pred <- expand.grid(x = xvals, y = xvals)
zmat <- matrix(predict(m, pred), 10, 10)

surface3d(xmat, zmat, color = "gold", alpha = 0.5, ymat - 5)


更新

要删除最高条上方的点,只需将它们设置为 NA。不过,您可能希望在执行此操作时提高分辨率:

xvals <- seq(-1.5, 1.5, len = 100)
xmat <- replicate(100, seq(1, 4, len = 100))
ymat <- t(xmat) - 5

pred <- expand.grid(x = xvals, y = xvals)
zmat <- matrix(predict(m, pred), 100, 100)

zmat[zmat > 2] <- NA

barplot3d(rows=3,cols=3, z=df$z, gap=0, alpha=0.4, phi = 45,
          topcolors = "gray", gridlines = TRUE)

surface3d(xmat, zmat, color = "gold", alpha = 0.5, ymat)

请注意,这会给图形带来不可避免的参差不齐的边缘

另一种方法是缩小计算 z 的 x、y 网格:

xvals <- seq(-1, 1, len = 100)
xmat <- replicate(100, seq(1.5, 3.5, len = 100))
ymat <- t(xmat) - 5

pred <- expand.grid(x = xvals, y = xvals)
zmat <- matrix(predict(m, pred), 100, 100)

barplot3d(rows=3,cols=3, z=df$z, gap=0, alpha=0.4, phi = 45,
          topcolors = "gray", gridlines = TRUE)

surface3d(xmat, zmat, ymat, color = "gray20", alpha = 0.5)


进一步更新

看来我们需要双击才能使 z 值正确:

barplot3d(rows=3,cols=3, z=df$z, gap=0, alpha=0.4, phi = 45,
          topcolors = "gray", sidecolors = "cyan", linecolors= "blue", 
          gridlines = FALSE, zlabels = FALSE)
surface3d(x = xmat, y = t(apply(t(apply(zmat, 1, rev)), 2, rev)),  
          z = ymat-5, color = "purple", alpha = 0.7)
axes3d()