Plot a 3D plume with a map underneath

In order to make an impressive plume plot of my oversampling work, I would like to plot a 3D plume (where hight represents the vertical  column density) with a map underneath. It turns out the RGL package, a powerful tool in making interactive 3D interface, is able to meet my need perfectly. After playing with this for almost a week, I find my way.

#=============================================================================
# Scripts updated from:
 
#----------------------------------------------------------------------------
# Load packages
#----------------------------------------------------------------------------
library(rgl)
library(gplots)
require(OpenStreetMap)
require(ggplot2)
 
#----------------------------------------------------------------------------
# Zoom in TX and LA
Houston_region <- c(-96.5,-94,29,31.5) # lon_left,lon_right,lat_down,lat_up
Res            <- 0.02                 # Resolution
Row_down       <- (Houston_region[3]-25 )/Res + 1
Row_up         <- (Houston_region[4]-25 )/Res
Col_left       <- (Houston_region[1]+110)/Res + 1
Col_right      <- (Houston_region[2]+110)/Res
 
#----------------------------------------------------------------------------
# Load HCHO data
#----------------------------------------------------------------------------
 
# OMI 2005-2008
filename       <- "data/Oversampling/OMI_oversampling_South_US_2005_2008_MJJAS_new.dat"
 
Overpass_data  <- read.table(file=filename,header=F)
Nlines         <- dim(Overpass_data)[1] # Total lines
head(Overpass_data)
 
Lat_grid       <- seq( min(Overpass_data$V3), max(Overpass_data$V3), by=0.02)
Lon_grid       <- seq( min(Overpass_data$V4), max(Overpass_data$V4), by=0.02)
 
NRows <- max(Overpass_data$V1)
NCols <- max(Overpass_data$V2)
#
VCD_fine_value_OMI <- array(0, dim=c(NRows,NCols))
 
for (line in 1:Nlines){
VCD_fine_value_OMI[Overpass_data$V1[line],Overpass_data$V2[line]] <- Overpass_data$V5[line]
}
 
HCHO <- VCD_fine_value_OMI[Row_down: Row_up, Col_left: Col_right]/5e+15 # Exaggerate the relief
HCHO <- t(HCHO)
 
#----------------------------------------------------------------------------
# Make a base map
#----------------------------------------------------------------------------
 
map3d <- function(map,...){
  if(length(map$tiles)!=1){
    stop("multiple tiles not implemented")
  }
 
  nx = map$tiles$xres
  ny = map$tiles$yres
 
  xmin = map$tiles$bbox$p1[1]
  xmax = map$tiles$bbox$p2[1]
  ymin = map$tiles$bbox$p1[2]
  ymax = map$tiles$bbox$p2[2]
 
  xc = seq(xmin,xmax,len=ny)
  yc = seq(ymin,ymax,len=nx)
  colours = matrix(map$tiles$colorData,ny,nx)
  m = matrix(min(HCHO),ny,nx)
  surface3d(xc,yc,m,col=colours)
}
 
#----------------------------------------------------------------------------
# Add 3D surface map
#----------------------------------------------------------------------------
 
map <- openmap(c(Houston_region[4], Houston_region[1]),c(Houston_region[3], Houston_region[2]),8,'stamen-toner')
map <- openproj(map)
map3d(map)
rgl.light(90,90) # Add additional light
 
Lat      <-  29+0.5*Res+Res*((1:nrow(HCHO))-1)
Lon      <- -96.5+0.5*Res+Res*((1:ncol(HCHO))-1)
HCHOlim  <- range(HCHO)
colorlut <- palette(rev(rich.colors(32)))
col      <-32-31*((HCHO-HCHOlim[1])/(HCHOlim[2]-HCHOlim[1]))
persp3d(Lon,Lat,HCHO,alpha=0.1,color=col,add=T)
axes3d(edges = "bbox", labels = TRUE, tick = TRUE, nticks = 5, box=FALSE, expand = 1.03,col="black",lwd=3)
 
# In case you want a rational view
#rgl.viewpoint(theta = 0, phi =-90)
#for(i in seq(0, 90, by = 0.15)) {
#  rgl.viewpoint(theta = 0, phi =-90+i)
#}

#=============================================================================

Here is what I got: