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,707), c(0, 286))

#Now we can plot it with rgl (do NOT use ModelGrid, it has the wrong projection)
#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:

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.

Downloading weather, sea ice, and wave model data with the rNOMADS package in R

The NOAA Operational Model Archive and Distribution System is a treasure trove of near real time and archived model outputs describing global and regional weather, sea ice, and wave data.  I developed the rNOMADS package about a year ago to make this available to R users. In this post, I’ll present some source code and a couple of figures showing a few of the useful things you can do with rNOMADS.

For detailed examples showing rNOMADS with GRIB file support (Linux only) see the vignette here, for a cross platform version of the same, see here.

1.  Getting wind speed at a specific point

library(rNOMADS)

#A location near my house
lat <- 35.828304
lon <- -79.107467

#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 nearest model nodes
lons <- seq(0, 359.5, by = 0.5)
lats <- seq(-90, 90, by = 0.5)
lon.diff <- abs(lon + 360 - lons)
lat.diff <- abs(lat - lats)
model.lon.ind <- which(lon.diff == min(lon.diff)) - 1 #NOMADS indexes at 0
model.lat.ind <- which(lat.diff == min(lat.diff)) - 1 

#Subset model
time <- c(0,0) #Model status at initialization
lon.inds <- c(model.lon.ind - 2, model.lon.ind + 2)
lat.inds <- c(model.lat.ind - 2, model.lat.ind + 2)
variable <- "ugrd10m" #First get east-west wind

u.grd.data <- DODSGrab(latest.model, latest.model.run, variable,
   time, lon.inds, lat.inds)

variable <- "vgrd10m" #Then get north-south wind

v.grd.data <- DODSGrab(latest.model, latest.model.run, variable,
   time, lon.inds, lat.inds)

#Reformat the data
u.grd.grid <- ModelGrid(u.grd.data, c(0.5, 0.5))
v.grd.grid <- ModelGrid(v.grd.data, c(0.5, 0.5))

#Interpolate it to the point of interest
u.point <- BuildProfile(u.grd.grid, lon, lat, TRUE)
v.point <- BuildProfile(v.grd.grid, lon, lat, TRUE)

#Present results!
print(paste("The East-West winds at Briar Chapel are going", sprintf("%.2f", u.point),
    "meters per second, and the north-south winds are going", sprintf("%.2f", v.point),
    "meters per second."))

#How did I know all these strange parameter names?
info <- GetDODSModelRunInfo(latest.model, latest.model.run)
print(info)

2. Getting a temperature profile from 0 to 40 km above a specific point

library(rNOMADS)

#A location near my house
lat <- 35.828304
lon <- -79.107467

#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 nearest model nodes
lons <- seq(0, 359.5, by = 0.5)
lats <- seq(-90, 90, by = 0.5)
lon.diff <- abs(lon + 360 - lons)
lat.diff <- abs(lat - lats)
model.lon.ind <- which(lon.diff == min(lon.diff)) - 1 #NOMADS indexes at 0
model.lat.ind <- which(lat.diff == min(lat.diff)) - 1 

#Subset model
time <- c(0,0) #Model status at initialization
lon.inds <- c(model.lon.ind - 2, model.lon.ind + 2)
lat.inds <- c(model.lat.ind - 2, model.lat.ind + 2)
levels <- c(0, 46) #All pressure levels
variable <- "tmpprs" #First get temperature

tmpprs.data <- DODSGrab(latest.model, latest.model.run, variable,
   time, lon.inds, lat.inds, levels)

variable <- "hgtprs" #Now get elevation of each pressure level
hgtprs.data <- DODSGrab(latest.model, latest.model.run, variable,
   time, lon.inds, lat.inds, levels)

#Reformat the data
tmpprs.grid <- ModelGrid(tmpprs.data, c(0.5, 0.5))
hgtprs.grid <- ModelGrid(hgtprs.data, c(0.5, 0.5))

#Interpolate it to the point of interest
tmpprs.point <- BuildProfile(tmpprs.grid, lon, lat, TRUE)
hgtprs.point <- BuildProfile(hgtprs.grid, lon, lat, TRUE)

#Plot it!
plot(tmpprs.point - 272.15, hgtprs.point, xlab = "Temperature (C)",
   ylab = "Elevation (m)", main = paste("Temperature above Chapel Hill, NC, at",
   tmpprs.grid$model.run.date))
Temperature profile of troposphere and stratosphere above Chapel Hill, NC, on May 28.

Temperature profile of troposphere and stratosphere above Chapel Hill, NC, on May 28.

2. A world map of surface temperature

library(GEOmap)
library(rNOMADS)

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)

time <- c(0,0) #Analysis model
lon <- c(0, 719) #All 720 longitude points
lat <- c(0, 360) #All 361 latitude points

tmp2m.data <- DODSGrab(latest.model, latest.model.run,
    "tmp2m", time, lon, lat)
atmos <- ModelGrid(tmp2m.data, c(0.5, 0.5), "latlon")
colormap <- rev(rainbow(500, start = 0 , end = 5/6))
image(atmos$x, sort(atmos$y), atmos$z[1,1,,], col = colormap,
    xlab = "Longitude", ylab = "Latitude",
    main = paste("World Temperature at Ground Level:", atmos$model.run.date))
plotGEOmap(coastmap, border = "black", add = TRUE,
    MAPcol = NA)
World temperature on May 28.

World temperature on May 28.

2. Wave heights in the northwest Atlantic ocean


library(rNOMADS)
library(GEOmap)
library(RCurl)

model.urls <- GetDODSDates("wave")
latest.model <- tail(model.urls$url, 1)
model.runs <- GetDODSModelRuns(latest.model) #West Atlantic waves
latest.model.run <- tail(model.runs$model.run, 1)
time <- c(0,0)
lon <- c(0, 274)
lat <- c(0, 202)

wave.data <- DODSGrab(latest.model, latest.model.run,
    "htsgwsfc", time, lon, lat)
wave.grid <- ModelGrid(wave.data, c(0.25, 0.25))
#Remove "no data" values
wave.grid$z[which(wave.grid$z>1e10, arr.ind=TRUE)] <- NA
colormap <- rainbow(500, start=0, end=5/6)
image(wave.grid$x, sort(wave.grid$y), wave.grid$z[1,1,,], col = colormap,
xlab = "Longitude", ylab = "Latitude",
main = paste("Wave Height:", wave.grid$model.run.date))
plotGEOmap(coastmap, border = "black", add = TRUE,
MAPcol = "black")

Wave heights in the north Atlantic on July 4th, 2014. High seas from Hurricane Arthur are visible just off the US East Coast.

Anyone seen a red lunch box with cameras attached?

Catastrophe!  Yesterday, while launching a hydrogen balloon in Chapel Hill, North Carolina, our satellite tracker fell out of the payload.  The camera box, however, cheerfully continued its ascent into the wild blue yonder.  The result?  There is a red lunch box with three cameras in it somewhere in North Carolina.

Detailed description of payload as it will appear on the ground:

Shreds of mylar and a balloon nozzle, followed by a 40′ string, then a plastic parachute (may look like a plastic bag from far away).  Another length of string approximately 10′ long, followed by a red lunch box with “Solar Balloon Payload” and my phone number written on it.  The lunch box has three cameras, one pointing down, one pointing out the side, and a smaller one pointing up.  A milk jug with a sticker saying “Soleil Multimedia” may also still be attached.

The projected landing zone was in north/northwest Harnett County.  I believe that is reasonably accurate.  So our payload, with spectacular footage, is somewhere out there.  I want it back.

Here’s a couple of potential flight trajectories:

Trajectory 1

Trajectory 2

Trajectory 3

If you have any information or know anyone in the region who might be able to help, please contact us.

Accessing global weather model data using the rNOMADS package in R

*** UPDATE ***

rNOMADS 2.0.0 was released yesterday (5-15-2014), and it now has cross-platform support (meaning you can use its data retrieval capabilities in Windows and Mac as well as Linux).

Stay tuned for more tutorials and usage examples in the coming days.

***

The rNOMADS package interfaces with the NOAA Operational Model Archive and Distribution System to provide access to 55 operational (i.e. real time and prediction) models describing the state of the ocean and the atmosphere. rNOMADS has been used to get wind and wave data for a real time sailing game, to quantify solar energy available for power plants in Europe, and to predict helium balloon flights. We look forward to continuing to hear about novel and unexpected uses for this spectacular data set.

In this post, we show how to use rNOMADS to do the following:
1. Plot global temperature and wind at different atmospheric pressure levels
2. Examine the distribution of wave heights in the Atlantic ocean
3. Produce a precise atmospheric temperature and wind profile for a specific point at a specific time
4. Generate simultaneous atmospheric profiles at over 100 locations on the Eastern Seaboard of the USA

Links to source code are provided below each set of figures. A link to rNOMADS installation instructions is provided at the end of the post.

Global Temperature and Wind Maps

The Global Forecast System (GFS) model provides weather data on a 0.5 x 0.5 degree grid for the entire planet.  It is run 4 times daily, and produces a prediction every three hours out to 180 hours.  The zero hour “analysis” forecast is the state of the atmosphere at the model run time, and it uses some observational data to increase its accuracy.  Here, we plot the temperature at 2 m above the surface, the wind at 10 m above the surface, and the winds at 300 mb (around 9 kilometers elevation).  The 300 mb plot often shows the northern and southern jet streams quite well.

Temperature at the Earth's surface determined using the Global Forecast System model.

Temperature at the Earth’s surface determined using the Global Forecast System model.

Winds at the surface of the Earth from the GFS model.  Note the little spot of high winds - that's Tropical Cyclone Gillian, a Category 3 storm when this image was generated.

Winds at the surface of the Earth from the GFS model. Note the little spot of high winds south of Indonesia – that’s Tropical Cyclone Gillian, a Category 3 storm when this image was generated.

Jet streams and Rossby waves are clearly visible in this image of the wind speeds of the upper troposphere/lower stratosphere.

Jet streams and Rossby waves are clearly visible in this image of the wind speeds of the upper troposphere/lower stratosphere.

Download the source code for these images here.

Wave Heights in the Atlantic Ocean

rNOMADS can also access the NOAA Wave Watch Model.  While I’m less familiar with this one, I was able to generate a plot of what (I think) are wave heights in the western Atlantic Ocean and the Caribbean.

Wave heights in the Atlantic ocean from the NOAA Wave Watch model.

Wave heights in the Atlantic ocean from the NOAA Wave Watch model.

Download the source code for this image here.

Instantaneous Atmospheric Profile over Sakura-Jima volcano, Japan

It’s important to know which direction the winds are going directly above active volcanoes, because eruptions can carry ash into air space and over inhabited areas.  One impetus for the continued development of rNOMADS was to provide a one-stop solution for generating high precision wind profiles over specific points, allowing ash distribution prediction as soon as an eruption begins.  Here, we have generated a spatially and temporally interpolated wind and temperature profile over Sakura-Jima volcano, Japan.  The profile is calculated for the exact time when the source code is run.

Instantaneous temperature profile above Sakura-jima volcano, Japan.

Instantaneous temperature profile above Sakura-jima volcano, Japan.

An eruption at this instant would produce ashfalls east of the volcano for a plume height of 15 km.  However, if a truly massive 30 km high plume was produced, ashfalls would occur both east and west of the volcano.

An eruption at this instant would produce ashfalls east of the volcano for a plume height of 15 km. However, if a truly massive 30 km high plume was produced, ashfalls would occur both east and west of the volcano.

Download the source code for these images here.

Wind and temperature profiles over the USArray seismic and infrasound network

The USArray is a massive grid of seismic and infrasound (low frequency sound) sensors currently occupying the Eastern seaboard of the United States.  Since infrasound propagation is strongly affected by temperature and wind, it’s important to know the weather conditions over each station.  Here, we’ve plotted the temperature and wind profiles for every single station in the USArray on the same plot.  The large variations in wind speed will result in very different infrasound propagation depending on where the stations are located.

Temperature profiles above USArray stations.

Temperature profiles above USArray stations.

Wind speed profiles over USArray stations.  It appears that wind speeds vary by over 100 km/hr at certain heights between stations.

Wind speed profiles over USArray stations. It appears that wind speeds vary by over 100 km/hr at certain heights between stations.

Download the source code for these images here.

Some of these scripts require the aqfig package in R to generate the colorbar legends.

Instructions on installing R, rNOMADS, and their dependencies are here.

Trajectory of the Qu8ke Rocket

Back in 2011, the Qu8ke rocket created quite a buzz with this amazing video from its flight:

One of the goals of this rocket was to win the Carmack Prize, an award for the first amateur rocket to obtains a GPS reading from over 100,000 ft. While this rocket did break the 100,000 ft mark, it unfortunately lost GPS lock during its high speed ascent. Analysis of the on board accelerometer data indicates that the rocket reached ~121,000 ft. In a previous post, we describe a method for determining the position a high altitude photograph was taken from using a set of landmarks, OpenCV, and information about the camera’s field of view. The video from the flight of the Qu8ke rocket provides an excellent opportunity to revisit this method and to confirm the accelerometer predicted elevation.

Step 1: Find the GoPro camera parameters

The Qu8ke rocket carried a GoPro camera, which recorded video of the flight in 1080p. To analyze this video with OpenCV, we first needed to find the camera matrix, which contains information about the field of view and radial distortion (fisheye) of the lens.  Following OpenCv’s  method for finding the camera matrix from a set of pictures with a known geometry (we used a chess board), we found that it was difficult to obtain a matrix that could reliably undistort an entire image taken by a GoPro. This may be due to the fact that OpenCV only uses three radial distortion coefficients, which are unable to fully describe the distortion in the image resulting from the 170 degree field of view lens. Despite this, the camera matrix is able to accurately describe mapping between points near the center of the image, as shown below.  We avoided this issue by only picking landmarks not close to the edges/corners of the image.

My beautiful picture

Original (left), Undistorted (Right)

Step 2: Extract a series of stills from Youtube video

We downloaded the source mp4 video from youtube with the 4K video downloader and extracted pictures using ffmpeg. For example the following command:

ffmpef -ss 100 -i qu8ke.mp4 -r 1 -qscale:v 2 %05.jpg

extracts relatively high quality jpegs from the movie every 1 s, starting 100 s into the video. ImageMagick’s convert tool was used to rotate, crop, and resize batches of images.

Step 3: Use OpenCV with a set of landmarks visible in the images to estimate the camera position:

Although the video from the flight is 1080p, the extracted still images aren’t very high quality due to motion blur and the melted plastic that accumulated on the lens (things get hot when you are going Mach 1+). A picture from the rocket near the apogee is shown below:

pnp_1The black circles are the human picked landmark locations, while the blue circles are the reprojected landmark locations given the predicted camera position and orientation. The fact that the circles overlap demonstrate that we have a self-consistent estimate of the camera pose. Adding small amounts of noise to the camera matrix and 2D/3D landmark locations, we are able to estimate the uncertainty in the camera pose. For most of the images we analyzed, the uncertainty in the position was below 6,000 ft. Using pictures throughout the flight like the one above, we reconstructed the trajectory of the rocket:

quake_flightOur predictions confirm that the Qu8ke rocket surpassed 100,000 ft in elevation, and indicate that the rocket went a bit higher (~13,000 ft) than predicted by the on board accelerometers. I am not an expert in interpreting accelerometer data, but I would guess that integrating the acceleration data over the course of the rocket flight becomes increasingly prone to error as time increases, which could account for the difference in the elevation predictions.

Balloon from Floral Arrangement Flies 135 Miles

Ever wondered how far a party balloon can fly if you let it go, and how long it stays in the air?  Turns out they can fly for over a hundred miles and for at least 8 hours or so.  I discovered this when I released an aluminized party balloon from the UNC Chapel Hill Campus with my email address attached.  This afternoon, I got an email from a resident of Pantego, North Carolina – a sleepy town of 200 people located right near the seashore, about 135 miles from Chapel Hill as the crow flies.  She found it in her front yard at about 10 AM today.

This is an example of the type of balloon I released. I think the inside’s coated with aluminum, as opposed to the stretchy ones you blow up with air.

I released it at about 8 AM yesterday morning, so it probably arrived sometime last night.  I’m actually surprised it travelled such a short distance given that it was in the air all day – but it could be that it did not take a direct flight path, stayed low to the ground (where winds are milder) or landed yesterday and she just didn’t see it.

She said that there were no rips in the balloon, so it must have lost its helium in transit.  This is rather surprising, since it sat in my lab for about 2 weeks and only lost about half its volume.  I assume that the increase in internal pressure compared to external pressure as the balloon rises increases the rate of helium loss.

The reason I had this balloon is because I recently passed an exam, so my mom sent me a floral arrangement with a balloon that said “Congratulations!  You did it!”  Initially I was curious to see how long it would take to lose its buoyancy.    Once it was clear that the end was near, I figured it was a pity to waste a perfectly good balloon, so I attached a message and let it go.

Launch: 8 AM, March 11th, at 35.908101 North, -79.052106 West

Landing: Sometime before 10 AM, March 12, at 35.589202 North, -76.662261 West