Background

When setting up a domain and forecast points, it can be useful to quickly visulaize and interact with domain and gage information. This vignette illustrates basic functionality for visualizing the WRF Hydro domain and channel files and potential forecast points. # Setup Load the rwrfhydro package.

library("rwrfhydro")

## To check rwrfhydro updates run: CheckForUpdates()

VisualizeDomain()

This is the path to the directory of WRF Hydro test cases that we recommend, but you may have to configure this to your machine.

 tcPath <- '~/wrfHydroTestCases/'

Setup: The specific test case and files should be the same once the above is set.

 fcPath <- paste0(tcPath,'Fourmile_Creek/')
 hydroFile<-paste0(fcPath,'/DOMAIN/Fulldom_hydro_OrodellBasin_100m.nc')

The following function plots the selected variable in the hydro file.

 GgMapFunction <- VisualizeDomain(hydroFile, "CHANNELGRID")

## Warning: bounding box given to google - spatial extent only approximate.

## converting bounding box to center/zoom specification. (experimental)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.040812,-105.442574&zoom=10&size=640x640&scale=2&maptype=hybrid&language=en-EN&sensor=false

Note that I know of no reliable function to translate bounding box to a google zoom level, so the default zoom may be very poor for other domains and will need to be set manually…. which we are about to describe how to do.

Notice that VisualizeDomain() returns a function. This is function is a closure, it is a function with data “inside” it (it’s own environment). The idea here is to stash the plot data inside the returned function while letting arguments to the function control aspects of the plot. (Below we also illustrate that the return value from the returned function/closure is actually a ggplot object that can be manipulated outside the closure for finer-grained control of plotting. The data can also be removed from the ggplot object as well.) Now we change the map type, the zoom, and the point attributes with arguments passed to the function.

  ggMap1 <- GgMapFunction(zoom=11, pointshape=15, pointsize=7, 
                         source="google", maptype="terrain")

## Warning: bounding box given to google - spatial extent only approximate.

## converting bounding box to center/zoom specification. (experimental)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.040812,-105.442574&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

For more information on the function returned by VisualizeDomain, see ?VisualizeDomain or to just list its arguments and their default values args(GgMapFunction).

Now add a streamflow gauge points to the domain. These points could be taken from frxst_points.txt or anywhere. Here we compare reality of gage location with that used in the model. This will illustrate an important point about WRF Hydro basins.

orodellLonLat <- data.frame(lon=c(254.6722259521484375, 254.67374999999998408)-360, 
                            lat=c(40.019321441650390625, 40.018666670000001773),
                            gauge=c('model','USGS'))
ggMap2 <- GgMapFunction(location=c(lon=orodellLonLat$lon[1], lat=orodellLonLat$lat[1]),
                         zoom=14, pointshape=15, pointsize=7, 
                         source="google", maptype="terrain", plot=TRUE) 

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.019321,-105.327774&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

Currently, the default plotting in VisualizeDomain always plots the whole domain. (You are free to improve this and do github pull requests on your work!) However, we can zoom in by using the the returned information from ggMapFunction. What is actually in ggMap2 above? (Avoiding the full output of str().)

class(ggMap2)

## [1] "gg"     "ggplot"

It’s a ggplot object that one can use directly with ggplot (instead of through GgMapFunction.) A little of customization with ggplot functions is used to limit the part of the domain plotted.

ggMap2 + geom_point(data=orodellLonLat, aes(x=lon,y=lat, shape=gauge)) +
          scale_x_continuous(limits=rev(orodellLonLat$lon+c( .01,-.01))) +
          scale_y_continuous(limits=rev(orodellLonLat$lat+c( .005,-.005))) 

## Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## Warning in loop_apply(n, do.ply): Removed 14518 rows containing missing values (geom_point).

Now we are zoomed in to the basin outlet. Note that the model point is above the gage and is not the last point in the basin. The gridded channel routing does not supply a flow into the last point in the domain, so this point cannot be used.

VisualizeChanNtwk()

Let’s use look at simulated flows to also see that the flow at the basin outlet is zero. Setup the path to a “CHRTOUT” data file.

 chrtFile <- paste0(fcPath,'/RUN.RTTESTS/OUTPUT_CHRT_DAILY/201308010000.CHRTOUT_DOMAIN1')

This is the basic function which shows the flow on the network

LocLinkFun<- VisualizeChanNtwk(chrtFile)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.037668,-105.437332&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

As with VisualizeDomain, VisualizeChanNtwk returns a closure. You can look at the closure function arguments with

args(LocLinkFun)

## function (location = c(lon = mean(range(linkDf$lon)), lat = mean(range(linkDf$lat))), 
##     zoom = 11, source = "google", maptype = "terrain", padPlot = 0.1, 
##     gaugeZoom = NULL, clickSelect = FALSE, linkShape = 5, gaugeShape = 4) 
## NULL

A main issue is knowing/finding the index of a given point. The click option to the closure lets you click and get the inded. Clicking at the outlet point to see that it has (q=) 0 flow and its index is 350. (We can’t show the interaction in a static document, so you have to try it yourself.)

LocLinkFun(click=TRUE)
## Please click on plot to select nearest point to your click...
## Selected point (in cyan on plot) data:
##  ind            lon           lat q
##  350 -105.325675964 40.0182723999 0 

Using the original function, any set of valid indices can be excluded from the data. After excluding the lowest point and getting a new LocLinkFun, clicking on the lowest point reveals that that index is 1.

LocLinkFun<-VisualizeChanNtwk(chrtFile, exclude=350) 
LocLinkFun(click=TRUE) 
## Please click on plot to select nearest point to your click...
## Selected point (in cyan on plot) data:
## ind            lon           lat               q
##   1 -105.327774048 40.0193214417 0.0453463643789

We want to add the real-life gauge points and find the nearest stream channel grid cells. You can supply gage points in the following format to the VisualizeChanNtwk function and the nearest neighbors on the channel network are returned.

gaugePts <-
   list(orodell   =data.frame(lon=254.67374999999998408,
                              lat=40.018666670000001773),
        loganMill =data.frame(lon=254.63508330000001934,
                              lat=40.042027779999997961),
        sunshine  =data.frame(lon=254.65122220000000652,
                              lat=40.05761110000000258)  )
LocLinkFun <- VisualizeChanNtwk(chrtFile, gaugePts=gaugePts, exc=350, plot=FALSE)

## ** orodell **********
##  system chanInd               q            lon           lat
##   gauge      NA              NA -105.326250000 40.0186666700
##   model       1 0.0453463643789 -105.327774048 40.0193214417
## ** loganMill **********
##  system chanInd               q            lon           lat
##   gauge      NA              NA -105.364916700 40.0420277800
##   model     259 0.0359468720853 -105.364463806 40.0423851013
## ** sunshine **********
##  system chanInd                 q            lon           lat
##   gauge      NA                NA -105.348777800 40.0576111000
##   model     260 0.000123511432321 -105.350837708 40.0423851013

Increase the accuracy of the lon/lat ouput and make the plot that was suppressed in the previous call.

LocLinkFun <- VisualizeChanNtwk(chrtFile, gaugePts=gaugePts, exc=350, plot=FALSE, gaugeAccuracy=17)

## ** orodell **********
##  system chanInd                    q                 lon                lat
##   gauge      NA                   NA -105.32625000000002 40.018666670000002
##   model       1 0.045346364378929138 -105.32777404785156 40.019321441650391
## ** loganMill **********
##  system chanInd                   q                 lon                lat
##   gauge      NA                  NA -105.36491669999998 40.042027779999998
##   model     259 0.03594687208533287 -105.36446380615234 40.042385101318359
## ** sunshine **********
##  system chanInd                      q                 lon                lat
##   gauge      NA                     NA -105.34877779999999 40.057611100000003
##   model     260 0.00012351143232081085 -105.35083770751953 40.042385101318359

LocLinkFun()

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.038193,-105.438381&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

Change the amount of padding around the domain and the shape of the gauge symbols.

 LocLinkFun(pad=.3, gaugeShape=16)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.038193,-105.438381&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

Zoom to the orodell gauge.

LocLinkFun(zoom=14, gaugeShape=16, gaugeZoom='orodell', pad=15)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.018994,-105.327012&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

## Warning in loop_apply(n, do.ply): Removed 327 rows containing missing values (geom_point).

## Warning in loop_apply(n, do.ply): Removed 4 rows containing missing values (geom_point).

Zoom to the Logan Mill gauge

 LocLinkFun(zoom=15, gaugeShape=16, gaugeZoom='loganMill', pad=15)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.042206,-105.36469&zoom=15&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

## Warning in loop_apply(n, do.ply): Removed 331 rows containing missing values (geom_point).

## Warning in loop_apply(n, do.ply): Removed 4 rows containing missing values (geom_point).

You can also click in the zoomed view.

 LocLinkFun(zoom=15, gaugeShape=16, gaugeZoom='loganMill', pad=15, click=TRUE)