Is my river’s snowmelt flood coming earlier as the mountains warm?

A spring flood on a snow-fed river is the mountains emptying out. Snow piles up all winter, then melts, runs downhill, and the river rises. Snow cover is just how much of the ground is white, and a satellite can see it from orbit, so you can watch the white retreat off your mountains every spring. As the climate warms, that melt tends to start sooner. So the real question is whether the snow over your headwaters is vanishing earlier than it used to.

Which data to use, and how

Use this dataset: MOD10A1 (short name MOD10A1), a daily, 500 m map of snow cover from NASA’s MODIS sensor that reaches back to 2000. It’s the one dataset you need for the core question. Everything else here is optional context.

A note on the number it reports. Each pixel carries an NDSI value (Normalized Difference Snow Index, a 0–100 score for how snowy that 500 m patch looks). 0–100 is real snow data. Anything above 100 is a flag for cloud, water, or missing data, which you throw away.

Verified locally. For the Colorado River headwaters (39.0–39.8 °N), the MOD10A1 NDSI_Snow_Cover field read 42.1% of the area still snow-covered (NDSI 40–100) on 1 April 2024. That’s a single daily snapshot of how full the snowpack still was that morning. Track that fraction day by day across a spring, and the date it crosses, say, 50% gives you a clean “melt-out timing” number you can compare year to year.

Then do four steps (the runnable code further down does all of this; this is just the idea):

1. Pick your basin, a small latitude/longitude box around your headwaters, and the spring you care about.
2. Snow vs no-snow each day. For one day’s map, count the snow-covered pixels (NDSI 40–100) and divide by the valid pixels (0–100, ignoring cloud/fill flags). That gives one number: the snow-covered fraction.
3. Walk through the spring. Repeat step 2 day by day and watch that fraction fall as the snow melts off.
4. Find the melt-out date. Note the day the fraction first drops below your threshold (e.g. 50%). That single date is your “snow disappearance date”. Line it up across 20+ springs to see if it’s drifting earlier.

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 headwaters is snow-covered on any given day, back to 2000: the fraction of pixels with NDSI_Snow_Cover between 40 and 100.
• When the snow melts out each spring: the date the snow-covered fraction drops below a threshold (e.g. 50%), the “snow disappearance date”.
• Whether melt-out is trending earlier. Line up that date across 20+ springs and look for a drift toward earlier dates as the mountains warm.
• How this spring compares to normal. Overlay this year’s melt curve on the multi-year average.
• Where rain is piling onto the melt. Add GPM IMERG rainfall, since spring rain falling on snow can speed up the flood.
• Whether the river downstream actually rose. Pair with SWOT river heights to connect the snowpack story to the channel.

What it can’t tell you

• How much water is in the snow. MOD10A1 sees snow area, not snow depth or snow water equivalent (the amount of water you’d get if you melted it). A thin late cover and a deep one can look identical from above. For water volume you need ground networks like SNOTEL or SNODAS, or airborne snow surveys.
• The actual flood peak or how fast the river is flowing. Snow cover is the driver, not the flow rate (discharge). For that you need a stream gauge (e.g. USGS), and SWOT gives water height, not a calibrated flow.
• What happens under clouds. MODIS is an optical (camera-like) sensor, so cloudy days leave gaps. Stormy melt periods are exactly when clouds are worst. You have to gap-fill or composite across days.
• Snow hidden by trees or in small patches. At 500 m a pixel is bigger than many features, so a dense forest canopy hides the snow underneath and tiny patches blur together. Melt timing in thick forest is the least reliable.
• Why it melted early. The data show that the timing shifted, not whether warm air, dust settling on the snow, or a light winter caused it.

Gotchas to watch for

• Throw away the flag values. Pixel values 0–100 are real snow percentages; values above 100 mean cloud, water, or no-data (fill). Average those in by accident and you corrupt the result, which is the single most common beginner mistake. Keep only 0 <= value <= 100.
• There’s no lat/lon grid to slice. MOD10A1 ships as a square tile (2400×2400 pixels) in the MODIS sinusoidal projection, a map grid measured in meters, not in latitude/longitude. So you can’t just ask for “my box”. You download a tile (the Colorado headwaters sit in tile h09v05) and convert your corner coordinates into row/column numbers to crop it. The code below does that math; just match the tile to your basin.
• Clouds leave holes. One cloudy day can blank out your basin and look like a fake “melt”. Smooth over several days or composite, and don’t trust the curve on heavily clouded dates.
• A “melt-out date” needs a threshold you choose. 50% is a convention, not a law of nature, and different thresholds shift the date you get. Pick one and use it for every year so the comparison is fair.

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, open each daily tile, clip it to your basin, count the snow fraction day by day across the spring, find the melt-out date, and plot it. It needs a free Earthdata Login.

Verified locally. MOD10A1, variable NDSI_Snow_Cover, is a daily HDF-EOS2/HDF4 tile in the MODIS sinusoidal grid. There is no global lat/lon index to slice, so download one granule, read the 2400×2400 array, and clip to your basin’s row/column window. Values 0–100 are percent snow; values above 100 are flags (cloud, water, fill), so keep only 0–100.

import os, re, math, tempfile, datetime, earthaccess, numpy as np
import matplotlib.pyplot as plt
from pyhdf.SD import SD, SDC

# load Earthdata creds from .env without `source` (passwords can break the shell)
for line in open(".env"):
    m = re.match(r'\s*(?:export\s+)?([A-Z0-9_]+)\s*=\s*(.*)\s*$', line)
    if m: os.environ.setdefault(m.group(1), m.group(2).strip().strip('"').strip("'"))
earthaccess.login(strategy="environment")   # free Earthdata Login

W, S, E, N = -106.9, 39.0, -105.9, 39.8     # your headwaters (Colorado River)
SPRING = ("2024-03-01", "2024-06-30")       # the melt season you care about
workdir = tempfile.mkdtemp()

# clip math for the sinusoidal tile (here h09v05): turn lon/lat corners into row/col
R, T, ncell = 6371007.181, 1111950.519667, 2400
cell = T / ncell
x0 = -20015109.354 + 9 * T          # tile h09 left edge (m)
y1 = 10007554.677 - 5 * T           # tile v05 top edge (m)
def rc(lon, lat):
    x = R * math.radians(lon) * math.cos(math.radians(lat))
    y = R * math.radians(lat)
    return (y1 - y) / cell, (x - x0) / cell
rows = [rc(lo, la)[0] for lo in (W, E) for la in (S, N)]
cols = [rc(lo, la)[1] for lo in (W, E) for la in (S, N)]
r0, r1 = max(0, int(min(rows))), min(ncell, int(max(rows)) + 1)
c0, c1 = max(0, int(min(cols))), min(ncell, int(max(cols)) + 1)

# search the whole spring, then read each daily tile and clip to the basin window
grans = earthaccess.search_data(short_name="MOD10A1", temporal=SPRING,
                                bounding_box=(W, S, E, N))
def day_of(path):                   # MOD10A1 files are tagged .AYYYYDDD. (year + DOY)
    yyyyddd = re.search(r"\.A(\d{7})\.", os.path.basename(path)).group(1)
    return datetime.datetime.strptime(yyyyddd, "%Y%j").date()

dates, frac = [], []
for p in earthaccess.download(grans, workdir):
    p = str(p)
    if not p.endswith(".hdf"): continue
    arr = SD(p, SDC.READ).select("NDSI_Snow_Cover").get()   # (2400, 2400) uint8
    sub = arr[r0:r1, c0:c1]
    snow  = (sub >= 40) & (sub <= 100)      # snow-covered pixels
    valid = (sub >= 0)  & (sub <= 100)      # exclude cloud/fill flags
    if valid.sum() < 0.5 * sub.size: continue   # skip heavily clouded days
    dates.append(day_of(p)); frac.append(100 * snow.sum() / valid.sum())

order = np.argsort(dates)
dates = np.array(dates)[order]; frac = np.array(frac)[order]

# melt-out date = first day the snow-covered fraction drops below 50%
below = np.where(frac < 50)[0]
meltout = dates[below[0]] if len(below) else None
print(f"melt-out date (fraction first < 50%): {meltout}"
      if meltout else "snow never dropped below 50% this window")

plt.plot(dates, frac, "-o", ms=3)
plt.axhline(50, ls="--", c="grey", label="50% threshold")
if meltout is not None:
    plt.axvline(meltout, ls=":", c="r", label=f"melt-out {meltout}")
plt.ylabel("snow-covered fraction (%)"); plt.xlabel("date")
plt.legend(); plt.tight_layout(); plt.show()

# To trust the date: line it up across 20+ springs (a single year is just weather),
# and remember snow *area* is not snow *depth* — cross-check water volume with SNOTEL.

Where the data comes from

The core dataset, MOD10A1, comes from NASA’s LP DAAC archive (the data center for land products) and arrives as HDF-EOS2 tiles from the MODIS sensor. The optional add-ons live in other NASA archives: GPM IMERG rainfall and SWOT river heights each come from their own data centers. All three reach through one free Earthdata Login. The same account unlocks every NASA archive.