This page was generated from nhd_demo.ipynb. Interactive online version:

NHD Data using HyRiver#

Author	Affiliation	Email
Taher Chegini	Purdue University	tchegini@purdue.edu
Dave Blodgett	USGS	dblodgett@usgs.gov

This notebook presents capabilities of HyRiver for accessing National Hydrography Dataset (NHDPlus MR and HR) and Watershed Boundary Dataset (WBD). For this purpose, we use the following web services:

NLDI (Network Linked Data Index)
Water Data
3DHP (3D Hydrography Program)
GeoConnex
NHDPlusHR

[1]:

import folium

from pygeohydro import NWIS, WBD
from pynhd import HP3D, NLDI, GeoConnex, NHDPlusHR, WaterData

We select the Wolf River at Langlade, WI (04074950) station for this demonstration.

[2]:

site_id = "04074950"

We start by instantiating the class for these web services. Note that each web service, we usually have to select a desired “layer”. Docstrings of these classes provide more information about the available layers.

[3]:

nldi = NLDI()
nhd_mr = WaterData("nhdflowline_network")
h4_wd = WaterData("wbd04")
h4_wbd = WBD("huc4")
nhd_hr = NHDPlusHR("flowline")
nwis = NWIS()
hp3d = HP3D("flowline")

We can use NLDI to get information about the station and also navigate the NHD MR network up to a certain distance and get the associated NHDPlus MR features.

[4]:

site_feature = nldi.getfeature_byid("nwissite", f"USGS-{site_id}")
upstream_network = nldi.navigate_byid(
    "nwissite", f"USGS-{site_id}", "upstreamMain", "flowlines", distance=9999
)

We use Folium to visualize the station and the network.

[5]:

m = upstream_network.explore()
folium.GeoJson(site_feature, tooltip=folium.GeoJsonTooltip(["identifier"])).add_to(m)
m

[5]:

Make this Notebook Trusted to load map: File -> Trust Notebook

We can also get the drainage basin for this station using NLDI. Additionally, instead of getting the flowlines that are only upstream of the station, we can use the basin’s geometry to obtain all flowlines that are within the bounds of the basin.

[6]:

basin = nldi.get_basins(site_id)
subset = nhd_mr.bygeom(basin.geometry.iloc[0], basin.crs)

[7]:

m = basin.explore(style_kwds={"fillColor": "gray"})
folium.GeoJson(subset, style_function=lambda _: {"color": "blue"}).add_to(m)
m

[7]:

Make this Notebook Trusted to load map: File -> Trust Notebook

Now, let’s use the NHDPlusHR web srvice from the National Map to the HR flowlines instead of the MR flowlines. As expected, the HR flowlines are more detailed than the MR flowlines.

[8]:

flw_hr = nhd_hr.bygeom(basin.geometry.iloc[0], basin.crs)

[9]:

m = basin.explore(style_kwds={"fillColor": "gray"})
folium.GeoJson(flw_hr, style_function=lambda _: {"color": "blue"}).add_to(m)
m

[9]:

Make this Notebook Trusted to load map: File -> Trust Notebook

We can the HR flowline that the station is on by getting all flowlines that are within a certain distance from the station, say 2 km, then use the total drainage area property of the NHDPlus HR to find the flowline that has the closest drainage area to the station’s. So, first we need to get the drainage area of the basin either by computing the drainage area from the baisn’s geometry or by using NWIS web service. Let’s compare the two approaches.

[10]:

site_info = nwis.get_info({"site": site_id}, expanded=True)
sqmi_to_sqkm = 2.59
area_sqkm = site_info["drain_area_va"].iloc[0] * sqmi_to_sqkm

We can see that they are close.

[11]:

area_sqkm, basin.to_crs(5070).area.iloc[0] * 1e-6

[11]:

(1199.1699999999998, 1205.2133095225754)

[12]:

potential_matches = nhd_hr.bygeom(site_info.geometry.iloc[0], site_info.crs, distance=2000)

[13]:

potential_matches.totdasqkm

[13]:

0    1147.879401
1       5.952500
Name: totdasqkm, dtype: float64

[14]:

match = potential_matches.iloc[[potential_matches.totdasqkm.sub(area_sqkm).abs().idxmin()]]

[15]:

m = match.explore()
folium.GeoJson(site_info, tooltip=folium.GeoJsonTooltip(["site_no"])).add_to(m)
m

[15]:

Make this Notebook Trusted to load map: File -> Trust Notebook

Next, we get the 3DHP flowlines within the basin and compare it with the NHDPlus HR flowlines.

[16]:

flw_3dhp = hp3d.bygeom(basin.unary_union, basin.crs)

We can see that there are some flowlines in 3DHP that are not in NHDPLus HR.

[17]:

m = basin.explore(style_kwds={"fillColor": "gray"})
folium.GeoJson(flw_3dhp, style_function=lambda _: {"color": "blue"}).add_to(m)
folium.GeoJson(flw_hr, style_function=lambda _: {"color": "red"}).add_to(m)
m

[17]:

Make this Notebook Trusted to load map: File -> Trust Notebook

GeoConnex provides access to similar datasets and some additional datasets. First, we instatiat its class but don’t provide any layer, so we can later on switch between layers, without having to create a new instance.

[18]:

gcx = GeoConnex()

If you print this object, you can see the available layers and their descriptions.

[19]:

gcx

[19]:

Available Endpoints:
'hu02': Two-digit Hydrologic Regions from USGS NHDPlus High Resolution
'hu04': Four-digit Hydrologic Subregion from USGS NHDPlus High Resolution
'hu06': Six-digit Hydrologic Basins from USGS NHDPlus High Resolution
'hu08': Eight-digit Hydrologic Subbasins from USGS NHDPlus High Resolution
'hu10': Ten-digit Watersheds from USGS NHDPlus High Resolution
'nat_aq': National Aquifers of the United States from USGS National Water Information System National Aquifer code list
'principal_aq': Principal Aquifers of the United States from 2003 USGS data release
'sec_hydrg_reg': Secondary Hydrogeologic Regions of the Conterminous United States from 2018 USGS data release
'gages': United States community contributed reference Stream Gage Monitoring Locations
'mainstems': United States community ontributed reference Mainstem Rivers
'dams': United States Community Contributed Reference Dams
'pws': Public Water Systems from United States EPA Safe Drinking Water Information System
'states': States from United States Census Bureau Cartographic Boundaries
'counties': Counties from United States Census Bureau Cartographic Boundaries
'aiannh': Federally recognized American Indian/Alaska Native Areas/Hawaiian Home Lands (AIANNH)
'cbsa': United States Metropolitan and Micropolitan Statistical Areas
'ua10': United States Urbanized Areas and Urban Clusters from the 2010 Census
'places': United States legally incorporated and Census designated places

Let’s get the gages and mainstems within the bounds of the basin.

[20]:

gcx.item = "gages"
wolf_gages = gcx.bygeometry(basin.unary_union, crs=basin.crs)
gcx.item = "mainstems"
wolf_mainstems = gcx.bygeometry(basin.unary_union, crs=basin.crs)

[21]:

m = basin.explore(style_kwds={"fillColor": "gray"})
folium.GeoJson(wolf_gages, tooltip=folium.GeoJsonTooltip(["provider_id"])).add_to(m)
folium.GeoJson(wolf_mainstems, style_function=lambda _: {"color": "blue"}).add_to(m)
m

[21]:

Make this Notebook Trusted to load map: File -> Trust Notebook

For retrieving the WBD, we have several options that give similar results. We can use GeoConnex, WaterData, or WBD classes. Let’s use the WBD class.

[22]:

wolf_huc4 = h4_wd.bygeom(basin.unary_union, basin.crs)

[23]:

m = wolf_huc4.explore(style_kwds={"fillColor": "gray"})
folium.GeoJson(basin, style_function=lambda _: {"color": "blue"}).add_to(m)
m

[23]:

Make this Notebook Trusted to load map: File -> Trust Notebook