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

Question

I'm trying to make a figure similar to this in R. Basically X and Y represent coordinates (e.g. lon/lat) and Z represents elevation. I want to plot a 3D bar chart with a elevation as the bar height, and then show a smooth curve of fitted Z values going through the bars. I am able to get a 3d barplot with barplot3d but I'm not sure if it's possible to add the fitted curve on top of that. Does anyone know how to do this? I have some example code below demonstrating what I've tried so far.

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)

Update: More complicated example

The curve doesn't seem to intersect the bars properly with current solution for a non-symmetrical curve (fixed in most updated answer).

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()

Answer 1

As you can probably tell, your 3d surface is rotated 90 degrees relative to where it should be. This is not your fault; it is just a difference between the way barplot3d is drawn compared to the other rgl shapes. You also need a bit of shifting and rescaling to get it to fit.

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)

---

Update

To remove the points above the highest bar, just set them to NA. You probably want to increase the resolution when you do this though:

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)

Note that this gives an unavoidably ragged edge to the graphic

An alternative is to shrink the x, y grids at which z is calculated:

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)

---

Further update

It looks as though we need a double flip to get the z values correct:

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()