# Plotting with nhdplusTools

#### dblodgett@usgs.gov

# Plotting with nhdplusTools

The goal of this vignette is to demonstrate a simple and lightweight approach to building maps with NHDPlus data.

## The `plot_nhdplus` function

`plot_nhdplus` is a work in progress. Not all inputs in the function have been implemented as of 11/18/2019 and additional functionality will be added later. Please leave feature requests and issues you find in an issue here.

`plot_nhdplus` is a function that makes getting a simple plot of NHDPlus data as easy as possible. It works with other functions from `nhdplusTools` for identifying and retrieving watershed outlet locations. See the `plot_nhdplus` documentation for more info.

If we pass `plot_nhdplus` a single NWIS site id, `nhdplusTools` uses web services to get data and we get a plot like this:

    plot_nhdplus("05428500")

If we want to add other watersheds, we can use any outlet available from the Network Linked Data Index. See “nldi” functions elsewhere in `nhdplusTools`.

    plot_nhdplus(list(list("nwissite", "USGS-05428500"),
                      list("huc12pp", "070900020602")))

    plot_nhdplus(list(list("nwissite", "USGS-05428500"),
                      list("huc12pp", "070900020602")))

If we don’t know a site id, we can just pass in one or more latitude / longitude locations.

    start_point <- sf::st_as_sf(data.frame(x = -89.36, y = 43.09),
                                coords = c("x", "y"), crs = 4326)

    plot_nhdplus(start_point)

`plot_nhdplus` also allows modification of streamorder (if you have data available locally) and styles. This plot request shows how to get a subset of data for a plot and the range of options. See documentation for more details.

    source(system.file("extdata/sample_data.R", package = "nhdplusTools"))

    plot_nhdplus(list(list("comid", "13293970"),
                      list("nwissite", "USGS-05428500"),
                      list("huc12pp", "070900020603"),
                      list("huc12pp", "070900020602")),
                 streamorder = 2,
                 nhdplus_data = sample_data,
                 plot_config = list(basin = list(lwd = 2),
                                    outlets = list(huc12pp = list(cex = 1.5),
                                                   comid = list(col = "green"))))

We can also plot NHDPlus data without an outlet at all.

    bbox <- sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816, xmax = -89.24681, ymax = 43.17192),
                        crs = "+proj=longlat +datum=WGS84 +no_defs")
    plot_nhdplus(bbox = bbox)

The plots above are all in the EPSG:3857 projection to be compatible with background tiles. Any data added to these plots must be projected to this coordinate system to be added to the plot.

## Getting Data

What follows shows how to use `nhdplusTools` to create plots without the `plot_nhdplus` function. While super convenient, we all know the “easy button” is never quite right, the description below should help get you started.

For this example, we’ll start from an outlet NWIS Site. Note that other options are possible with `discover_nhdplus_id` and `dataRetrieval::get_nldi_sources`.

    library(sf)
    library(nhdplusTools)
    nwissite <- list(featureSource = "nwissite",
                         featureID = "USGS-05428500")

    flowline <- navigate_nldi(nwissite,
                              mode = "upstreamTributaries",
                              data_source = "flowlines")

    nhdplus <- subset_nhdplus(comids = as.integer(flowline$UT$nhdplus_comid),
                              output_file = file.path(work_dir, "nhdplus.gpkg"),
                              nhdplus_data = "download",
                              overwrite = TRUE, return_data = FALSE)
    #> All intersections performed in latitude/longitude.
    #> Reading NHDFlowline_Network
    #> Spherical geometry (s2) switched off
    #> https://api.water.usgs.gov/geoserver/wmadata/ows
    #> <?xml version="1.0" encoding="UTF-8"?>
    #> <wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:gml="http://www.opengis.net/gml" service="WFS" version="1.1.0" outputFormat="application/json" xsi:schemaLocation="http://www.opengis.net/wfs                         http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
    #>   <wfs:Query xmlns:feature="https://api.water.usgs.gov/wmadata" typeName="feature:nhdflowline_network" srsName="EPSG:4269">
    #>     <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
    #>       <ogc:Or>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13293634</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13293750</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13294298</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13294300</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13294304</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13294308</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13294310</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13294312</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13294314</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>         <ogc:PropertyIsEqualTo>
    #>           <ogc:PropertyName>comid</ogc:PropertyName>
    #>           <ogc:Literal>13294368</ogc:Literal>
    #>         </ogc:PropertyIsEqualTo>
    #>       </ogc:Or>
    #>     </ogc:Filter>
    #>   </wfs:Query>
    #> </wfs:GetFeature>
    #> Spherical geometry (s2) switched on
    #> Writing NHDFlowline_Network

    flowline <- read_sf(nhdplus, "NHDFlowline_Network")

    upstream_nwis <- navigate_nldi(nwissite,
                                   mode = "upstreamTributaries",
                                   data_source = "nwissite")

    basin <- get_nldi_basin(nwissite)
    #> https://labs-beta.waterdata.usgs.gov/api/nldi/linked-data/nwissite/USGS-05428500/basin?simplify=true&splitCatchment=false

Now we have a file at the path held in the variable `nhdplus` and three `sf` `data.frame`s with contents that look like:

    st_layers(nhdplus)
    #> Driver: GPKG
    #> Available layers:
    #>            layer_name geometry_type features fields crs_name
    #> 1 NHDFlowline_Network   Line String       10    137    NAD83
    names(flowline)
    #>   [1] "comid"        "fdate"        "resolution"   "gnis_id"      "gnis_name"
    #>   [6] "lengthkm"     "reachcode"    "flowdir"      "wbareacomi"   "ftype"
    #>  [11] "fcode"        "shape_length" "streamleve"   "streamorde"   "streamcalc"
    #>  [16] "fromnode"     "tonode"       "hydroseq"     "levelpathi"   "pathlength"
    #>  [21] "terminalpa"   "arbolatesu"   "divergence"   "startflag"    "terminalfl"
    #>  [26] "dnlevel"      "uplevelpat"   "uphydroseq"   "dnlevelpat"   "dnminorhyd"
    #>  [31] "dndraincou"   "dnhydroseq"   "frommeas"     "tomeas"       "rtndiv"
    #>  [36] "vpuin"        "vpuout"       "areasqkm"     "totdasqkm"    "divdasqkm"
    #>  [41] "tidal"        "totma"        "wbareatype"   "pathtimema"   "hwnodesqkm"
    #>  [46] "maxelevraw"   "minelevraw"   "maxelevsmo"   "minelevsmo"   "slope"
    #>  [51] "elevfixed"    "hwtype"       "slopelenkm"   "qa_ma"        "va_ma"
    #>  [56] "qc_ma"        "vc_ma"        "qe_ma"        "ve_ma"        "qa_01"
    #>  [61] "va_01"        "qc_01"        "vc_01"        "qe_01"        "ve_01"
    #>  [66] "qa_02"        "va_02"        "qc_02"        "vc_02"        "qe_02"
    #>  [71] "ve_02"        "qa_03"        "va_03"        "qc_03"        "vc_03"
    #>  [76] "qe_03"        "ve_03"        "qa_04"        "va_04"        "qc_04"
    #>  [81] "vc_04"        "qe_04"        "ve_04"        "qa_05"        "va_05"
    #>  [86] "qc_05"        "vc_05"        "qe_05"        "ve_05"        "qa_06"
    #>  [91] "va_06"        "qc_06"        "vc_06"        "qe_06"        "ve_06"
    #>  [96] "qa_07"        "va_07"        "qc_07"        "vc_07"        "qe_07"
    #> [101] "ve_07"        "qa_08"        "va_08"        "qc_08"        "vc_08"
    #> [106] "qe_08"        "ve_08"        "qa_09"        "va_09"        "qc_09"
    #> [111] "vc_09"        "qe_09"        "ve_09"        "qa_10"        "va_10"
    #> [116] "qc_10"        "vc_10"        "qe_10"        "ve_10"        "qa_11"
    #> [121] "va_11"        "qc_11"        "vc_11"        "qe_11"        "ve_11"
    #> [126] "qa_12"        "va_12"        "qc_12"        "vc_12"        "qe_12"
    #> [131] "ve_12"        "lakefract"    "surfarea"     "rareahload"   "rpuid"
    #> [136] "vpuid"        "enabled"      "geom"
    names(upstream_nwis)
    #> [1] "origin"      "UT_nwissite"
    names(basin)
    #> [1] "geometry"
    class(st_geometry(flowline))
    #> [1] "sfc_LINESTRING" "sfc"
    class(st_geometry(upstream_nwis$UT_nwissite))
    #> [1] "sfc_POINT" "sfc"
    class(st_geometry(basin))
    #> [1] "sfc_POLYGON" "sfc"

Our file has four layers: network flowlines, simplified catchments, nhd area features, and nhd waterbody features.

The flowlines have a large set of attributes from the NHDPlus dataset. And the nwis sites have a few attributes that came from the NLDI. Attributes for NWIS sites can be found using the `dataRetrieval` package.

See the NHDPlus user guide linked here for more on what these layers are and what the flowline attributes entail.

## Base R Plotting

In order to maximize flexibility and make sure we understand what’s going on with coordinate reference systems, the demonstration below shows how to use base R plotting with the package `mapsf` and `maptiles`.

In this example, we have to plot just the geometry, extracted with `st_geometry` and we need to project the geometry into the plotting coordinate reference system, EPSG:3857 also known as “web mercator”. The reason we have to make this transformation is that practically all basemap tiles are in this projection and reprojection of pre-rendered tiles doesn’t look good. We do this with a simple `prep_layer` function.

The `mapsf`, `maptiles`, and base `plot` commands put data onto the R plotting device so the first to be plotted is on the bottom. A couple hints here. `lwd` is line width. `pch` is point style. `cex` is an expansion factor. Colors shown below are basic R colors. the `rgb` function is handy for creating colors with transparency if that’s of interest.

    prep_layer <- function(x) st_geometry(st_transform(x, 3857))

    bb <- sf::st_as_sfc(sf::st_bbox(prep_layer(basin)))

    tiles <- maptiles::get_tiles(bb,
                                 zoom = 11, crop = FALSE,
                                 verbose = FALSE,
                                 provider = "Esri.NatGeoWorldMap")

    mapsf::mf_map(bb, type = "base", col = NA, border = NA)

    maptiles::plot_tiles(tiles, add = TRUE)

    mapsf::mf_map(bb, type = "base", add = TRUE, col = NA, border = NA)
    mapsf::mf_arrow(adjust = bb)
    #> Warning: 'adjust' is deprecated.
    #> Use 'align' instead.
    mapsf::mf_scale()

    plot(prep_layer(basin),
         lwd = 2, add = TRUE)

    plot(prep_layer(flowline),
         lwd = 1.5, col = "deepskyblue", add = TRUE)

    plot(prep_layer(dplyr::filter(flowline, streamorde > 2)),
         lwd = 3, col = "darkblue", add = TRUE)

    us_nwis_layer <- prep_layer(upstream_nwis$UT_nwissite)

    plot(us_nwis_layer,
         pch = 17, cex = 1.5, col = "yellow", add = TRUE)

    label_pos <- st_coordinates(us_nwis_layer)

    text(label_pos[,1],label_pos[,2],
         upstream_nwis$identifier,
         adj = c(-0.2, 0.5), cex = 0.7)

## Plotting with ggplot2

Below is a very similar example using `ggmap` and `ggplot2` `geom_sf`. Note that ggmap takes case of projections for us, which should either make you happy because it just works or very nervous because it just works.

    library(ggmap)
    library(ggplot2)

    ggmap_bbox <- setNames(sf::st_bbox(basin), c("left", "bottom", "right", "top"))
    ggmap_bbox
    #>      left    bottom     right       top
    #> -89.60465  43.03507 -89.20378  43.36607

    upstream_nwis <- dplyr::bind_cols(upstream_nwis$UT_nwissite,
                               dplyr::rename(dplyr::as_tibble(sf::st_coordinates(upstream_nwis$UT_nwissite)),
                                             lat = Y, lon = X))

    # ggmap now requires api keys
    # basemap_toner <- get_map(source = "stamen", maptype = "toner",
    #                          location = ggmap_bbox, zoom = 11, messaging = FALSE)
    # basemap_terrain <- get_map(source = "stamen", maptype = "terrain-lines",
    #                            location = ggmap_bbox, zoom = 11, messaging = FALSE)
    # toner_map <- ggmap(basemap_toner)
    # terrain_map <- ggmap(basemap_terrain)
    #
    # toner_map

    ggplot() + geom_sf(data = basin,
                            inherit.aes = FALSE,
                            color = "black", fill = NA) +
      geom_sf(data = flowline,
              inherit.aes = FALSE,
              color = "deepskyblue") +
      geom_sf(data = dplyr::filter(flowline, streamorde > 2),
              inherit.aes = FALSE,
              color = "darkblue") +
      geom_sf(data = upstream_nwis, inherit.aes = FALSE, color = "red") +
      geom_text(data = upstream_nwis, aes(label = identifier, x = lon, y = lat),
                hjust = 0, size=2.5, nudge_x = 0.02, col = "black")

Hopefully these examples give a good head start to plotting NHDPlus data. Please submit questions via github issues for more!! Pull requests on this vignette are more than welcome if you have additions or suggestions.