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

What you can answer

Snow-covered area (SCA) for any AOI, daily at 500 m, from MODIS/VIIRS NDSI snow cover
Seasonal snow accumulation and melt curve — when snow peaks (typically ~April 1 in the Sierra) and when it disappears
Snow-cover duration trends — is the snow-on/snow-off season getting shorter across the record?
Snow-line elevation shifts — by combining NDSI snow cover with a DEM, where the persistent snow line sits each year
Relative year-over-year comparison — this year’s peak SCA and melt-out date vs the 2000–present climatology

What you can NOT answer with these datasets alone

Snow water equivalent (SWE) in mm directly from MODIS/VIIRS — optical sensors see extent, not depth or water content. SWE needs passive microwave (AMSR2/SMAP, coarse ~10–25 km), Daymet (model-based, North America only), airborne lidar (ASO), or ground stations (SNOTEL).
Snow under dense forest canopy — optical NDSI is blinded by canopy; microwave and modeled products fill this gap with their own large uncertainties.
Cloud-obscured days — MODIS/VIIRS see nothing through cloud. You must gap-fill (Terra+Aqua merge, temporal compositing) and honestly flag interpolated days.
The actual deliverable runoff volume — converting SWE to streamflow needs a hydrologic model plus sublimation, soil moisture, and routing; snow data is an input, not the answer.

Code template (Python, cloud-direct)

import earthaccess
import xarray as xr
import numpy as np
import rioxarray  # for NDSI raster reads via GDAL/HDF-EOS

earthaccess.login(strategy="netrc")

# AOI: Sierra Nevada (California snow-fed basins) — [W, S, E, N]
aoi = (-120.6, 36.4, -118.2, 39.3)
water_year = 2025  # Oct 2024 – Sep 2025

# 1. MODIS Terra daily snow cover (the SCA workhorse, 2000+)
mod = earthaccess.search_data(
    short_name="MOD10A1",            # Terra; pair with MYD10A1 (Aqua) to cut cloud gaps
    bounding_box=aoi,
    temporal=(f"{water_year-1}-10-01", f"{water_year}-09-30"),
)
# Open granules cloud-direct; read the NDSI_Snow_Cover band (0–100 = snow, 200–255 = flags)
files = earthaccess.open(mod)
# For each day: mask cloud/fill flags, count pixels with NDSI_Snow_Cover >= 40
#   -> snow-covered area (km^2) and fractional SCA for the AOI

# 2. VIIRS continuity (2012+) — extends the record as MODIS Terra ages out
vnp = earthaccess.search_data(
    short_name="VNP10A1",
    bounding_box=aoi,
    temporal=(f"{water_year-1}-10-01", f"{water_year}-09-30"),
)

# 3. SWE context (NOT from optical) — Daymet over North America
swe = earthaccess.search_data(
    short_name="Daymet_Daily_V4R1",  # use the 'swe' variable
    bounding_box=aoi,
    temporal=(f"{water_year-1}-10-01", f"{water_year}-09-30"),
)
# Sum/peak SWE over the AOI -> modeled basin water stored as snow (mm -> volume)

# 4. Build the seasonal curve + compare to climatology
#    - daily SCA time series for this water year
#    - 2000-present mean SCA curve (per day-of-water-year)
#    - peak SWE date and melt-out date vs the long-term mean
#    - linear trend in snow-cover duration (days/decade)

Expected output

Seasonal SCA curve: fractional snow-covered area by day-of-water-year, this year overlaid on the 2000–present mean (with min/max envelope)
Peak-snow map: maximum NDSI snow-cover extent for the AOI near April 1, with the climatological snow line marked
SWE context panel (Daymet): modeled peak basin SWE in mm and equivalent stored volume — clearly labeled as model-based, not optical
Trend readout: snow-cover duration change in days/decade and melt-out-date shift, with the caveat that ~25 years is a short hydrologic record

Caveats

NDSI threshold matters. The MOD10A1 default snow/no-snow cut is NDSI ≥ 0.4 (band value ≥ 40); document whichever you use.
Terra+Aqua merge cuts cloud gaps but introduces a ~3-hour viewing-time difference; for melt-rate estimates note it.
Collection 6.1 is current (MOD10A1.061). Don’t mix C6 and C6.1 NDSI values in one trend.
Microwave SWE saturates above ~150–200 mm SWE and fails in wet/melting snow — don’t read AMSR2/SMAP SWE as truth in deep mountain snowpack.
Daymet is model output, not an observation; its SWE is a reanalysis-style estimate and North America only.