1 Background

USGS streamflow observations are a primary source of hydrologic information and often used for validation and calibration of hydrlogic models. Web services have been developed at NWIS and the dataRetrieval R package has emerged to make it easy to get USGS data using R.

This vignette demonstrates how to download USGS streamflow observations using the dataRetrieval package, import modeled streamflow from WRF-Hydro output files using functionality from rwrfhydro, and plot both the observed and modeled time series. Note that there are different ways to read the simulated streamflow depending on which option we have turned on in the hdyro.namelist.

1.1 Load packages and scripts

Load required packages and scripts.

> library(dataRetrieval)
> library(rwrfhydro)
To check rwrfhydro updates run: CheckForUpdates()
> library(data.table)
> library(ggplot2)
> library(foreach)

2 Download USGS streamflow observations

Download daily USGS streamflow observations using a function from the dataRetrieval package. And the feature_id to it.

> gageID <- c("01374559", "01374581", "0137462010")
> obsDF <- readNWISuv(siteNumbers=gageID,
+                     parameterCd="00060",
+                     startDate="2011-08-25",
+                     endDate="2011-09-03")

Rename columns and create a column for streamflow with units of cubic meters per second.

> colnames(obsDF) <- c("agency","site_no","dateTime","q_cfs","qc", "tz_cd")
> 
> obsDF$dateTime <- as.POSIXct(obsDF$dateTime) 
> obsDF$streamflow <- obsDF$q_cfs/35.31

Read in the reference file, and add the feature_id to the observations

> referenceFile <- "~/wrf-hydro-training/example_case/Gridded/croton_frxst_pts.csv"
> refDF <- read.csv(referenceFile, colClasses = c("integer", "numeric", "numeric", "character", "character"))
> print(refDF)
  FID       LON      LAT STATION    Site_No
1   1 -73.75900 41.46811       A   01374559
2   2 -73.73707 41.44807       B   01374581
3   3 -73.69425 41.41023       C 0137462010
> obsDF <- merge(obsDF,  refDF[, c("FID", "Site_No")], by.x = "site_no", by.y = "Site_No")
> obsDF$feature_id <- obsDF$FID

3 Import WRF-Hydro simulated streamflow at the gage location

One can read the simulated streamflow values different ways depending on what output files are turned on in the hydro.namelist. The easiest option is to read it from frxst_pts_out.txt file.

> # specify the path to the frxst_pts_out.txt file
> modelOutputPath <- "~/wrf-hydro-training/output/lesson4/run_gridded_baseline/frxst_pts_out.txt"
> 
> simQ <- read.csv(modelOutputPath, 
+                  col.names =  c("TimeFromStart", "Time", "station", "longitude", "latitude", "q_cms", "q_cfs", "stage"),
+                  colClasses = c("integer", "character", "character", rep("numeric", 5)), header = FALSE)
> head(simQ)
  TimeFromStart                Time      station longitude latitude q_cms
1          3600 2011-08-26 01:00:00            3 -73.69508 41.41042 2.278
2          3600 2011-08-26 01:00:00            2 -73.73784 41.44910 1.340
3          3600 2011-08-26 01:00:00            1 -73.75716 41.46555 0.699
4          7200 2011-08-26 02:00:00            3 -73.69508 41.41042 2.278
5          7200 2011-08-26 02:00:00            2 -73.73784 41.44910 1.328
6          7200 2011-08-26 02:00:00            1 -73.75716 41.46555 0.693
   q_cfs stage
1 80.440 0.400
2 47.333 0.467
3 24.690 0.437
4 80.441 0.400
5 46.882 0.465
6 24.490 0.435

The other option is to read it from the CHANOBS files like what you did in the training lessons.

> chnObsFiles <- list.files("~/wrf-hydro-training/output/lesson4/run_gridded_baseline/",
+                           glob2rx("2011*CHANOBS*"), full.names = TRUE)
> 
> simQ <- foreach(f = chnObsFiles, .combine = rbind.data.frame) %do% {
+   a <- rwrfhydro::GetNcdfFile(f, variables = c("feature_id", "streamflow"), quiet = TRUE)
+   a$dateTime <- as.POSIXct(substr(basename(f), 1, 12), format = "%Y%m%d%H%M", tz = "UTC")
+   return(a)
+ }

You could also read this information from the CHRTOUT_DOMAIN files and CHRTOUT_GRID, refer to the manual page of ReadChrtout

3.1 Plot hydrographs

Plot observed and simulated hydrographs using ggplot2.

> obsDF$run <- "Observation"
> simQ$run <- "Gridded Baseline"
> 
> merged <- rbind(simQ, obsDF[, names(simQ)])
> 
> ggplot(data = merged) + geom_line(aes(dateTime, streamflow, color = run)) + facet_wrap(~feature_id, ncol = 1)

4 Import WRF-Hydro simulated streamflow at the gage location for NWM

The simulated streamflow can be found in CHANOBS, CHRTOUT_DOMAN and CHRTOUT_GRID files. Here is an example on how to read it from CHANOBS files.

> chnObsFiles <- list.files("~/wrf-hydro-training/output/lesson5/run_nwm",
+                           glob2rx("2011*CHANOBS*"), full.names = TRUE)
> 
> simQnwm <- foreach(f = chnObsFiles, .combine = rbind.data.frame) %do% {
+   a <- rwrfhydro::GetNcdfFile(f, variables = c("feature_id", "streamflow"), quiet = TRUE)
+   a$dateTime <- as.POSIXct(substr(basename(f), 1, 12), format = "%Y%m%d%H%M", tz = "UTC")
+   return(a)
+ }

The relationship between the links/feature_id and the gage ID can be found in the routelink file.

> routelinkFile <- "~/wrf-hydro-training/example_case/NWM/DOMAIN/Route_Link.nc"
> 
> rtlink <- rwrfhydro::GetNcdfFile(routelinkFile, quiet = TRUE)
> head(rtlink)
  BtmWdth ChSlp Kchan Length MusK MusX NHDWaterbodyComID Qi    So    alt
1       2  0.05     0   2203 3600  0.2             -9999  0 0.070 255.40
2       2  0.05     0   2248 3600  0.2             -9999  0 0.076 271.49
3       2  0.05     0   1943 3600  0.2             -9999  0 0.002  69.14
4       2  0.05     0   2233 3600  0.2             -9999  0 0.089 283.30
5       2  0.05     0   2339 3600  0.2             -9999  0 0.028 159.00
6       2  0.05     0   4445 3600  0.2             -9999  0 0.005 116.03
  ascendingIndex from           gages      lat    link       lon    n
1            183    0                 41.55200 6212218 -73.72379 0.06
2            155    0                 41.55030 6212238 -73.73109 0.06
3            184    0                 41.54454 6212242 -73.83585 0.06
4            166    0                 41.51133 6212272 -73.78367 0.06
5            156    0                 41.50933 6212276 -73.79778 0.06
6            168    0                 41.54963 6212528 -73.75257 0.06
  order      to
1     1       0
2     1       0
3     1       0
4     1 6212262
5     1 6212262
6     1 6212228
> # let s subset to only those that have a gage on them
> rtlink <- subset(rtlink, trimws(gages) != "")
> rtlink
     BtmWdth ChSlp Kchan Length MusK MusX NHDWaterbodyComID Qi    So
200 2.286175  0.05     0    397 3600  0.2             -9999  0 0.001
221 5.500000  0.05     0    505 3600  0.2             -9999  0 0.002
240 5.500000  0.05     0   3166 3600  0.2             -9999  0 0.015
       alt ascendingIndex from           gages      lat    link       lon
200 180.48             53    0        01374559 41.47079 6226948 -73.76059
221 164.32             58    0        01374581 41.44981 6227008 -73.73565
240 147.61             77    0      0137462010 41.40192 6227150 -73.68741
        n order      to
200 0.055     3 6228462
221 0.055     4 6227012
240 0.055     4 6227160

Now you can merge the name of the gage to th simulated streamflow.

> simQnwm <- merge(simQnwm, rtlink[, c("link", "gages")], by.x = "feature_id", by.y = "link")

Concatenate the observation and plot it.

> names(simQnwm)[which(names(simQnwm) == "gages")] <- "site_no"
> simQnwm$run <- "NWM"
> 
> merged <- rbind(simQnwm, obsDF[, names(simQnwm)])
> 
> ggplot(data = merged) + geom_line(aes(dateTime, streamflow, color = run)) + facet_wrap(~trimws(site_no), ncol = 1)