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

GeoConnex

from pathlib import Path

import geopandas as gpd
import pandas as pd

import pynhd
from pynhd import GeoConnex

GeoConnex is a web service developed by Internet Of Water that provides access to many hydro-linked datasets. PyNHD provides access to all GeoConnex endpoints.

For single queries it’s recommended to use the pynhd.geoconnex function and for multiple queries, it’s more efficient to directly use the underlying pynhd.GeoConnex class that the pynhd.geoconnex function calls under-the-hood.

Let’s start by the geoconnex function. If you run the function without any arguments, it will print out a list of available endpoints. Moreover, if you run the function with item but no query, it will print out the description, queryable fields, and extent of the selected endpoint (item).

pynhd.geoconnex()

For example, let’s select the hu02 item and find all its queryable fields.

pynhd.geoconnex("hu02")

Now that we know all the queryable fields of hu02, we can decide what’s the best way to pull our desired data. For example, let’s get the data for “HUC02” using HUC2 field.

huc2 = pynhd.geoconnex(item="hu02", query={"HUC2": "02"})
ax = huc2.plot(figsize=(6, 6))
ax.set_axis_off()

GeoConnex, also supports spatial query. At the moment, it only accepts requesting features within a bounding box. For example, let’s get all the stations within HUC02 that we previously retrieved.

gages = pynhd.geoconnex(item="gages", query={"geometry": huc2.geometry.iloc[0]})

ax = huc2.plot(facecolor="none", edgecolor="black", figsize=(6, 6))
gages.plot(ax=ax, color="red", marker="*")
ax.set_axis_off()
ax.figure.savefig(Path("_static", "geoconnex.png"), dpi=300, bbox_inches="tight", facecolor="w")

Next, we want to retrieve the available data for a couple of USGS stations using the GeoConnex class.

gcx = GeoConnex()
gcx.item = "gages"  # instead we can instantiate the class like so GeoConnex("gages")
station_ids = ["01031450", "01031500", "01031510"]
gages = gpd.GeoDataFrame(
    pd.concat((gcx.query({"provider_id": sid}) for sid in station_ids), ignore_index=True),
    crs="epsg:4326",
)
gages

Let’s make a bit more advanced spatial query. Assume that we want to get all the stations within a 5-km radius of each the stations that we obtained above. First, we need to buffer the points by 5 km to create a multipolygon. Then we pass it as a query to geoconnex.

buff = gages.to_crs(5070).buffer(5e3).to_crs(4326)
gages = pynhd.geoconnex(item="gages", query={"geometry": buff.unary_union})

ax = buff.plot(facecolor="none", edgecolor="black", figsize=(6, 6))
gages.plot(ax=ax, color="red", marker="*")
ax.set_axis_off()