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

#### 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)

#### 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")