Accessing ensemble weather models with rNOMADS

Downloading a weather model and using it to predict stuff is all well and good – but what are the sources of error?  Specifically, how accurate is the atmospheric parametrization?   Luckily, NOMADS provides sets of models with slightly different initialization values.  The result is a set of realizations that captures the variability of the model product.  Here’s how to get that product, and examples of what you can do with it.


library(rNOMADS)

#Get the latest ensemble model run
model.urls <- GetDODSDates("gens")
latest.model <- tail(model.urls$url, 1)
model.runs <- GetDODSModelRuns(latest.model)
model.run <- tail(model.runs$model.run[grepl("all",
    model.runs$model.run)], 1)

#Define region of interest: Chapel Hill, NC
lon <- -79.052104
lat <- 35.907553

#Set up DODS model grid indices
#Use GetDODSModelRunInfo() to figure out how to do this for other products
lons <- seq(0, 359, by = 1)
lats <- seq(-90, 90, by = 1)

lon.diff <- abs(lon + 360 - lons)
lat.diff <- abs(lat - lats)

model.lon.ind <- which(lon.diff == min(lon.diff)) - 1
model.lat.ind <- which(lat.diff == min(lat.diff)) - 1

#Set up call to NOMADS GrADS-DODS server
#Analysis(?) model only
time <- c(0, 0)
#Longitude grid
node.lon  <- c(model.lon.ind - 2, model.lon.ind + 2)
#Latitude grid
node.lat  <- c(model.lat.ind - 2, model.lat.ind + 2)
#Temperature, wind speeds, and geopotential height
variables <- c("tmpprs", "ugrdprs", "vgrdprs", "hgtprs")
levels    <- c(0, 25) #All available levels
ensemble <- c(0, 20)  #All available ensembles

model.data <- DODSGrab(latest.model, model.run, variables,
    time, node.lon, node.lat, levels = levels, ensemble = ensemble)

#Plot temperature variability
plot(c(-65, 25), c(0, 35), xlab = "Temperature (C)",
    ylab = "Geopotential Height (km)", type = 'n')

for(k in ((ensemble[1]:ensemble[2] + 1))) {
    model.data.sub <- SubsetNOMADS(model.data, ensemble = c(k),
        variables = c("hgtprs", "tmpprs"))
    profile <- BuildProfile(model.data.sub, lon, lat)
    points(profile[[1]]$profile.data[,
        which(profile[[1]]$variables == "tmpprs"), 1] - 272.15,
        profile[[1]]$profile.data[,
        which(profile[[1]]$variables == "hgtprs"), 1] / 1000)
}

#Plot winds
zonal.wind <- NULL
merid.wind <- NULL
height     <- NULL

for(k in ((ensemble[1]:ensemble[2] + 1))) {
    model.data.sub <- SubsetNOMADS(model.data,
        ensemble = c(k), variables = c("hgtprs",
        "ugrdprs", "vgrdprs"))
    profile <- BuildProfile(model.data.sub, lon, lat)
    hgt     <- profile[[1]]$profile.data[,
        which(profile[[1]]$variables == "hgtprs"),]
    ugrd    <- profile[[1]]$profile.data[,
        which(profile[[1]]$variables == "ugrdprs"),]
    vgrd    <- profile[[1]]$profile.data[,
        which(profile[[1]]$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)
   merid.wind[[k]] <- vgrd.spline(synth.hgt)
   height[[k]] <- synth.hgt
}

dev.new()
PlotWindProfile(zonal.wind, merid.wind, height, lines = TRUE,
    points = FALSE, elev.circles = c(0, 15000, 30000),
    elev.labels = c(0, 15, 30), radial.lines = seq(45, 360, by = 45),
    colorbar = TRUE, invert = FALSE,
    point.cex = 2, pch = 19, lty = 1, lwd = 1,
    height.range = c(0, 30000), colorbar.label = "Wind Speed (m/s)")

Here’s atmospheric profiles from GFS ensemble runs on November 24, 2016 for Chapel Hill, NC. Temperature is quite variable at low altitudes (several degrees Celsius) and at 5 km above sea level (a couple degrees). Ensemble runs are more consistent for the rest of the atmosphere.

chapel_hill_temperature

Temperature profile above the Geology Department at UNC Chapel Hill for November 24, 2016 at 18:00 UTC.

Wind speeds are pretty consistent, at least at the resolution of the color scale I’m using.  Wind direction varies by a few degrees below 15 km, but becomes quite a bit less reliable between 15 and 30 km altitude.  This is one reason that our solar balloon trajectories are so hard to predict.  The other reason is because we still don’t really understand the physics of solar balloon flight, but that’s a subject for another post.

chapel_hill_wind_profile

Extracting North American Mesoscale (NAM) Model Data with rNOMADS using DODS

The (relatively) recent upgrade of rNOMADS to version 2.0.0 allows cross platform support for downloading atmospheric and oceanic operational model data into R via the GrADS-DODS system.  However, users have been experiencing difficulty using DODS to get North American Mesoscale (NAM) real time data.  The script below shows an example of getting the NAM 12z model for North America using DODS.  Below the script and the figure, I’ll talk about why getting NAM data via DODS is different than, say, GFS.

#Get data for NAM subsets using DODS
library(rNOMADS)
library(rgl)
#Individual NAM models are listed separately in GRIB but all together in DODS.
#So we get the dates of the NAM runs first, then we have to sort out which "type" of NAM we want.

all.nam.dates <- GetDODSDates("nam", request.sleep = 1)

#Here, we pull the model runs.  This is when we can get individual models out.
#Get the most recent one:
all.nam.mr <- GetDODSModelRuns(tail(all.nam.dates$url, 1))

#For example, let's get North America "nam_na" model, 12z model run.
namna.model.runs <- all.nam.mr$model.run[grepl("nam_na_12z", all.nam.mr$model.run)]

#We can get info on what that model contains by inspecting the output from this:
info <- GetDODSModelRunInfo(tail(all.nam.dates$url, 1), tail(namna.model.runs, 1))

#Get data on surface solar flux across North America
sf.data <- DODSGrab(tail(all.nam.dates$url, 1), tail(namna.model.runs, 1), "csdsfsfc", 0, c(0,706), c(0, 285))

#Now we can plot it with rgl 
#First, set the "missing data" values to -1 (or whatever)

#Credit to Dr. David Forrest for this display code:
valScrubbed<-as.numeric(sf.data$value)
valScrubbed[valScrubbed==9.999e+20]<- -1
with(sf.data,plot3d(lon,lat,valScrubbed,col=heat.colors(100)[cut(valScrubbed,breaks=100)],asp=c(1/cos(35*pi/180),1,.2)))

Solar radiation flux over North America from the NAM NA model.

Solar radiation flux over North America from the NAM NA model.

The reason we can’t just get a specific NAM subset model using (for example)

GetDODSDates("nam_na", request.sleep = 1)

is because the DODS server is set up differently than the GRIB server for NAM. It turns out that all the NAM models are stored together in one directory (in GRIB, each one has its own directory), and so we actually have to sort out individual models at the model run level, not the date level. So the way the above script deals with this is:
1. Figure out all dates when NAM was run.
2. Figure out the model runs for a specific date.
3. Pull out the specific NAM model subset you want (e.g. nam_na, nam_ak, whatever).
4. Get the model data.

If you have questions, don’t hesitate to contact me: