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: