How much of my region’s water supply is locked in snow, and is the snowpack shrinking?

Snow is a savings account for water. In many places the rain and snow of winter pile up as snowpack in the mountains, then slowly melt through spring and summer to feed rivers and reservoirs. Satellites can’t weigh that snow directly, but they can see where snow is sitting from orbit. So you can watch the white patch grow each winter, shrink each summer, and check whether it’s getting smaller year over year.

Which data to use, and how

Use this dataset: MODIS MOD10A1 (short name MOD10A1), daily snow-cover maps at 500 m going back to 2000. It’s the workhorse for this question and the easiest to start with. Later you can extend the record with VIIRS VNP10A1 (VNP10A1, 2012+), the newer sensor that takes over as the old satellite ages out.

A quick word on what these see. The satellite computes an NDSI (Normalized Difference Snow Index), a brightness ratio that lights up where snow is. It tells you snow is there, not how deep it is.

Then do four steps. The runnable code further down does all of this; this is just the idea.

1. Pick your basin (a latitude/longitude box around your mountains) and a water year (the hydrologist’s calendar: Oct 1 to Sep 30, so one full snow season stays in one bucket).
2. Count the snow each day. For every daily map, count the pixels flagged as snow (NDSI band value ≥ 40, the standard snow/no-snow cut) and turn that into snow-covered area (SCA, how many km² are white).
3. Draw the season’s curve. Plot SCA day by day. It rises through winter, peaks (around April 1 in the Sierra Nevada), then melts away. The date it hits zero is the melt-out date.
4. Compare to normal and check the trend. Overlay this year on the 2000-present average, and fit a line through all the years to see if the snow season is getting shorter (fewer snow-covered days per decade). That’s the fingerprint of a shrinking snowpack.

That’s the whole method. Below, “How a scientist answers this” names the exact tools, and the Run it cell does it live on synthetic data so you can see each step work.

What you can find out

• How much of your region is under snow on any day, daily, at 500 m.
• The accumulation-and-melt curve: when the snow peaks (typically ~April 1 in the Sierra) and when it finally disappears.
• Whether the snow season is shrinking. Is the snow-on/snow-off window getting shorter across the 2000-present record?
• Where the snow line is creeping. Combine the snow maps with a terrain map (a DEM, a Digital Elevation Model = a grid of ground heights) to find the lowest elevation that stays snow-covered each year.
• This year vs a normal year: this winter’s peak coverage and melt-out date against the long-term average.

What it can’t tell you

• How much actual water is in the snow. Optical satellites see extent (the white patch), not depth or water content. The water amount, called snow water equivalent (SWE) (the depth of water you’d get if you melted the snow flat), needs a different tool: passive microwave (AMSR2 / SMAP, but coarse ~10–25 km), a model like Daymet (North America only), airborne lidar surveys, or ground stations (SNOTEL).
• Snow hidden under thick forest. Tree canopy blinds the optical view; it can’t see the snow on the ground beneath. Microwave and modeled products try to fill this gap, but with large uncertainties of their own.
• Anything on a cloudy day. These satellites see nothing through cloud. You have to gap-fill (merge the two satellites, or stitch nearby days together) and honestly flag the days you guessed.
• The actual river flow you’ll get. Turning snow into streamflow needs a full hydrology model (accounting for evaporation, soil moisture, and how water routes downhill). Snow data is an input to that, not the final answer.

Gotchas to watch for

• The snow/no-snow line is a choice. “Snow” here means NDSI ≥ 0.4 (band value ≥ 40), the standard cut. A different threshold gives a different snow area, so always write down the one you used so your numbers are reproducible.
• Merging the two satellites fills cloud holes, but they fly at different times. Pairing MODIS MYD10A1 (MYD10A1, the Aqua satellite) with MOD10A1 (the Terra satellite) plugs gaps left by clouds, but the two pass overhead ~3 hours apart. That time gap matters if you’re measuring how fast snow melts within a day, so note it.
• Don’t mix data versions in one trend. The current data release is Collection 6.1 (MOD10A1.061). Older Collection 6 values differ slightly, so blending C6 and C6.1 in a single long trend can fake a change that’s really just a version switch.
• Microwave SWE lies in deep or wet snow. The AMSR2 / SMAP microwave estimate saturates above ~150–200 mm of water and breaks down in melting snow, so don’t trust it as truth in a deep mountain snowpack.
• Daymet is a model, not a measurement. Its SWE is a reanalysis-style estimate (a computer’s best reconstruction), and it only covers North America. Label it as model-based, never as observed.
• 25 years is short for water. The satellite record starts in 2000. That’s a real trend, but a short one for hydrology, so say so when you report a slope.

The real-data code

The Run it cell above runs this method on synthetic data with no login. Below is the whole recipe against the real archive: search MOD10A1, open each daily granule, count the snow pixels into a snow-covered-area curve, compare this year to the 2000-present normal, fit the duration trend, and plot it. It needs a free Earthdata Login.

import earthaccess, numpy as np, rioxarray
from scipy import stats
import matplotlib.pyplot as plt

earthaccess.login(strategy="netrc")

aoi = (-120.6, 36.4, -118.2, 39.3)   # Sierra Nevada (CA snow-fed basins) as (W, S, E, N)
years = range(2001, 2026)            # water years (Oct of year-1 -> Sep of year)

# 1. One snow-covered-area number per day: open each MOD10A1 granule, keep the
#    NDSI_Snow_Cover band, count pixels >= 40 (the standard snow cut), convert
#    500 m pixels to km^2. Flags (200-255: cloud, water, fill) fall out of the >=40 cut.
PIX_KM2 = 0.463 * 0.463   # MODIS sinusoidal 500 m pixel area
def daily_sca(y):
    grans = earthaccess.search_data(
        short_name="MOD10A1", version="61", bounding_box=aoi,
        temporal=(f"{y-1}-10-01", f"{y}-09-30"))
    out = {}
    for g in earthaccess.open(grans):
        da = rioxarray.open_rasterio(g, variable="NDSI_Snow_Cover").squeeze()
        sca = float(((da >= 40) & (da <= 100)).sum()) * PIX_KM2
        doy = int(g.granule["umm"]["TemporalExtent"]["RangeDateTime"]
                  ["BeginningDateTime"][:10].split("-")[2])  # day within month is enough to order
        out[g.granule["meta"]["native-id"]] = sca
    return np.array(list(out.values()))

# Build a year x day-of-water-year matrix of SCA (NaN-pad short years to align days).
rows = [daily_sca(y) for y in years]
W = max(len(r) for r in rows)
sca = np.full((len(rows), W), np.nan)
for i, r in enumerate(rows):
    sca[i, :len(r)] = r

# 2. The season curve + the long-term normal (per day-of-water-year mean).
this_year = sca[-1]
normal = np.nanmean(sca, axis=0)

# 3. Snow-cover duration per year (days with SCA above a small basin-fraction
#    floor), then a least-squares trend in days/decade across 2000-present.
basin_km2 = (aoi[2]-aoi[0]) * (aoi[3]-aoi[1]) * 111 * 111
duration = np.nansum(sca > 0.05 * basin_km2, axis=1).astype(float)
yr = np.array(list(years))
slope, intercept, r, p, se = stats.linregress(yr, duration)
print(f"Snow-cover duration trend: {slope*10:+.1f} days/decade, "
      f"p={p:.3f}", "(significant)" if p < 0.05 else "(not distinguishable from noise)")

# 4. Plot this year's SCA curve against the normal, plus the duration trend.
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].plot(normal, "k--", label="2000-present normal"); ax[0].plot(this_year, label=f"WY{years[-1]}")
ax[0].set_xlabel("day of water year"); ax[0].set_ylabel("snow-covered area (km²)"); ax[0].legend()
ax[1].plot(yr, duration, "o"); ax[1].plot(yr, intercept + slope*yr, "r-")
ax[1].set_xlabel("water year"); ax[1].set_ylabel("snow-cover days"); plt.tight_layout(); plt.show()

# Before you trust it: cross-check the trend against SNOTEL ground stations (or modeled
# Daymet SWE), and remember NDSI is snow *extent*, not depth — and 25 years is short for water.

Where the data comes from

Everything comes through one free login (Earthdata Login). The snow-cover maps (MOD10A1, MYD10A1, VNP10A1) live at NASA’s snow-and-ice archive (NSIDC); the modeled Daymet SWE comes from a separate NASA archive (ORNL). The earthaccess library reaches both with the same single login, so you don’t juggle accounts.

Sources

• MODIS MOD10A1 (Terra daily snow cover): https://nsidc.org/data/mod10a1
• MODIS MYD10A1 (Aqua daily snow cover): https://nsidc.org/data/myd10a1
• VIIRS VNP10A1 (daily snow cover): https://nsidc.org/data/vnp10a1
• Daymet V4 (includes SWE): https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=2129
• NSIDC snow & ice data overview: https://nsidc.org/learn/parts-cryosphere/snow