Note
Go to the end to download the full example code.
SNOTEL and snow courses (USDA NRCS AWDB)#
The USDA NRCS Air and Water Database holds the SNOTEL and SNOLITE telemetry
sites, the SCAN soil-climate network and the manual snow courses of the
western United States and parts of western Canada. SNOTEL sites report SWE,
snow depth, precipitation and temperature every day (hourly for most);
snow courses are measured by hand around the first of the month through the
winter. Both kinds sit in the inventory, and daily_only keeps the first.
There is one route, source="awdb", the AWDB REST API v1. It needs no
account. The API answers in inches and degrees Fahrenheit; the client converts
to centimetres and degrees Celsius before anything reaches the Dataset.
The figures show the stations around a volcano in the Cascades, coloured by kind of site; one water year of SWE and snow depth at the daily sites; and five water years of SWE at one site on a single October-to-September axis.
import matplotlib.pyplot as plt
import easysnowdata as esd
aoi = (-122.1, 46.6, -121.3, 47.1)
box = esd.parse_aoi(aoi)
print(box, box.utm_crs) # the map below is drawn in the AOI's UTM zone
AOI(bounds=(-122.1000, 46.6000, -121.3000, 47.1000), clip=True) EPSG:32610
Where the product comes from.
for src in esd.catalog.get("awdb-stations").sources:
print(f"{src.id:20} {src.title:45} {', '.join(src.requires) or 'no account'}")
awdb USDA NRCS AWDB REST API v1 no account
Every AWDB station in the box, then only those whose daily record the
archive pipeline has actually retrieved. The network_code column is
AWDB’s own kind-of-site code: SNTL for SNOTEL, SNOW for a course.
everything_gdf = esd.stations.inventory(aoi, networks="awdb")
daily_gdf = esd.stations.inventory(aoi, networks="awdb", daily_only=True)
print(f"{len(everything_gdf)} AWDB stations, {len(daily_gdf)} with a daily record")
print(everything_gdf["network_code"].value_counts().to_string())
24 AWDB stations, 11 with a daily record
network_code
SNOW 13
SNTL 11
The snow courses sit at the same passes as the SNOTEL sites, often a few hundred metres apart: many of them predate the telemetry and were kept.
ax = esd.plotting.points(
everything_gdf.to_crs(box.utm_crs),
column="network_code",
legend_label="AWDB site kind",
title="AWDB stations around Mount Rainier",
)
One water year of SWE and snow depth at the daily sites, in centimetres.
obs_ds = esd.stations.load(daily_gdf, variables=["swe", "snwd"], time="2023-10/2024-09")
print(obs_ds)
<xarray.Dataset> Size: 85kB
Dimensions: (station: 11, time: 366)
Coordinates: (12/30)
* station (station) <U12 528B '375_WA_SNTL' ... '863_WA_SNTL'
network (station) <U4 176B 'awdb' 'awdb' ... 'awdb' 'awdb'
name (station) <U17 748B 'Bumping Ridge' ... 'White Pass...
client (station) <U4 176B 'awdb' 'awdb' ... 'awdb' 'awdb'
network_code (station) <U4 176B 'SNTL' 'SNTL' ... 'SNTL' 'SNTL'
latitude (station) float64 88B 46.81 47.04 ... 46.64 46.64
... ...
metadata_fetched_at (station) <U19 836B '2026-11-07 00:00:00' ... '2026...
station_id (station) <U12 528B '375:WA:SNTL' ... '863:WA:SNTL'
hemisphere (station) <U8 352B 'northern' ... 'northern'
* time (time) datetime64[us] 3kB 2023-10-01 ... 2024-09-30
water_year (time) int64 3kB 2024 2024 2024 ... 2024 2024 2024
dowy (time) int64 3kB 1 2 3 4 5 6 ... 362 363 364 365 366
Data variables:
snwd (station, time) float64 32kB 0.0 0.0 0.0 ... 0.0 0.0
swe (station, time) float64 32kB 0.0 0.0 0.0 ... 0.0 0.0
Attributes:
source: awdb
source_id: awdb
source_title: USDA NRCS AWDB REST API v1
source_url: https://wcc.sc.egov.usda.gov/awdbRestApi/services/v1
product_id: awdb-stations
title: SNOTEL, SCAN, and snow-course observations (USDA N...
data_citation: USDA Natural Resources Conservation Service, Natio...
license: Public domain (US government data)
easysnowdata_version: 0.3.9.dev1+gd97689c2f
interval: daily
The legend labels come from the name and elevation_m coordinates,
so the series read as stations rather than as codes.
fig, axes = plt.subplots(2, 1, figsize=(9, 7), sharex=True)
esd.plotting.timeseries(obs_ds["swe"], ax=axes[0], title="Water year 2024")
esd.plotting.timeseries(obs_ds["snwd"], ax=axes[1], title="", legend=False)
fig.tight_layout()
Snow depth peaks weeks before SWE does: the pack keeps taking on water as it settles and densifies, so the deepest day and the day of peak SWE are not the same day.
peak_snwd_da = obs_ds["snwd"].idxmax(dim="time").dt.strftime("%Y-%m-%d")
peak_swe_da = obs_ds["swe"].idxmax(dim="time").dt.strftime("%Y-%m-%d")
for station in obs_ds["station"].values:
print(
f"{str(obs_ds['name'].sel(station=station).values):16}"
f" deepest {str(peak_snwd_da.sel(station=station).values)}"
f" peak SWE {str(peak_swe_da.sel(station=station).values)}"
)
Bumping Ridge deepest 2024-03-13 peak SWE 2024-03-30
Burnt Mountain deepest 2024-03-05 peak SWE 2024-04-09
Cayuse Pass deepest 2024-03-12 peak SWE 2024-04-10
Corral Pass deepest 2024-03-13 peak SWE 2024-04-11
Huckleberry Creek deepest 2024-01-12 peak SWE 2024-01-22
Morse Lake deepest 2024-03-12 peak SWE 2024-04-03
Mowich deepest 2024-01-11 peak SWE 2024-01-11
Paradise deepest 2024-03-13 peak SWE 2024-05-10
Pigtail Peak deepest 2024-03-13 peak SWE 2024-05-11
Skate Creek deepest 2024-03-13 peak SWE 2024-03-15
White Pass E.S. deepest 2024-03-13 peak SWE 2024-03-30
Five water years at Paradise, overlaid on one October-to-September axis.
paradise_ds = esd.stations.load("679_WA_SNTL", variables="swe", time="2019-10/2024-09")
ax = esd.plotting.timeseries(
paradise_ds["swe"],
by_water_year=True,
title=f"{str(paradise_ds['name'].values[0])}"
f" ({float(paradise_ds['elevation_m'].values[0]):.0f} m): five water years",
)
