Find a place in the network to start

Lets say we have an observation lakeobs1 that was taken at some point. We want to see if there is a lake at that point, and get the information, including ID, of that lake.

Lo and behold, this point is over top of Lake Tahoe. We can get further information about this lake from the mapping layer we used to link. National Hydrography Dataset Highres (NHDH) has some other useful information included.

#>   permanent_      fdate resolution  gnis_id  gnis_name areasqkm elevation
#> 1   44560536 2012-03-12          2 01654975 Lake Tahoe 498.1268      6229
#>        reachcode ftype fcode shape_leng shape_area visibility centroid_x
#> 1 16050101000339   390 39009   1.624393 0.05186653          0   -2037514
#>   centroid_y MATCH_ID
#> 1    2044681 lakeobs1

Now, lets generate map of the linked lake (Lake Tahoe). Using hydrolinks, we don’t need to dig through shapefiles. We can just use the ID we got from linking.

lake_poly = get_shape_by_id(linked_wb$permanent_, dataset = 'nhdh', feature_type = 'waterbody')
#Simple quick viz of polygon
plot(st_geometry(lake_poly), main='Lake Tahoe', col='Dodgerblue')

The polygon returned is fully functional, so it could be used for quantiative analysis of the lake shape or location. We use the sf package throughout for spatial data handling.

Network functionality

The latest and greatest feature allows for the traversal of the hydrologic network. Buildling on what we did above, we can quickly grab very useful information from the hydrologic network.

Upstream traversal

For example, lets grab and plot all inflows into Lake Tahoe. We will start from the lake and traverse up the hydrologic network. We will specify a max traversal distance of 50km to prevent the traversal from exploding (can happen with unlimited bounds going up the network, e.g., Mississippi).


upstream = traverse_flowlines(50, linked_wb$permanent_, direction = 'in')
upstream_shp = get_shape_by_id(upstream$permanent_, dataset = 'nhdh', feature_type = 'flowline')
plot(st_geometry(upstream_shp), col='palegreen')
plot(st_geometry(lake_poly), main='Lake Tahoe', col='Dodgerblue', add=TRUE)

Boom! We have Tahoe and all its input tributaries. Again, the shape data returned are the same as from the underlying hydrologic network datasets, so they can be used in spatial and other analyses.

Note: Upstream traversal can quickly blow up if you are not careful. The river network is very dendritic and so the traversal quickly becomes expontential the further up any stream network you go (unless you are near the headwater streams).

Downstream traversal

We can also do downstream traversal. We’ll start at Lake Tahoe above and see where the water goes. One of the unique aspects of this is we have combined the flow network for both lakes and streams. This means when we do the traversal, we can examine the stream and lake portion of the flow network independently. Below is an example using Tahoe, which, I did not earlier know, doesn’t drain to the ocean.


downstream = traverse_flowlines(2000, linked_wb$permanent_, direction = 'out')
downstream_shp = get_shape_by_id(downstream$permanent_, dataset = 'nhdh', feature_type = 'flowline')
downstream_lk_shp = get_shape_by_id(downstream$permanent_, dataset = 'nhdh', feature_type = 'waterbody')
plot(st_geometry(downstream_shp), col='palegreen')
plot(st_geometry(downstream_lk_shp), main='Lake Tahoe', col='Dodgerblue', add=TRUE)
plot(st_geometry(lake_poly), main='Lake Tahoe', col='Dodgerblue', add=TRUE)

Another Example

Here is another example of Lake Mendota, Wisconsin. This does an upstream traversal. You can see the multiple small impoundments and other ponds/lakes in the upstream network.

id = link_to_waterbodies(43.112449, -89.429409, 'mendota')

hmm = traverse_flowlines(100, id$permanent_, "in")

fls = get_shape_by_id(hmm$permanent_, feature_type = 'flowline', dataset='nhdh')
wbs = get_shape_by_id(hmm$permanent_, feature_type = 'waterbody', dataset='nhdh')

wb   = get_shape_by_id(id$permanent_, feature_type = 'waterbody')

plot(st_geometry(fls), col='green')
plot(st_geometry(wb), col='orange', add=TRUE)
plot(st_geometry(fls), add=TRUE, col='green')
plot(st_geometry(wbs), add=TRUE, col='blue')


111