🚨 This package is under heavy development and breaking changes are likely to happen. 🚨
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 ato_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.
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
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)
Contributions are very welcomed. Please read CODE_OF_CONDUCT.rst and CONTRIBUTING.rst files for instructions.