Note
Go to the end to download the full example code.
Basins: HUC, HydroBASINS and GRDC#
Four basin products answer the same question, “which drainage is this in?”, at different scales and for different parts of the world. The USGS Watershed Boundary Dataset (HUC) nests the United States from 2-digit regions to 16-digit sub-watersheds. HydroBASINS is the global equivalent, Pfafstetter levels 1 to 12. The GRDC layers are coarse and whole-world: 405 major river basins, and the WMO basins and sub-basins with their region numbers.
Every loader returns a GeoDataFrame in EPSG:4326 holding the features that
intersect the AOI, whole and uncut. HUC defaults to the public USGS ArcGIS
REST service and needs no account; source="gee" reads the 2017 snapshot on
Earth Engine instead. HydroBASINS defaults to the 2.7 GB BasinATLAS
geodatabase (with its several hundred attributes, read with mask pushdown),
or, as here, source="hydrosheds" fetches one ~300 MB regional zip with the
geometry and Pfafstetter fields only. The GRDC layers are single zips. Nothing
executed on this page needs credentials.
The first figure fills the HUC12 sub-watersheds inside the HUC8 subbasins, then lays the HUC8 divides over HydroBASINS level 8: two independent delineations of the same drainage. The second nests level 6 inside the WMO basins and steps out to the GRDC major basin the whole AOI drains to.
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
import easysnowdata as esd
aoi = (-121.94, 46.72, -121.54, 46.99) # Mount Rainier, WA
box = esd.parse_aoi(aoi)
utm = box.utm_crs # the local maps are drawn in the AOI's UTM zone
print(box, utm)
AOI(bounds=(-121.9400, 46.7200, -121.5400, 46.9900), clip=True) EPSG:32610
Where each product comes from, and what the route asks for.
for product_id in ("huc", "hydrobasins", "grdc-major-river-basins", "grdc-wmo-basins"):
for src in esd.catalog.get(product_id).sources:
print(
f"{product_id:24} {src.id:20} {src.title:28} {', '.join(src.requires) or 'no account'}"
)
huc usgs-wbd USGS WBD (ArcGIS REST) no account
huc gee Earth Engine earthengine
hydrobasins figshare-basinatlas figshare BasinATLAS gdb no account
hydrobasins hydrosheds HydroSHEDS regional zips no account
hydrobasins gee Earth Engine earthengine
grdc-major-river-basins world-bank World Bank zip no account
grdc-wmo-basins grdc GRDC zip no account
HUC8 subbasins and HUC12 sub-watersheds from the REST service. The loader pages the service in small blocks, because one unpaged HUC12 query over a large box times out and comes back as HTML rather than GeoJSON.
huc8_gdf = esd.hydro.basins.huc(aoi, level=8)
huc12_gdf = esd.hydro.basins.huc(aoi, level=12)
print(huc8_gdf[["name", "huc8", "areasqkm", "states"]].to_string(index=False))
print(f"{len(huc12_gdf)} HUC12 sub-watersheds intersect the box")
name huc8 areasqkm states
Nisqually 17110015 1994.47 WA
Puyallup 17110014 2533.97 WA
Upper Cowlitz 17080004 2655.30 WA
20 HUC12 sub-watersheds intersect the box
The same call with source="gee" returns the same columns from the 2017
WBD snapshot on Earth Engine, for users who already hold an account:
huc8_gee_gdf = esd.hydro.basins.huc(aoi, level=8, source="gee")
The geometries differ slightly from the REST service’s, which is current and rounds coordinates to six decimals.
HydroBASINS, here from the HydroATLAS copy on Earth Engine. The HydroSHEDS
regional zip (source="hydrosheds", no account, cached after the first
download) serves the same polygons, but that server refuses cloud-runner
addresses such as the one this page is built on, so the Earth Engine route
is used for the build; on your own machine either works. PFAF_ID is the
Pfafstetter code: every digit is one level of nesting, so a level-8 basin’s
level-6 parent is its first six digits. Level 8 here is 1 500 to 2 600 km²
per basin, the size of a HUC8.
hybas_source = "hydrosheds"
try:
hybas8_gdf = esd.hydro.basins.hydrobasins(aoi, level=8, source=hybas_source)
except Exception as exc: # the HydroSHEDS server blocks some networks
print(
f"HydroSHEDS route unavailable here ({type(exc).__name__}); using Earth Engine"
)
hybas8_gdf = esd.hydro.basins.hydrobasins(aoi, level=8, source="gee")
print(hybas8_gdf[["HYBAS_ID", "PFAF_ID", "SUB_AREA", "UP_AREA"]].to_string(index=False))
0%| | 0.00/317M [00:00<?, ?B/s]
1%|▎ | 2.59M/317M [00:00<00:12, 25.9MB/s]
4%|█▌ | 12.7M/317M [00:00<00:04, 70.2MB/s]
7%|██▊ | 23.3M/317M [00:00<00:03, 86.5MB/s]
11%|████ | 33.8M/317M [00:00<00:03, 93.8MB/s]
14%|█████▎ | 44.4M/317M [00:00<00:02, 98.3MB/s]
17%|██████▊ | 55.0M/317M [00:00<00:02, 101MB/s]
21%|████████ | 65.5M/317M [00:00<00:02, 102MB/s]
24%|█████████▍ | 76.1M/317M [00:00<00:02, 104MB/s]
27%|██████████▋ | 86.7M/317M [00:00<00:02, 104MB/s]
31%|███████████▉ | 97.3M/317M [00:01<00:02, 105MB/s]
34%|█████████████▋ | 108M/317M [00:01<00:01, 105MB/s]
37%|██████████████▉ | 118M/317M [00:01<00:01, 105MB/s]
41%|████████████████▎ | 129M/317M [00:01<00:01, 105MB/s]
44%|█████████████████▋ | 140M/317M [00:01<00:01, 106MB/s]
47%|██████████████████▉ | 150M/317M [00:01<00:01, 106MB/s]
51%|████████████████████▎ | 161M/317M [00:01<00:01, 106MB/s]
54%|█████████████████████▋ | 171M/317M [00:01<00:01, 106MB/s]
58%|███████████████████████ | 182M/317M [00:01<00:01, 106MB/s]
61%|████████████████████████▎ | 193M/317M [00:01<00:01, 106MB/s]
64%|█████████████████████████▋ | 203M/317M [00:02<00:01, 106MB/s]
68%|███████████████████████████ | 214M/317M [00:02<00:00, 106MB/s]
71%|████████████████████████████▎ | 225M/317M [00:02<00:00, 106MB/s]
74%|█████████████████████████████▋ | 235M/317M [00:02<00:00, 106MB/s]
78%|███████████████████████████████ | 246M/317M [00:02<00:00, 106MB/s]
81%|████████████████████████████████▍ | 256M/317M [00:02<00:00, 106MB/s]
84%|█████████████████████████████████▋ | 267M/317M [00:02<00:00, 106MB/s]
88%|███████████████████████████████████ | 278M/317M [00:02<00:00, 106MB/s]
91%|████████████████████████████████████▍ | 288M/317M [00:02<00:00, 106MB/s]
94%|█████████████████████████████████████▊ | 299M/317M [00:02<00:00, 106MB/s]
98%|███████████████████████████████████████ | 310M/317M [00:03<00:00, 106MB/s]
0%| | 0.00/317M [00:00<?, ?B/s]
100%|███████████████████████████████████████| 317M/317M [00:00<00:00, 1.56TB/s]
HYBAS_ID PFAF_ID SUB_AREA UP_AREA
7080371080 78210029 1520.0 1520.3
7080015820 78310380 1878.7 1878.7
7080015900 78310400 2637.9 2637.9
Left: HUC12s coloured by name inside the HUC8 outlines. Right: the HUC8
divides over HydroBASINS level 8. The two products were drawn from
different elevation models by different rules, and their main divides still
agree closely; the units are not one-to-one, though: 78210029 covers
1 520 km² of the 2 655 km² Upper Cowlitz HUC8, the rest falling in level-8
units that do not touch the box and so were not returned. Vector layers
draw with geopandas, here reprojected to the AOI’s UTM zone; finish_map
adds the basemap, graticule and scale bar.
huc8_utm_gdf, huc12_utm_gdf, hybas8_utm_gdf = (
basin_gdf.to_crs(utm) for basin_gdf in (huc8_gdf, huc12_gdf, hybas8_gdf)
)
fig, axes = plt.subplots(1, 2, figsize=(12, 5.6))
huc12_utm_gdf.plot(
ax=axes[0],
column="name",
cmap="tab20",
alpha=0.55,
edgecolor="white",
linewidth=0.5,
)
huc8_utm_gdf.boundary.plot(ax=axes[0], color="black", linewidth=1.6)
for _, row in huc8_utm_gdf.iterrows():
axes[0].annotate(
row["name"],
row.geometry.centroid.coords[0],
ha="center",
fontsize=9,
fontweight="bold",
)
axes[0].set_title("HUC12 sub-watersheds inside the HUC8 subbasins (USGS WBD)")
palette8 = dict(
zip(hybas8_gdf["PFAF_ID"], ("#fdc086", "#beaed4", "#7fc97f"), strict=True)
)
hybas8_utm_gdf.plot(
ax=axes[1],
color=hybas8_utm_gdf["PFAF_ID"].map(palette8),
alpha=0.55,
edgecolor="white",
linewidth=0.8,
)
huc8_utm_gdf.boundary.plot(ax=axes[1], color="black", linewidth=1.4, linestyle="--")
axes[1].legend(
handles=[
*[
Patch(facecolor=c, alpha=0.55, label=f"HydroBASINS {k}")
for k, c in palette8.items()
],
Line2D(
[],
[],
color="black",
linewidth=1.4,
linestyle="--",
label="HUC8 (USGS WBD)",
),
],
fontsize=8,
loc="lower right",
)
axes[1].set_title("HUC8 divides over HydroBASINS level 8")
x0, y0, x1, y1 = box.to_crs(utm).total_bounds
for ax in axes:
ax.set_xlim(x0 - 30_000, x1 + 30_000)
ax.set_ylim(y0 - 30_000, y1 + 30_000)
esd.plotting.finish_map(ax, utm, basemap=True)
fig.tight_layout()
The GRDC layers. The WMO basins carry region numbers and names; the major
river basins are one polygon each. The shipped area columns of the major
basins layer (ACRES, LAEA_HA) are zero in this release, so the area
below is taken from the geometry in an equal-area projection.
wmo_gdf = esd.hydro.basins.grdc_wmo(aoi)
print(
wmo_gdf[["WMOBB", "WMOBB_NAME", "REGNAME", "SUM_SUB_AREA"]].to_string(index=False)
)
major_gdf = esd.hydro.basins.grdc_major(aoi)
area_km2 = major_gdf.to_crs("ESRI:102008").area.iloc[0] / 1e6
print(
f"{major_gdf.iloc[0]['NAME']} basin: shipped ACRES={major_gdf.iloc[0]['ACRES']}, geometry {area_km2:,.0f} km2"
)
0%| | 0.00/379M [00:00<?, ?B/s]
0%| | 8.19k/379M [00:00<1:24:42, 74.5kB/s]
0%| | 72.7k/379M [00:00<18:02, 350kB/s]
0%| | 169k/379M [00:00<11:28, 550kB/s]
0%| | 378k/379M [00:00<06:16, 1.01MB/s]
0%| | 769k/379M [00:00<03:30, 1.79MB/s]
0%| | 1.17M/379M [00:00<02:45, 2.29MB/s]
0%|▏ | 1.63M/379M [00:00<02:17, 2.75MB/s]
1%|▎ | 3.21M/379M [00:00<01:02, 5.96MB/s]
1%|▍ | 4.67M/379M [00:01<00:47, 7.80MB/s]
2%|▋ | 7.27M/379M [00:01<00:31, 11.9MB/s]
2%|▉ | 8.77M/379M [00:01<00:31, 11.9MB/s]
3%|█▎ | 12.8M/379M [00:01<00:20, 17.6MB/s]
4%|█▌ | 15.4M/379M [00:01<00:19, 18.5MB/s]
5%|█▋ | 17.2M/379M [00:01<00:21, 17.2MB/s]
5%|█▉ | 19.5M/379M [00:01<00:20, 17.8MB/s]
6%|██▏ | 21.3M/379M [00:01<00:21, 16.8MB/s]
6%|██▍ | 24.3M/379M [00:02<00:18, 18.9MB/s]
7%|██▌ | 26.1M/379M [00:02<00:19, 17.7MB/s]
8%|███ | 30.1M/379M [00:02<00:15, 21.9MB/s]
9%|███▍ | 34.2M/379M [00:02<00:13, 25.1MB/s]
10%|███▊ | 38.2M/379M [00:02<00:12, 27.2MB/s]
11%|████▏ | 42.3M/379M [00:02<00:11, 28.7MB/s]
12%|████▋ | 46.3M/379M [00:02<00:11, 29.5MB/s]
13%|█████ | 50.3M/379M [00:02<00:10, 30.2MB/s]
14%|█████▍ | 54.4M/379M [00:03<00:10, 30.7MB/s]
15%|█████▊ | 58.4M/379M [00:03<00:10, 31.1MB/s]
16%|██████▎ | 62.5M/379M [00:03<00:10, 31.4MB/s]
18%|██████▋ | 66.5M/379M [00:03<00:09, 31.7MB/s]
19%|███████ | 70.6M/379M [00:03<00:09, 31.8MB/s]
20%|███████▍ | 74.6M/379M [00:03<00:09, 31.9MB/s]
21%|███████▉ | 78.7M/379M [00:03<00:09, 32.0MB/s]
22%|████████▎ | 82.7M/379M [00:03<00:09, 32.0MB/s]
23%|████████▋ | 86.8M/379M [00:04<00:09, 32.1MB/s]
24%|█████████ | 90.9M/379M [00:04<00:08, 32.1MB/s]
25%|█████████▌ | 95.0M/379M [00:04<00:08, 32.2MB/s]
26%|█████████▉ | 99.0M/379M [00:04<00:08, 32.1MB/s]
27%|██████████▌ | 103M/379M [00:04<00:08, 32.0MB/s]
28%|███████████ | 107M/379M [00:04<00:08, 32.0MB/s]
29%|███████████▍ | 111M/379M [00:04<00:08, 32.0MB/s]
30%|███████████▊ | 115M/379M [00:05<00:08, 31.9MB/s]
31%|████████████▎ | 119M/379M [00:05<00:08, 31.9MB/s]
33%|████████████▋ | 123M/379M [00:05<00:08, 31.9MB/s]
34%|█████████████ | 127M/379M [00:05<00:07, 31.9MB/s]
35%|█████████████▌ | 131M/379M [00:05<00:07, 32.0MB/s]
36%|█████████████▉ | 135M/379M [00:05<00:07, 32.0MB/s]
37%|██████████████▎ | 139M/379M [00:05<00:07, 32.0MB/s]
38%|██████████████▊ | 143M/379M [00:05<00:07, 32.1MB/s]
39%|███████████████▏ | 148M/379M [00:06<00:07, 32.1MB/s]
40%|███████████████▌ | 152M/379M [00:06<00:07, 32.1MB/s]
41%|████████████████ | 156M/379M [00:06<00:06, 32.1MB/s]
42%|████████████████▍ | 160M/379M [00:06<00:06, 32.0MB/s]
43%|████████████████▊ | 164M/379M [00:06<00:06, 32.0MB/s]
44%|█████████████████▎ | 168M/379M [00:06<00:06, 32.0MB/s]
45%|█████████████████▋ | 172M/379M [00:06<00:06, 32.1MB/s]
46%|██████████████████ | 176M/379M [00:06<00:06, 32.1MB/s]
48%|██████████████████▌ | 180M/379M [00:07<00:06, 32.0MB/s]
49%|██████████████████▉ | 184M/379M [00:07<00:06, 32.1MB/s]
50%|███████████████████▎ | 188M/379M [00:07<00:05, 32.0MB/s]
51%|███████████████████▊ | 192M/379M [00:07<00:05, 31.8MB/s]
52%|████████████████████▏ | 196M/379M [00:07<00:05, 31.7MB/s]
53%|████████████████████▌ | 200M/379M [00:07<00:05, 31.4MB/s]
54%|████████████████████▉ | 203M/379M [00:07<00:06, 29.2MB/s]
55%|█████████████████████▎ | 207M/379M [00:07<00:05, 29.7MB/s]
56%|█████████████████████▋ | 211M/379M [00:08<00:05, 30.5MB/s]
57%|██████████████████████ | 214M/379M [00:08<00:05, 29.1MB/s]
58%|██████████████████████▍ | 218M/379M [00:08<00:05, 30.1MB/s]
59%|██████████████████████▉ | 222M/379M [00:08<00:05, 30.7MB/s]
60%|███████████████████████▎ | 226M/379M [00:08<00:05, 29.9MB/s]
61%|███████████████████████▋ | 230M/379M [00:08<00:04, 30.4MB/s]
62%|████████████████████████ | 234M/379M [00:08<00:04, 31.0MB/s]
63%|████████████████████████▌ | 238M/379M [00:08<00:04, 31.3MB/s]
64%|████████████████████████▉ | 242M/379M [00:09<00:04, 31.6MB/s]
65%|█████████████████████████▎ | 246M/379M [00:09<00:04, 31.7MB/s]
66%|█████████████████████████▊ | 250M/379M [00:09<00:04, 31.8MB/s]
67%|██████████████████████████▏ | 254M/379M [00:09<00:04, 31.0MB/s]
68%|██████████████████████████▌ | 258M/379M [00:09<00:03, 31.2MB/s]
69%|██████████████████████████▉ | 262M/379M [00:09<00:03, 31.5MB/s]
70%|███████████████████████████▍ | 266M/379M [00:09<00:03, 31.7MB/s]
71%|███████████████████████████▊ | 270M/379M [00:09<00:03, 29.4MB/s]
72%|████████████████████████████▏ | 274M/379M [00:10<00:03, 30.2MB/s]
73%|████████████████████████████▋ | 278M/379M [00:10<00:03, 30.7MB/s]
75%|█████████████████████████████ | 282M/379M [00:10<00:03, 31.2MB/s]
76%|█████████████████████████████▍ | 286M/379M [00:10<00:02, 31.4MB/s]
77%|█████████████████████████████▉ | 291M/379M [00:10<00:02, 31.7MB/s]
78%|██████████████████████████████▎ | 295M/379M [00:10<00:02, 31.8MB/s]
79%|██████████████████████████████▊ | 299M/379M [00:10<00:02, 31.9MB/s]
80%|███████████████████████████████▏ | 303M/379M [00:10<00:02, 31.8MB/s]
81%|███████████████████████████████▌ | 307M/379M [00:11<00:02, 31.9MB/s]
82%|███████████████████████████████▉ | 311M/379M [00:11<00:02, 31.8MB/s]
83%|████████████████████████████████▍ | 315M/379M [00:11<00:02, 31.9MB/s]
84%|████████████████████████████████▊ | 319M/379M [00:11<00:01, 32.0MB/s]
85%|█████████████████████████████████▏ | 323M/379M [00:11<00:01, 31.9MB/s]
86%|█████████████████████████████████▋ | 327M/379M [00:11<00:01, 32.0MB/s]
87%|██████████████████████████████████ | 331M/379M [00:11<00:01, 32.0MB/s]
88%|██████████████████████████████████▍ | 335M/379M [00:12<00:01, 32.0MB/s]
90%|██████████████████████████████████▉ | 339M/379M [00:12<00:01, 32.0MB/s]
91%|███████████████████████████████████▎ | 343M/379M [00:12<00:01, 32.1MB/s]
92%|███████████████████████████████████▊ | 347M/379M [00:12<00:00, 32.1MB/s]
93%|████████████████████████████████████▏ | 351M/379M [00:12<00:00, 32.1MB/s]
94%|████████████████████████████████████▌ | 355M/379M [00:12<00:00, 32.1MB/s]
95%|█████████████████████████████████████ | 359M/379M [00:12<00:00, 32.0MB/s]
96%|█████████████████████████████████████▍ | 363M/379M [00:12<00:00, 32.0MB/s]
97%|█████████████████████████████████████▊ | 367M/379M [00:13<00:00, 31.9MB/s]
98%|██████████████████████████████████████▏| 371M/379M [00:13<00:00, 31.9MB/s]
99%|██████████████████████████████████████▋| 375M/379M [00:13<00:00, 31.9MB/s]
0%| | 0.00/379M [00:00<?, ?B/s]
100%|███████████████████████████████████████| 379M/379M [00:00<00:00, 1.74TB/s]
WMOBB WMOBB_NAME REGNAME SUM_SUB_AREA
415 COLUMBIA (ex 416) North America, Central America and the Caribbean 390679.1
446 North Pacific (ex 341, 405, 406, 445, 455, 456, 457, 471, 474, 477) North America, Central America and the Caribbean 386726.3
Columbia basin: shipped ACRES=0.002, geometry 670,165 km2
Left: level-6 HydroBASINS (fill) inside the WMO basins (outline: the Columbia, and the North Pacific coastal basins), with the AOI in red. Right: the GRDC major basin, the Columbia, with the same AOI as a dot, to show the scale jump from a sub-watershed to a continental basin.
hybas6_gdf = esd.hydro.basins.hydrobasins(aoi, level=6, source=hybas_source)
fig, axes = plt.subplots(1, 2, figsize=(12, 5.6))
palette6 = dict(
zip(hybas6_gdf["PFAF_ID"], ("#66c2a5", "#e6f598", "#8da0cb"), strict=True)
)
hybas6_gdf.plot(
ax=axes[0],
color=hybas6_gdf["PFAF_ID"].map(palette6),
alpha=0.6,
edgecolor="white",
linewidth=0.8,
)
wmo_gdf.boundary.plot(ax=axes[0], color="black", linewidth=1.4)
esd.aoi.parse_aoi(aoi).footprint.boundary.plot(ax=axes[0], color="red", linewidth=1.5)
handles = [
Patch(facecolor=c, alpha=0.6, label=f"PFAF {k}") for k, c in palette6.items()
]
handles += [
Line2D([], [], color="black", linewidth=1.4, label="GRDC / WMO basins"),
Line2D([], [], color="red", linewidth=1.5, label="AOI"),
]
axes[0].legend(handles=handles, fontsize=8, loc="upper right")
x0, y0, x1, y1 = hybas6_gdf.total_bounds
axes[0].set_xlim(x0 - 0.2, x1 + 0.2)
axes[0].set_ylim(y0 - 0.2, y1 + 0.2)
axes[0].set_title("HydroBASINS level 6 inside the GRDC / WMO basins")
major_gdf.plot(
ax=axes[1], facecolor="#cfe3f7", edgecolor="#1f4e79", linewidth=1.2, alpha=0.8
)
axes[1].plot(
*esd.aoi.parse_aoi(aoi).geometry.centroid.coords[0],
marker="o",
color="red",
markersize=6,
)
axes[1].set_title(f"GRDC major basin: {major_gdf.iloc[0]['NAME']}, {area_km2:,.0f} km²")
for ax in axes:
esd.plotting.finish_map(ax, hybas6_gdf.crs, basemap=True)
fig.tight_layout()
Total running time of the script: (0 minutes 36.739 seconds)