I uploaded the newest version of rNOMADS to CRAN yesterday. This one has a new plotting function for wind altitude azimuth, and magnitude (see below for plot and source code). I also added a function for reading GRIB inventories, fixed a few typos, and tweaked a few functions. In other news, the rNOMADS subversion repository has migrated to R-Forge: http://rnomads.r-forge.r-project.org/.

I’ve also set up an rNOMADS mailing list; subscribe here if you need to ask for help or want to hear about new releases.

library(rNOMADS) download.file("http://www.unc.edu/~haksaeng/rNOMADS/myTA.RDATA", destfile = "myTA.RDATA") load("myTA.RDATA") #Find the latest Global Forecast System model run model.urls <- GetDODSDates("gfs_hd") 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) & lons >= min(myTA$lon + 360)) - 1 lat.ind <- which(lats <= max(myTA$lat) & lats >= min(myTA$lat)) - 1 levels <- c(0, 46) time <- c(0, 0) #Get data hgt.data <- DODSGrab(latest.model, latest.model.run, "hgtprs", time, c(min(lon.ind), max(lon.ind)), c(min(lat.ind), max(lat.ind)), levels) ugrd.data <- DODSGrab(latest.model, latest.model.run, "ugrdprs", time, c(min(lon.ind), max(lon.ind)), c(min(lat.ind), max(lat.ind)), levels) vgrd.data <- DODSGrab(latest.model, latest.model.run, "vgrdprs", time, c(min(lon.ind), max(lon.ind)), c(min(lat.ind), max(lat.ind)), levels) #Reformat the data hgt.grid <- ModelGrid(hgt.data, c(0.5, 0.5)) ugrd.grid <- ModelGrid(ugrd.data, c(0.5, 0.5)) vgrd.grid <- ModelGrid(vgrd.data, c(0.5, 0.5)) #Build profiles zonal.wind <- c() meridional.wind <- c() height >- c() for(k in 1:length(myTA$lon)) { hgt.profile <- BuildProfile(hgt.grid, myTA$lon[k], myTA$lat[k], spatial.average = TRUE) ugrd.profile <- BuildProfile(ugrd.grid, myTA$lon[k], myTA$lat[k], spatial.average = TRUE) vgrd.profile <- BuildProfile(vgrd.grid, myTA$lon[k], myTA$lat[k], spatial.average = TRUE) synth.hgt <- seq(min(hgt.profile), max(hgt.profile), length.out = 1000) ugrd.spline <- splinefun(hgt.profile, ugrd.profile, method = "natural") vgrd.spline <- splinefun(hgt.profile, vgrd.profile, method = "natural") zonal.wind[[k]] <- ugrd.spline(synth.hgt) meridional.wind[[k]] <- vgrd.spline(synth.hgt) height[[k]] <- synth.hgt print(paste("Finished", k, "of", length(myTA$lon), "profiles!")) } #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)", elev.labels.az = 15)