NCEP/NCAR relative humidity

## =============================================================================
## Plotting 3D data - assume files have been downloaded
## implemented by Karline Soetaert
## =============================================================================

# example data files from 
url <- "http://www.unidata.ucar.edu/software/netcdf/examples/rhum.2003.nc"

require(RNetCDF)
require(OceanView)

#NCAR S-band dual-polarized (S-Pol) radar netCDF sweep file format. Range Height Indicator (RHI) files are vertical scans, Plan Position Indicator (PPI) files are horizontal scans. Each file contains one sweep (azimuth (up/down) for RHI or elevation (horizontal) for PPI).

#NCEP/NCAR Reanalysis 1 data for relative humidity at multiple levels, daily averages for 2003 
Nc <- open.nc("rhum.2003.nc")
lat <- var.get.nc(Nc, "lat") #dgN
lon <- var.get.nc(Nc, "lon") #dgE
time <- var.get.nc(Nc, "time") #hours since 1/1/1
level<- var.get.nc(Nc, "level") #millibar
rhum <- var.get.nc(Nc, "rhum") * 0.01 + 302.66 #lon, lat, level, time

ii <- c(1, 4, 7)

# Show data from all levels
image2D(x = lon, y = lat, z = rhum[,, ,1], shade = 0.1, 
  xlab = "dgE", ylab = "dgN", subset = ii, 
  main = paste ("level ",level[ii]), 
  clim = c(0, 100), colkey = FALSE)
#colkey(clim = c(0, 100), clab = c("humidity", "%"))

# Volumetric data
M <- mesh(lon, lat, level)

image3D(x = lon, y = lat, z = 1000, colvar = rhum[,,1,1], 
  alpha = 0.5, plot = FALSE, clim = c(0, 100), zlim = range(level),
  main = "NCEP/NCAR relative humidity",
  clab = "%", xlab = "dgE", ylab = "dgN", zlab = "level, mB",
  ticktype = "detailed")

# Isosurface 
isosurf3D(x = lon, y = lat, z = level, 
  colvar = rhum[,,,1], level = 20, 
  col = jet.col()[20], alpha = 0.3, plot = FALSE, add = TRUE)  
image3D(x = lon, y = lat, z = 300, colvar = rhum[,,8,1], 
  alpha = 0.5, add = TRUE, plot = FALSE, clim = c(0, 100), 
  colkey = FALSE)
plotdev()

# try:
#plotrgl()