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

How to install rNOMADS with GRIB file support on Windows

Two years ago, I wrote a software package for R called “rNOMADS” that interfaces with online weather and sea ice model repositories to gather data in real time, for free.  The data are delivered in two ways: a simple, pure R, cross platform interface using GrADS-DODS, and binary files in GRIB format.  The one issue with GRIB is that this format can’t be read directly into R; it requires the external program “wgrib2.”  Installing rNOMADS with GRIB support for Linux is covered in this post (Mac OS is probably similar).  I thought GRIB support for Windows was impossible until the guy who runs this blog casually told me he’d figured it out.  So, I finally got around to trying it myself, and I’m happy to say I got it to work!  Here’s how:

Step 1:

Download the most recent wgrib2.exe version and all DLL files from Wesley Ebisuzaki’s web site:  Click here.

Step 2:

Make a directory somewhere on your computer.   I chose:  C:\Program Files\wgrib2

Copy the .exe and the .dll files to this directory.

Step 3:

Append the directory path to the Windows PATH variable.  Find out how to do this here.

Edit: a user points out that sometimes this doesn’t work.  If not, change the path variable in R using Sys.setenv():

Sys.setenv("PATH" = "C:/Program Files/wgrib2")

 

Step 4:  If R is open, close and reopen it.  Then try the following command in the R interpreter:
system("wgrib2")

If you get a bunch of text that looks like the image below, you’ve succeeded.  You can start reading GRIB files with rNOMADS!

Output of wgrib2 call from R.

Output of wgrib2 call from R.

rNOMADS 2.0.1 released

The latest version of rNOMADS is now available on CRAN.  This update resolves several minor bugs and one major one involving multiple variable/level selections when using the ModelGrid function.  I have also added support for two more models on NOMADS:  Climate Forecast System Flux Products and Climate Forecast System 3D Pressure Products.  This brings the total number of real time model products supported by rNOMADS to 57.  A quick test of the Climate Forecast System Flux Products produced the following image of surface albedo for 1200 GMT on July 4th, 2014:

Albedo of the Earth's surface at 1200 GMT, July 4th, 2014.

Albedo of the Earth’s surface at 1200 GMT, July 4th, 2014.

In other news, data curators at the University Corporation for Atmospheric Research (UCAR) are currently investigating rNOMADS as an open source software package to read their GRIB archive.  UCAR hosts a series of model archive data, some of which is freely available to the public.  Their specific interest is in using rNOMADS to read their archive of the FNL model.  They have gotten wgrib2 to work on MAC OS X, allowing rNOMADS GRIB functionality to be extended to Macintosh machines.