pynhd 0.1.0

Process and access NLDI and NHDPlus V2 databases from web services.

🚨 This package is under heavy development and breaking changes are likely to happen. 🚨

Features

PyNHD is a part of Hydrodata software stack and provides an interface to access WaterData and NLDI web services. These two web services can be used to navigate and extract vector data from NHDPlus V2 database such as catchments, HUC8, HUC12, GagesII, flowlines, and water bodies.

Additionally, PyNHD offers some extra utilities for processing the flowlines:

• prepare_nhdplus: For cleaning up the dataframe by, for example, removing tiny networks, adding a to_comid column, and finding a terminal flowlines if it doesn’t exist.

• topoogical_sort: For sorting the river network topologically which is useful for routing and flow accumulation.

• vector_accumulation: For computing flow accumulation in a river network. This function is generic and any routing method can be plugged in.

These utilities are developed based on an R package called nhdplusTools. Moreover, requests for additional functionalities can be submitted via issue tracker.

Installation

You can install PyNHD using pip after installing libgdal on your system (for example, in Ubuntu run sudo apt install libgdal-dev):

$ pip install pynhd

Alternatively, PyNHD can be installed from the conda-forge repository using Conda:

$ conda install -c conda-forge pynhd

Quickstart

Let’s explore the capabilities of NLDI. We need to instantiate the class first:

from pynhd import NLDI

nldi = NLDI()
station_id = "USGS-01031500"
UT = "upstreamTributaries"
UM = "upstreamMain"

We can get the basin geometry for the USGS station number 01031500:

basin = nldi.getfeature_byid("nwissite", station_id, basin=True)

NLDI offers navigating a river network from any point in the network in the upstream or downstream direction. We can also limit the navigation distance (in km). The navigation can be done for all the valid NLDI sources which are comid, huc12pp, nwissite, wade, WQP. For example, let’s find all the USGS stations upstream of 01031500, throught the tributaries, and then limit the navigation to only the main channel.

args = {
    "fsource": "nwissite",
    "fid": station_id,
    "navigation": UM,
    "source": "nwissite",
    "distance": None,
}

st_main = nldi.navigate_byid(**args)

args["distance"] = 20  # km
st_d150 = nldi.navigate_byid(**args)

args.update({"distance": None, "navigation": UT})
st_trib = nldi.navigate_byid(**args)

We can set the source to huc12pp to get HUC12 pour points.

args["source"] = "huc12pp"
pp = nldi.navigate_byid(**args)

NLDI provides only comid and geometry of the flowlines which can further be used to get the other available columns in the NHDPlus database. Let’s see how we can combine NLDI and WaterData to get the NHDPlus data for our station.

wd = WaterData("nhdflowline_network")

args.update({"source" : None, "navigation": UM})
comids = nldi.navigate_byid(**args).nhdplus_comid.tolist()
flw_main = wd.byid("comid", comids)

args["navigation"] = UT
comids = nldi.navigate_byid(**args).nhdplus_comid.tolist()
flw_trib = wd.byid("comid", comids)

Other feature sources in the WaterData database are nhdarea, nhdwaterbody, catchmentsp, gagesii, huc08, huc12, huc12agg, and huc12all. For example, we can get the contributing catchments of the flowlines using catchmentsp.

wd = WaterData("catchmentsp")
catchments = wd.byid("featureid", comids)

The WaterData class also has a method called bybox to get data from the feature sources within a bounding box.

wd = WaterData("nhdwaterbody")
wb = wd.bybox((-69.7718, 45.0742, -69.3141, 45.4534))

Next, lets clean up the flowlines and use it to compute flow accumulation. For simplicity, we assume that the flow in each river segment is equal to the length of the segment. Therefore, the accumulated flow at each point should be equal to the sum of the lengths of all its upstream river segments i.e., arbolatesu column in the NHDPlus database. We can use this to validate the flow accumulation result.

import pynhd as nhd

flw = nhd.prepare_nhdplus(flw_trib, 1, 1, 1, True, True)

def routing(qin, q):
    return qin + q

qsim = nhd.vector_accumulation(
    flw[["comid", "tocomid", "lengthkm"]], routing, "lengthkm", ["lengthkm"],
)
flw = flw.merge(qsim, on="comid")
diff = flw.arbolatesu - flw.acc

print(diff.abs().sum() < 1e-5)

Contributing

Contributions are very welcomed. Please read CODE_OF_CONDUCT.rst and CONTRIBUTING.rst files for instructions.