Wind velocity from flight

## ======================================================================================
## Velocity vectors 
## implemented by Karline Soetaert
## ======================================================================================

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

require(RNetCDF)
require(OceanView)  

# Aircraft track files used by Zebra and IDV from BAMEX field project. 
# A flight is broken down into separate hourly files (based on clock time not flight time).
Nc <- open.nc("WMI_Lear.nc")
print.nc(Nc)
WMI <- data.frame(time = var.get.nc(Nc, "time"),
                  altitude = var.get.nc(Nc, "altitude"),    #km
                  latitude = var.get.nc(Nc, "latitude"),    #dgN
                  longitude = var.get.nc(Nc, "longitude"),  #dgE
                  pressure  = var.get.nc(Nc, "pressure"),   #Hpa
                  tdry  = var.get.nc(Nc, "tdry"),       #temperature
                  dp  = var.get.nc(Nc, "dp"),           #dewpoint temperature
                  mr  = var.get.nc(Nc, "mr"),           #mixing ratio
                  wspd  = var.get.nc(Nc, "wspd"),       #wind speed
                  wdir  = var.get.nc(Nc, "wdir"))       #wind direction
Wind.x <- WMI$wspd * cos(WMI$wdir/360*2*pi)
Wind.y <- WMI$wspd * sin(WMI$wdir/360*2*pi)

par(mfrow = c(1, 1), oma = c(2, 0, 0, 0))
vectorplot(Wind.x, Wind.y, x = WMI$time - WMI$time[1], colvar = WMI$altitude,
  xlab = "time", ylab = "m/s", #clab = "height, km", 
  main = "wind velocity from flight", colkey = FALSE)
colkey (side = 3, length = 0.5, width = 0.5, dist = -0.71, shift = -0.2, 
  clim = range(WMI$altitude), add = TRUE, cex.axis = 0.8, 
  mgp = c(1, 0.5, 0), clab = "height, km")

par(new = TRUE)
par(fig = c(0.5, 1.0, 0.45, 1.0))

range(WMI$latitude)
range(WMI$longitude)

require(RgoogleMaps)
destfile = "WMI.png"
MetaInfo <- GetMap(center = c( lat = mean(range(WMI$latitude)),  #size = c(300, 300), 
  lon = mean(range(WMI$longitude))), zoom = 9, maptype = "roadmap", destfile = destfile)

require(png)
img <- as.raster(readPNG(destfile))

x <- seq(MetaInfo$BBOX$ll[2], MetaInfo$BBOX$ur[2], length.out = MetaInfo$size[2])
y <- seq(MetaInfo$BBOX$ll[1], MetaInfo$BBOX$ur[1], length.out = MetaInfo$size[1])
xlim <- range(x)
ylim <- range(y)
image2D(x, y, z = img[nrow(img):1,], 
       xlab = expression(""^o~E), ylab = expression(""^o~N), 
       xaxs = "i", yaxs = "i", cex.lab = 0.6, cex.axis = 0.6)
with (WMI, scatter2D(y = latitude, x = longitude, colvar = altitude, add = TRUE,
      type = "l", lwd = 3, colkey = FALSE))