Land surface model pixels

We begin by taking a look at the spatial information file for the LSM. First we specify the path to the ‘Fourmile_Creek’ test case. We pick the “geo” file, which holds the geospatial information for the LSM. The following commands plot the topographic height for each cell over a map.

geoFile <- paste0(fcPath,'/run.ChannelRouting/DOMAIN/geo_em_d01.Fourmile1km.nlcd11.nc')
coordsProj <- GetDomainCoordsProj(geoFile)
ClosureGeo <- VisualizeDomain(geoFile, plotVar='HGT_M', plot=FALSE, plotDf=coordsProj)

The step of getting coordsProj could be done internally to VisualizeDomain but in some cases (e.g. plotting many variables for a big domain) can speed up the visualization.

For illustration purposes, we set plot=FALSE in the call to VisualizeDomain to skip the default plot output. Typically this first, default plot serves as a basis for tweaking the options to the returned function, named ClosureGeo in this example, which help achieve the desired graphic representation.

Note carefully that VisualizeDomain returns a function. A returned function is termed a closure: a function with data inside it. (Because our style convention is first-letter capitals for functions, ClosureGeo begins with a capital letter.) The returned function can also return objects (perhpas another closure even). The arguments to this closure allow most of the main plot attributes to be tweaked when called. Because the graphics in this function are handled using the ggplot2 package, which has a signifcant learning curve, we design these closures to give the user the 90% solution without th need to learn ggplot2. The return values of the closures give the data and the ggplot objects which can befurther manipulated as necessary, which can be useful if you do know ggplot2.

Let’s examine the structure and the arguments of ClosureGeo and it’s return value, so we know what we have to work with.

class(ClosureGeo)

## [1] "function"

args(ClosureGeo)

## function (location = as.numeric(bbox), zoom = 10, extent = "panel", 
##     source = "google", maptype = "hybrid", pointsize = 3, pointshape = 15, 
##     alpha = 1, gradNColors = c("red", "green", "blue"), subsetRange = NULL, 
##     xlim = NULL, ylim = NULL, plot = TRUE, returnComponents = FALSE) 
## NULL

closureReturnGeo <- ClosureGeo(plot=FALSE)

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

## converting bounding box to center/zoom specification. (experimental)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040369,-105.442619&zoom=10&size=640x640&scale=2&maptype=hybrid&language=en-EN

class(closureReturnGeo)

## [1] "list"

names(closureReturnGeo)

## [1] "plotDf" "ggObj"

We see ClosureGeo is a function and all its arguments. The closureRreturnGeo object is a list (the output of str was much too long to put here) and it has two top-level objects. Note below, we can reference some of these objects in calls to ClosureGeo. In calling ClosureGeo to tweak the plot, zoom and pointsize depend on our domain and on the size of the graphics device. These were adjusted to get the output in the vignette and may not match well to own graphics device.

mapMargin <- .05*c(-1,1)
closureRtnGeo <-
  ClosureGeo(zoom=11, pointsize=9, 
             grad=topo.colors(15), alpha=.6, maptype='terrain',
             subsetRange=range(plotDf$value),
             xlim=range(plotDf$long)+mapMargin*1.5,
             ylim = range(plotDf$lat)+mapMargin)

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

## converting bounding box to center/zoom specification. (experimental)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040369,-105.442619&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN

We could work directly with the data and the plot object from closureRtnGeo instead of using ClosureReturn. The returnComponents argument to the closure function gives even more fine grained output. We’ll illustrate working with the return object below.

Routing pixels

A similar set of function calls allows visualization of the hydrologic grids. We look at the gridded channel network file. This time we skip the call to GetDomainCoordsProj.

hydroFile <- paste0(fcPath,'/run.ChannelRouting/DOMAIN/Fulldom_hires_hydrofile.Fourmile100m.nc')
ClosureHydro <- VisualizeDomain(hydroFile, plotVar='CHANNELGRID', plot=FALSE)
closureRtnHydro <-
  ClosureHydro(zoom=11, pointsize=1,
               ## can reference the internal plotDf (or other variables internal to the closure)
               location=c(lon=mean(plotDf$long), lat=mean(plotDf$lat)), alpha=.2,
               grad=c('white','blue'), maptype='terrain')

## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040415,-105.442563&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN

We highlight use of several of the available arguments to the closure. To see the full set of arguments to the closure, we could use the args function or the print (default) method on the closure by “calling” the function without its paraentheses. Printing the closure will return the entire code and the arguments are shown at the top with their default values (e.g. argument=default). Reading the coode will reveal two variables which actually come from the VisualizeDomain function: bbox and plotDf. (Because the returned function encloses this data, it is refered to as a closure). Several of the arguments to the closure can reference these variables: location, xlim, ylim, subsetRange.

While the above command shows all the hydrologic grid pixels, many of them are missing values (-9999) because they are not on the gridded channel network. We use the subsetRange keyword to limit the range of values represented in the plot. We again use the xlim and ylim keywords to zoom in. We also limit the color scale to blue.

closureRtnHydro <-
  ClosureHydro(zoom=14, pointsize=1,
               ## can reference the internal plotDf (or other variables internal to the closure)
               location=c(lon=mean(plotDf$long), lat=mean(plotDf$lat)),
               subsetRange=c(0),
               grad='blue', maptype='terrain',
               xlim=mean(plotDf$long)+mapMargin/2,
               ylim =mean(plotDf$lat)+mapMargin/2/1.5)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.040415,-105.442563&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN

## Warning: Removed 258 rows containing missing values (geom_point).

To illustrate use of the return values of the above closures, we show how to combine the two plots.

names(closureRtnGeo)

## [1] "plotDf" "ggObj"

names(closureRtnHydro)

## [1] "plotDf" "ggObj"

closureRtnGeo$ggObj +
  geom_point(data=closureRtnHydro$plotDf,
             aes(x=long, y=lat), size=.5, shape=21, fill='white', color='darkblue') +
               ggtitle('Fourmile Creek, CO - Elevation and Stream Channel')

The returnComponents argument to the closures will return more fine-grained objects from the plot in case those are useful.