navass11 / hydrodata Goto Github PK

Hydrodata is a python library designed to aid in watershed analysis. It provides easy and consistent access to a handful of hydrology and climatology databases with some helper functions for visualization. Currently, the following data retrieval services are supported:

• NLDI and NHDPlus V2 for vector river network, catchments, and other NHDPlus data.
• Daymet for climatology data, both single pixel and gridded
• SSEBop for daily actual evapotranspiration, both single pixel and gridded
• NLCD 2016 for land cover, land use (some utilities are available for analysing and plotting the cover data)
• NWIS for daily streamflow observations
• HCDN 2009 for identifying sites where human activity affects the natural flow of the watercourse
• 3DEP from National Map service for getting data such as Digital Elevation Model, slope, and aspect.

Additionally, the following functionalities are offered:

• Interactive map for exploring USGS stations within a bounding box,
• Efficient vector-based flow accumulation in a stream network,
• Computing Potential Evapotranspiration (PET) using Daymet data based on FAO-56,
• High level APIs for easy access to all ArcGIS RESTful-based services as well as WMS- and WFS-based services,
• Helpers for plotting land cover data based on official NLCD cover legends,
• A roughness coefficients lookup table for each land cover type which is useful for overland flow routing.

Requests for additional databases or functionalities can be submitted via issue tracker.

Learn more about Hydrodata in its official documentation at https://hydrodata.readthedocs.io.

You can install Hydrodata using pip after installing libgdal (for example libgdal in Conda environment or libgdal-dev in Ubuntu) on your system or environment:

$ pip install hydrodata

With just a few lines of code, Hydrodata provides easy access to a handful of databases. We can start by exploring the available USGS stations within a bounding box:

import hydrodata.datasets as hds

hds.interactive_map([-70, 44, -69, 46])

Then, we can either specify a station ID or coordinates to the Station function and gathers the USGS site information such as name, contributing drainage area, and watershed geometry.

from hydrodata import Station

wshed = Station(coords=(-69.32, 45.17), dates=('2000-01-01', '2010-01-21'))

The generated wshed object has a property that shows whether the station is in HCDN database i.e., whether it's a natural watershed or is affected by human activity. For this watershed wshed.hcdn is True, therefore, this is a natural watershed. Moreover, using the retrieved information, datasets module provides access to other databases. For example, we can get the main river channel and the tributaries of the watershed, the USGS stations upstream (or downstream) of the main river channel (or tributatires) up to a certain distance, say 150 km or all the stations:

from hydrodata import NLDI

tributaries = NLDI.tributaries(wshed.station_id)
main = NLDI.main(wshed.station_id)
stations = NLDI.stations(wshed.station_id)
stations_m150 = NLDI.stations(wshed.station_id, navigation="upstreamMain", distance=150)

For demonstrating the flow accumulation function, lets assume the flow in each river segment is equal to the length of the river segment. Therefore, it should produce the same results as the arbolatesu variable in the NHDPlus database.

from hydrodata import utils

flw = utils.prepare_nhdplus(NLDI.flowlines('11092450'), 0, 0, purge_non_dendritic=False)

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

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

We can check the validity of the results using diff.abs().sum() = 5e-14. Furthermore, DEM, slope, and aspect can be retrieved for the station's contributing watershed at 30 arc-second (~1 km) resolution as follows:

from hydrodata import NationalMap

nm = NationalMap(wshed.geometry, resolution=30)
dem, slope, aspect = nm.get_dem(), nm.get_slope(), nm.get_aspect()

The climate data and streamflow observations for a location of interest can be retrieved as well. Note the use of pet flag for computing PET:

variables = ["tmin", "tmax", "prcp"]
clm_p = hds.daymet_byloc(wshed.lon, wshed.lat,
                         start=wshed.start, end=wshed.end,
                         variables=variables, pet=True)
clm_p['Q (cms)'] = hds.nwis_streamflow(wshed.station_id, wshed.start, wshed.end)

Other than point-based data, we can get data from gridded databases. The retrieved data are masked with the watershed geometry:

clm_g = hds.daymet_bygeom(wshed.geometry,
                          start='2005-01-01', end='2005-01-31',
                          variables=variables, pet=True)
eta_g = hds.ssebopeta_bygeom(wshed.geometry, start='2005-01-01', end='2005-01-31')

All the gridded data are returned as xarray datasets that has efficient data processing tools. Additionally, Hydrodata has a plot module that plots five hydrologic signatures graphs in one plot:

from hydrodata import plot

plot.signatures(clm_loc['Q (cms)'], wshed.drainage_area, prcp=clm_loc['prcp (mm/day)'], title=wshed.name)

Some example plots are shown below:

The services module can be used to access some other web services as well. For example, we can access Los Angeles GeoHub RESTful service, NationalMap's 3D Eleveation Program via WMS and FEMA National Flood Hazard Layer via WFS as follows:

from hydrodata import ArcGISREST, WFS, services
import geopandas as gpd

la_wshed = Station('2005-01-01', '2005-01-31', '11092450')

url_rest = "https://maps.lacity.org/lahub/rest/services/Stormwater_Information/MapServer/10"
s = ArcGISREST(url_rest, outFormat="json")
s.get_featureids(la_wshed.geometry)
storm_pipes = s.get_features()

url_wms = "https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer"
hillshade = services.wms_bygeom(
    url_wms,
    geometry=wshed.geometry,
    version="1.3.0",
    layers={"hillshade": "3DEPElevation:GreyHillshade_elevationFill"},
    outFormat="image/tiff",
    resolution=1
)

url_wfs = "https://hazards.fema.gov/gis/nfhl/services/public/NFHL/MapServer/WFSServer"
wfs = WFS(
    url_wfs,
    layer="public_NFHL:Base_Flood_Elevations",
    outFormat="esrigeojson",
    crs="epsg:4269",
)
r = wfs.getfeature_bybox(wshed.geometry.bounds, in_crs="epsg:4326")
flood = utils.json_togeodf(r.json(), "epsg:4269", "epsg:4326")

Hydrodata offers some limited statistical analysis. It could be more useful to the watershed modeling community to integrate more data exploratory capabilities to the package. Additionally, adding support for more databases such as water quality, phenology, and water level, are very welcome. If you are interested please get in touch. You can find information about contributing to hydrodata at our Contributing page.

This package was created with Cookiecutter and the audreyr/cookiecutter-pypackage project template.