# GeoConnex#

[ ]:

from pathlib import Path

import geopandas as gpd
import pandas as pd

import pynhd
from pynhd import GeoConnex


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()