if (FALSE) {
download.file("http://www.gmail.com/~haksaeng/rNOMADS/myTA.RDATA",
destfile = "myTA.RDATA")
load("myTA.RDATA")
#Find the latest Global Forecast System model run
model.urls <- GetDODSDates("gfs_0p50")
latest.model <- tail(model.urls$url, 1)
model.runs <- GetDODSModelRuns(latest.model)
latest.model.run <- tail(model.runs$model.run, 1)
#Get model nodes
lons <- seq(0, 359.5, by = 0.5)
lats <- seq(-90, 90, by = 0.5)
lon.ind <- which(lons <= (max(myTA$lon + 360) + 1) & lons >= (min(myTA$lon + 360) - 1))
lat.ind <- which(lats <= (max(myTA$lat) + 1) & lats >= (min(myTA$lat) - 1))
levels <- c(0, 46)
time <- c(0, 0)
#Get data
variables <- c("hgtprs", "ugrdprs", "vgrdprs")
model.data <- DODSGrab(latest.model, latest.model.run,
variables, time, c(min(lon.ind), max(lon.ind)),
c(min(lat.ind), max(lat.ind)), levels)
#Build profiles
profile <- BuildProfile(model.data, myTA$lon + 360, myTA$lat,
spatial.average = FALSE)
#Build profiles
zonal.wind <- NULL
meridional.wind <- NULL
height <- NULL
for(k in 1:length(profile)) {
hgt <- profile[[k]]$profile.data[, which(profile[[k]]$variables == "hgtprs"),]
ugrd <- profile[[k]]$profile.data[, which(profile[[k]]$variables == "ugrdprs"),]
vgrd <- profile[[k]]$profile.data[, which(profile[[k]]$variables == "vgrdprs"),]
synth.hgt <- seq(min(hgt),
max(hgt), length.out = 1000)
ugrd.spline <- splinefun(hgt, ugrd, method = "natural")
vgrd.spline <- splinefun(hgt, vgrd, method = "natural")
zonal.wind[[k]] <- ugrd.spline(synth.hgt)
meridional.wind[[k]] <- vgrd.spline(synth.hgt)
height[[k]] <- synth.hgt
}
#Plot them all
PlotWindProfile(zonal.wind, meridional.wind, height, lines = TRUE,
points = FALSE, elev.circles = c(0, 25000, 50000), elev.labels = c(0, 25, 50),
radial.lines = seq(45, 360, by = 45), colorbar = TRUE, invert = FALSE,
point.cex = 2, pch = 19, lty = 1, lwd = 1,
height.range = c(0, 50000), colorbar.label = "Wind Speed (m/s)")
}
Run the code above in your browser using DataLab