Add a google map at x,y,0 to a 3d map

I’m using the R plot_ly library. I have generated a 3d surface of ozone concentrations in my neighborhood (Grand Junction, CO). I’d like to put a map of the region ‘underneath’ it so the casual observer will have some feel for where the data are located in the real world.

Here is the code for the 3d surface

x <- c("-108.6125","-108.5114","-108.805","-108.4014","-108.5615","-108.8349","-108.225","-108.3139","-108.5568","-108.4968")
y <- c(“39.02205”,“39.22255”,“39.598”,“38.89478”,“39.06429”,“39.27625”,“39.03”,“39.1306”,“39.14823”,“38.89795”)
z <- c(“60.7735”,“56.45783”,“49.65”,“60.15”,“50”,“53.95417”,“50.825”,“56”,“55.843”,“38.73333”)
df <- data.frame(x = as.numeric(x),y = as.numeric(y),z = as.numeric(z))

s = interp(x = df$x, y = df$y, z = df$z)
p <- plot_ly(x = s$x, y = s$y, z = s$z) %>% add_surface()

And some code that gets a google map that is relatively close (the ggmap code doesn’t generate quite what I want yet but it’s close)

get google map of the grand valley

loc <- c(-109, 38.8, -108, 39.8)
tx_map_gmaps <- get_map(location=loc, zoom = 9, source=“google”, maptype=“terrain”)
gg <- ggmap(tx_map_gmaps)

Is there some way to put the google map at x,y,0?

You can plot the map as a scatter3d plot. Create xx and yy arrays with the same dimensions as gg, and with the correct geocoordinates for the map you have retrieved. Then use add_trace(p, x=xx, y=yy, z=0, type=“scatter3d”, color=gg). If you are mapping a large area, gg will be large and you may need to subsample a lower-resolution version to get it to plot correctly.

Here’s my code (I’ve used “landmap” instead of “gg” and “MapLocation” instead of “loc”, and have used “ind1” and “ind2” to reduce the resolution of my map):

landmap <- get_map(location=MapLocation, crop=TRUE)
bb <- as.numeric(attributes(landmap)$bb)
landy <- seq(bb[3], bb[1], length.out = dim(landmap)[1])
landx <- seq(bb[2], bb[4], length.out = dim(landmap)[2])
landxx <- (array(landx, dim=dim(t(landmap))))
landyy <- t(array(landy, dim=dim(landmap)))

ind1 <- seq(1, dim(landxx)[1], 4)
ind2 <- seq(1, dim(landxx)[2], 4)

p <- add_trace(p,  x=(c(landxx[ind1, ind2])), 
	       y=(c(landyy[ind1, ind2])), 
               marker=list(color=c(landmap[ind1, ind2]), size=4),