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 see the vignette here, for a cross platform version of the same using the GrADS-DODS framework, 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_0p50")
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)
variables <- c("ugrd10m", "vgrd10m") #E-W and N-S wind

wind.data <- DODSGrab(latest.model, latest.model.run, variables,
time, lon.inds, lat.inds)

profile <- BuildProfile(wind.data, lon, lat, spatial.average = TRUE, points = 4)

#Present results!
print(paste("At", profile[[1]]$forecast.date, "the East-West winds at Briar Chapel were going", sprintf("%.2f", profile[[1]]$profile.data[1, which(profile[[1]]$variables == "ugrd10m"), 1]),
"meters per second, and the north-south winds were going", sprintf("%.2f", profile[[1]]$profile.data[1, which(profile[[1]]$variables == "vgrd10m"), 1]),
"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_0p50")
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
variables <- c("tmpprs", "hgtprs") #First get temperature

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

#Interpolate to the point of interest
profile <- BuildProfile(model.data, lon, lat,
spatial.average = TRUE, points = 4)

#Plot it!
tmp <- profile[[1]]$profile.data[,which(profile[[1]]$variables == "tmpprs"),] - 272.15
hgt <- profile[[1]]$profile.data[,which(profile[[1]]$variables == "hgtprs"),]

plot(tmp, hgt, xlab = "Temperature (C)",
ylab = "Geopotential Height", main = paste("Temperature above Chapel Hill, NC, at",
profile[[1]]$forecast.date))

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_0p50")
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))
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.

2. Wave heights in the northwest Atlantic ocean

library(rNOMADS)
library(GEOmap)

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.