In mikejohnson51/NWM: Provides access to the NOAA National Water Model output

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6, 
  fig.height = 4
)

devtools::load_all(".")
library(leaflet)

In August 2016 NOAA made the National Water Model operational aimed towards impact-based forecasting of weather and water events. The NWM runs on NOAA’s centralized Weather and Climate Operational Supercomputing System (WCOSS) utilizing the community WRF-Hydro modeling system. It delivers streamflow forecasts for the 2.7 million USGS NHDPlusv2 river reaches as well as 1-kilometer gridded analyses for a range of hydrologic variables across the CONUS.

The nwm package offers methods to access and subset this data in a quick fashion and returns processed data ready for analysis in the R environment.

What are the challenges?

Given the size of CONUS and the resolution of both point and grid data, the sheer size of output is enormous. Each day, approximately 400 GB of data relating to streamflow, terrain and atmospheric processes is produced. This data is stored for a 48 hour rolling window on the NOAA NOMADS server, and for a 40 day rolling window on the HydroShare THREDDS server hosted at RENCI. These outputs come from four unique model configurations that cycle on different intervals, using different forcing data, and generate different output.

The first challenge this package hopes to address is how people can better to understand and interact with the NWM output.

This is addressed with the look function which prints the meta data and description of each model configuration and output type.

The second challenge is how to this data can be queried and subset to an Area of Interest.

The is addressed through a dependency on the AOI package to define an area of interest; the getFilelist function to define the paths to the desired data sets, and the downloadNWM function to get the data and process it for analysis.

The third challenge is how to make these processes fast, efficient, and useful in the R environment.

This is handeled on the backed of the package making use of the parralization capabilities of foreach and doParallel packages along with the subsetting capabilities of the OPeNDAP framework

This document will walk users through the use of this package, but first we need to install nwm!

Challenge #1:

Understanding the configurations, output types, data, units, and nomenclature.

The NWM runs on four unique configurations - each cycling on a different time intervals, producing forecasts out to different lengths, and including varying suites of parameters.

Specifically these are analysis & assimulation, short range, medium range, and long range. To get a better understanding of these configurations users can call on the look() function. Below is an example for the short_range configuration.

look("short_range")

The output of this call is sent to the R console and provides (1) a description of the configuration, (2) the output types offered, (3) the time when these forcasts are made (t), and (4) how far out - or forward - each goes (f).

For example a forecast made at noon for 3PM would have a t value of 12 and an f value of 3.

At the most general level, each model configuration generates three sets of outputs.

Channel output are point values generated at the outlets of each of the 2.7 million NHDPusV2 reaches in the CONUS

land outputs are gridded data set to a 1 km grid

forcing data are those parameters used to generate the channel and land forecasts, these are also 1 km gridded data.

look("medium_range", 'channel')

Great so now using the look function you have the tools to check the parameters, cycle time and forecast duration of any configuration, type combination - now we can actually get some data!

Challenge #2:

Querying data and subsetting to an AOI

Defining Area of Interest in the nwm package is based on the AOI package which was developed to help find, define, and refine AOI spatial objects. AOIs can be defined by a state, county or clip area generated from a location and bounding box. To learn more about this package please consult the package documentation.

Here is a basic example querying a 100 square mile AOI centered on 'Colorado Springs, Colorado'. The AOI check() function can be used to view its boundaries:

getAOI(clip = list("Colorado Springs", 10, 10)) %>% check()

We can use our AOI to subset and download the NHD reaches:

AOI = getAOI(clip = list("Colorado Springs", 10, 10))

system.time({
  nhd <-getNHD(AOI)
})

head(nhd)

Here we see that that there are 91 flow lines (so 91 points) where channel forecasts can be found from the NWM. To view these flow lines we can chain leaflet::addPolylines() to the check(AOI) call. We can also quickly label each river by its COMID using the popup options, and symbolize each reach by its order using weight.

check(AOI)$m %>% addPolylines(data = nhd, popup = paste0("COMID: ", nhd$comid), weight  = nhd$streamorde)

With an area selected, we next need to define the file list we need to subset. This requires defining a model configuration (config), the date(s), the time of forecast(t), and the time forward from each t, (f).

In this example we are getting the medium range, channel, forecasts made on the July 12, 2018 at 12AM, 6AM, and 12PM UTC time.

Remember the look() function can be useful in paramatrizing this query:

files = getFilelist(config = "medium_range", date = "20180712", type = "channel", t = c(0,6,12), f = 3)

print(files)

With a known AOI, and a list of files we can use the downloadNWM function to specifiy the parameter we want to grab from these files subset to the region. In this basic example

xx = system.time({
  flows <-downloadNWM(AOI, files, "streamflow")
})

This operation took r xx[3] seconds to extract dim(flows$streamflow)[1] flow recods for the AOI>

Alternitvly we could grab the velocity values from the same file list simple by changing the parameter value to 'velocity':

xx = system.time({
  vel <-downloadNWM(AOI, files, "velocity")
})

head(flows$streamflow, 8)

max = flows$streamflow[order(flows$streamflow$streamflow, decreasing = T),]
head(max)

{par(mfrow = c(1,2))
vizFlows(AOI = AOI, data = flows$streamflow, num = 10, max = TRUE)
vizFlows(AOI = AOI, nhd = nhd,  data = flows$streamflow, num = 10, max = FALSE)}

library(RColorBrewer)

files = getFilelist(config = "medium_range", date = "20180711", type = "land", t = 6, f = seq(3,30,3))

system.time({
  d <-downloadNWM(AOI, files, "accet")
})

rasterVis::levelplot(d$accet,
                     main = "National Water Model: Accumulated ET", 
                     names.attr = as.character(nwm::getGridTime(d$accet)),
                     col.regions = colorRampPalette(brewer.pal(9,"YlOrRd")))

check(AOI)$m %>% addRasterImage(x = d$accet[[1]], opacity = .8)

library(ggplot2)
loc = "Garden of the Gods, Colorado Springs"

pt = AOI::geocode(loc, pt = T)

vals = data.frame(date = getGridTime(d$accet), accet = c(raster::extract(d$accet, y =  pt$pt)))

plot(x = vals$date, 
     y = vals$accet, 
     type = 'o', 
     pch = 16, 
     ylab = "ET", 
     xlab = "Time (UTC)",
     main = loc)