How to install wgrib2 in OSX

Prompted by both my own struggles with wgrib2 compilation and a plea on the rNOMADS email listserv, I’m going to describe how to compile and install wgrib2 on Mac OS.

First of all, some background: wgrib2 is an excellent utility written by Wesley Ebisuzaki at NOAA.  It allows for a number of swift and stable operations on GRIB2 files (a common file format for weather and climate data).  It is also a requirement for using grib files in rNOMADS (the function ReadGrib() in particular).

So here’s how to install it on Mac OS.

  1.  Get Command Line Tools for Xcode (search for it on Duckduckgo or your engine of choice) and also make sure gcc is installed. If you don’t know what gcc is, stop now and find someone who does (it will save you a lot of time).
  2. Download wgrib2 here (note download links are pretty far down the page).
  3. Untar the tarball somewhere, and roll up your sleeves.  cd into the resulting wgrib directory.
  4.  In the makefile, uncomment the lines

    #export cc=gcc
    #export FC=gfortran
    Also search for makefile.darwin in the makefile and uncomment the line containing it.  You’ll see instructions to this effect in the makefile anyway.
  5. Now we have to edit the included libpng package, since it is untarred by the makefile and doesn’t inherit our compiler specifications in step 4. Ensuring that we’re in the wgrib directory:

    tar -xvf libpng-1.2.57.tar.gz
    cd libpng-1.2.57/scripts/

    now edit the makefile.darwin file, changing

    CC=cc

    to

    CC=gcc

    Now, return to the wgrib directory, and re-tar libpng!

    tar -cf libpng-1.2.57.tar libpng-1.2.57
    gzip libpng-1.2.57.tar

    If it asks you if you want to replace the original tar.gz file, say “yes”.  What we’ve done here is edited libpng to make sure it uses the right compiler.
  6. Finally, type make to build wgrib2.

If you are still having problems (for example, libaec complaining that there is no c compiler), make sure that all the compiler commands (gcc, cc, etc.) all point to something other than clang (the default compiler that comes with OSX).  You may have to edit your bash_profile file to ensure this.

As always, contact me on the form below if you’re having unresolvable issues.

Advertisements

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 package has been published in the journal Computers & Geosciences

I’m happy to announce that my article “Near Real Time Weather and Ocean Model Data Access with rNOMADS” was recently accepted for publication in the journal Computers & Geosciences. The editor has kindly provided me with a link to the article that’s free until late April, check it out here.  I’m excited to have the package published somewhere official and I’m looking forward to seeing who ends up using it.

In other news, I’m still waiting for the tropospheric jet to point anywhere besides East so we can launch a tracked, camera bearing solar balloon from Chapel Hill.   I wrote an rNOMADS script that generates wind profiles from the 0600 GMT model run on today’s date out to three days from now, and put it on my computer’s crontab.  Every day, I look at the winds, and wait.  Here’s the profiles for today, let’s see if it updates daily on wordpress too:

Wind profiles above Chapel Hill, NC, USA

You can have a look at the script that generates this image here.

Note that the “montage” system command uses the Ubuntu imagemagick package to combine separate PNGs.  If you’re on windows or mac, this will not work for you.

New updates to the rNOMADS package and big changes in the GFS model

I rolled out a big update to the rNOMADS package in R about two weeks ago.  Now, the list of real time weather, ocean, and sea ice models available through rNOMADS updates automatically by scraping the NOMADS web site.  This way, changes in model inventories will be instantly reflected in rNOMADS without the need for a new version release.

Keep abreast of future updates to rNOMADS by subscribing to the mailing list here.  Feel free to ask for help or make comments on this list as well.

In other news, NOAA just updated the Global Forecast System to provide 0.25 x 0.25 degree output – doubling the resolution of the model!  Check out this crystal clear views of surface temperatures across the planet (source code below the image):

World temperature at 2 m above ground using the 0.25 x 0.25 degree output of the Global Forecast System model.

World temperature at 2 m above ground using the 0.25 x 0.25 degree output of the Global Forecast System model.

 


library(rNOMADS)
library(GEOmap)

#Get dates of model output
model.urls <- GetDODSDates(“gfs_0p25”)

#Find day of most recent model run
latest.model <- tail(model.urls$url, 1)

#Find most recent model run on that day
model.runs <- GetDODSModelRuns(latest.model)

#Get the most recent model (excluding analysis only)
latest.model.run <- tail(model.runs$model.run[which(grepl(“z$”, model.runs$model.run))], 1)

#Define model domain
time <- c(0,0) #Analysis model
lon <- c(0, 1439) #All longitude points
lat <- c(0, 720) #All latitude points
variables <- c(“tmp2m”) #Temperature 2 m above ground

#Get data from NOMADS real time server
tmp.data <- DODSGrab(latest.model, latest.model.run,
variables, time, lon, lat, display.url = FALSE)

#Reformat it
tmp.grid <- ModelGrid(tmp.data, c(0.25, 0.25))

#Define color scale
colormap <- rev(rainbow(500, start = 0 , end = 5/6))

#Plot it
image(x = tmp.grid$x, y = sort(tmp.grid$y), z = tmp.grid$z[1,1,,], col = colormap,
xlab = “Longitude”, ylab = “Latitude”,
main = paste(“World Temperature at Ground Level:”,
tmp.grid$fcst.date))

plotGEOmap(coastmap, border = “black”, add = TRUE,
MAPcol = NA)

rNOMADS 2.0.2 released

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.

Wind profiles at each infrasound station in the Transportable Array (currently on the US East Coast). Ground surface is at the center of the plot, the top of the stratosphere is at the outer radius. Wind magnitudes are denoted by color, and azimuth by position around the circle.

Wind profiles at each infrasound station in the Transportable Array (currently on the US East Coast). Ground surface is at the center of the plot, the top of the stratosphere is at the outer radius. Wind magnitudes are denoted by color, and azimuth by position around the circle.


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_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, 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)")


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: