This page was generated from pygeoapi.ipynb.
Interactive online version:
Split Catchment with PyGeoAPI#
[1]:
import geopandas as gpd
import shapely.geometry as sgeom
import pynhd
PyGeoAPI service provides four functionalities:
flow_trace
: Trace flow from a starting point to up/downstream direction.split_catchment
: Split the local catchment of a point of interest at the point’s location.elevation_profile
: Extract elevation profile along a flow path between two points.cross_section
: Extract cross-section at a point of interest along a flow line.
The pygeoapi
function in PyNHD requires two inputs, a geopandas.GeoDataFrame
that must contain all the required inputs corresponding to target services. Let’s take a look at them in an example:
[2]:
gdf = gpd.GeoDataFrame(
{
"direction": [
"none",
]
},
geometry=[sgeom.Point((1774209.63, 856381.68))],
crs="ESRI:102003",
)
trace = pynhd.pygeoapi(gdf, "flow_trace")
In the split_catchment
service we can set upstream
to True
to get the entire upstream catchments not just the split local catchments.
[3]:
gdf = gpd.GeoDataFrame(
{
"upstream": [
False,
]
},
geometry=[sgeom.Point((-73.82705, 43.29139))],
crs="epsg:4326",
)
split = pynhd.pygeoapi(gdf, "split_catchment")
[4]:
ax = split.plot(figsize=(8, 8), facecolor="lightgrey", edgecolor="black", alpha=0.8)
trace.plot(ax=ax, color="b", linewidth=3.0)
ax.axis("off")
ax.set_title("Split Catchment")
ax.margins(x=0)
ax.figure.set_dpi(100)
ax.figure.savefig("_static/split_catchment.png", bbox_inches="tight", facecolor="w", dpi=100)

The elevation_profile
service returns the elevation profile along a flow path between two points at a given number of points and a specified resolution for DEM.
[5]:
coords = [(-103.801086, 40.26772), (-103.80097, 40.270568)]
gdf = gpd.GeoDataFrame(
{
"numpts": [
101,
],
"dem_res": [
1,
],
},
geometry=[sgeom.MultiPoint(coords)],
crs="epsg:4326",
)
profile = pynhd.pygeoapi(gdf, "elevation_profile")
[6]:
ax = profile[["distance", "elevation"]].plot(
x="distance", y="elevation", color="b", linewidth=3.0, figsize=(12, 4)
)
ax.set_ylim(1280, 1320)
ax.legend(["Elevation Profile"])
ax.margins(x=0)

The cross_section
service extracts the cross-section within a buffer distance from a point of interest along a flow line and at a given number of points.
[7]:
gdf = gpd.GeoDataFrame(
{
"numpts": [
101,
],
"width": [
1000,
],
},
geometry=[sgeom.Point((-103.80119, 40.2684))],
crs="epsg:4326",
)
section = pynhd.pygeoapi(gdf, "cross_section")
[8]:
ax = section[["distance", "elevation"]].plot(
x="distance", y="elevation", color="b", linewidth=3.0, figsize=(10, 5)
)
ax.legend(["Cross Section"])
ax.margins(x=0)
