MODIS snow cover (MOD10A1, MOD10A1F, MOD10A2)

The MODIS snow products from Terra (MOD) and Aqua (MYD) at 500 m on the sinusoidal grid, Collection 6.1: the daily MOD10A1 NDSI snow cover, the cloud-gap-filled daily MOD10A1F, and the 8-day maximum snow extent MOD10A2. Since Collection 6 the daily product reports the NDSI itself (0-100) rather than a fixed-threshold snow map, so the threshold is yours; values above 100 are sentinels (cloud, night, water, fill) carried in the CF flag_values / flag_meanings attributes of the array.

Two routes. source="nsidc" (the default) downloads the HDF-EOS2 granules through Earthdata Login and needs a GDAL with the HDF4 driver; it is the only route with the cloud-gap-filled product and the Aqua siblings. source="planetary-computer" reads the MOD10A1 / MOD10A2 COG mirrors with no credentials.

The figures show a clear day of NDSI and its binary snow map, the daily clear-sky and snow fractions of the box, an 8-day maximum extent, the cloud-gap-filled product on a day the raw product saw only cloud, and the difference between the two routes.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

import easysnowdata as esd

aoi = (-121.94, 46.72, -121.54, 46.99)  # Mount Rainier, WA
week = "2023-03-15/2023-03-22"
clear_day = "2023-03-18"
cloudy_day = "2023-03-19"

box = esd.parse_aoi(aoi)
grid = box.to_geobox(crs="utm", resolution=500)  # every map below is drawn on this grid
print(box)
print(grid)

AOI(bounds=(-121.9400, 46.7200, -121.5400, 46.9900), clip=True)
GeoBox((62, 63), Affine(500.0, 0.0, 580500.0,
       0.0, -500.0, 5205500.0), CRS('EPSG:32610'))

Where the product comes from, and what each route asks for.

for src in esd.catalog.get("modis-snow").sources:
    print(f"{src.id:22} {src.title:45} {', '.join(src.requires) or 'no account'}")

nsidc                  NSIDC (Earthdata)                             earthdata
planetary-computer     Planetary Computer                            no account

A week of the daily product from the credential-free mirror. The loader keeps the raw byte with its CF flags rather than deciding what you want, and keeps the native sinusoidal grid, which at this longitude is sheared: north is not up and a lat/lon box is a parallelogram on it. The maps are therefore drawn on the AOI’s UTM grid, with nearest-neighbour resampling so the coded bytes survive; the box is loaded with a 1 km margin so that grid is fully covered. The native array is kept for the byte comparison at the end.

daily_ds = esd.snow.modis.load(
    box.buffer(1000), week, product="MOD10A1", source="planetary-computer"
)
native_da = daily_ds["NDSI_Snow_Cover"].compute()
ndsi_da = native_da.odc.reproject(grid, resampling="nearest")
print(ndsi_da.attrs["flag_values"])
print(ndsi_da.attrs["flag_meanings"])

[200, 201, 211, 237, 239, 250, 251, 252, 254, 255]
missing_data no_decision night inland_water ocean cloud cloud cloud detector_saturated fill

The NDSI on a clear day, and the binary snow map at the heritage threshold of 40 %. valid_da drops the sentinels, so cloud and water are NaN, not 0.

scene_da = ndsi_da.sel(time=clear_day)
valid_da = scene_da <= 100
snow_da = (scene_da >= 40).where(valid_da)
snow_da = esd.processing.set_flags(
    snow_da,
    [0, 1],
    ["no snow", "snow"],
    ["#a6611a", "#2166ac"],
    long_name="snow (NDSI ≥ 40 %)",
)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
esd.plotting.map(
    scene_da.where(valid_da),
    ax=axes[0],
    cmap="Blues",
    vmin=0,
    vmax=100,
    title=f"NDSI, {clear_day}",
)
esd.plotting.categorical(snow_da, ax=axes[1], title=f"binary snow, {clear_day}")
fig.tight_layout()

How much of the box each day could see, and how much of what it saw was snow. Cloud is the story in March in the Cascades: one day saw nothing.

clear_da = (ndsi_da <= 100).mean(dim=["y", "x"])
snowy_da = ((ndsi_da >= 40) & (ndsi_da <= 100)).sum(dim=["y", "x"]) / (
    ndsi_da <= 100
).sum(dim=["y", "x"])
fractions_da = xr.concat(
    [clear_da, snowy_da],
    dim=pd.Index(
        ["clear-sky fraction", "snow fraction of clear pixels"], name="quantity"
    ),
)
fractions_da.attrs.update({"long_name": "fraction of the box", "units": "1"})

ax = esd.plotting.timeseries(fractions_da, marker="o")
ax.set_ylim(0, 1.05)
ax.figure.tight_layout()

The 8-day maximum snow extent is a classified product, so it is drawn from its class table. Its periods start on fixed days of the year (1, 9, 17, …), so the composite dated 14 March covers 14-21 March.

extent_ds = esd.snow.modis.load(
    box.buffer(1000), week, product="MOD10A2", source="planetary-computer"
)
composite_da = (
    extent_ds["Maximum_Snow_Extent"]
    .isel(time=0)
    .odc.reproject(grid, resampling="nearest")
)
print(f"composites in the window: {extent_ds['time'].values.astype('datetime64[D]')}")

ax = esd.plotting.categorical(composite_da, title="MOD10A2 maximum snow extent")
ax.figure.tight_layout()

composites in the window: ['2023-03-14' '2023-03-22']

The cloud-gap-filled product from NSIDC carries the last clear observation forward through cloudy days. On the day the raw product saw only cloud it still has a snow map.

short = f"{cloudy_day}/{cloudy_day}"
filled_ds = esd.snow.modis.load(box.buffer(1000), short, product="MOD10A1F")
cgf_da = (
    filled_ds["CGF_NDSI_Snow_Cover"]
    .sel(time=cloudy_day)
    .odc.reproject(grid, resampling="nearest")
    .compute()
)
raw_da = ndsi_da.sel(time=cloudy_day)
print(f"raw MOD10A1 valid pixels on {cloudy_day}: {int((raw_da <= 100).sum())}")
print(f"MOD10A1F valid pixels on {cloudy_day}: {int((cgf_da <= 100).sum())}")

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
esd.plotting.map(
    raw_da.where(raw_da <= 100),
    ax=axes[0],
    cmap="Blues",
    vmin=0,
    vmax=100,
    title=f"MOD10A1, {cloudy_day}",
    cbar_label="NDSI [%]",
)
esd.plotting.map(
    cgf_da.where(cgf_da <= 100),
    ax=axes[1],
    cmap="Blues",
    vmin=0,
    vmax=100,
    title=f"MOD10A1F, {cloudy_day}",
    cbar_label="NDSI [%]",
)
fig.tight_layout()

/home/runner/work/easysnowdata/easysnowdata/.pixi/envs/docs/lib/python3.13/site-packages/earthaccess/results.py:348: FutureWarning: As of version 1.0, `DataGranule.size` will be accessed as an attribute; e.g. use `DataCollection.size` **not** `DataCollection.size()`
  self["size"] = self.size()
/home/runner/work/easysnowdata/easysnowdata/.pixi/envs/docs/lib/python3.13/site-packages/earthaccess/store.py:838: FutureWarning: As of version 1.0, `DataGranule.size` will be accessed as an attribute; e.g. use `DataCollection.size` **not** `DataCollection.size()`
  total_size = round(sum(granule.size() for granule in granules) / 1024, 2)
/home/runner/work/easysnowdata/easysnowdata/.pixi/envs/docs/lib/python3.13/site-packages/rasterio/__init__.py:367: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dataset = DatasetReader(path, driver=driver, sharing=sharing, thread_safe=thread_safe, **kwargs)
raw MOD10A1 valid pixels on 2023-03-19: 0
MOD10A1F valid pixels on 2023-03-19: 3906

The two routes, same granule: NSIDC’s HDF-EOS2 file against the mirror’s COG on the clear day. On the native sinusoidal grid they land on the same pixels to within a few micrometres, so the bytes can be compared directly and should be identical; the comparison uses the native arrays and only the difference is reprojected for drawing.

archive_ds = esd.snow.modis.load(
    box.buffer(1000), f"{clear_day}/{clear_day}", product="MOD10A1"
)
nsidc_da = archive_ds["NDSI_Snow_Cover"].sel(time=clear_day).compute()
mirror_da = native_da.sel(time=clear_day)
nsidc_da = nsidc_da.reindex_like(mirror_da, method="nearest", tolerance=1)
identical = float((nsidc_da.values == mirror_da.values).mean()) * 100
print(f"bytes identical between NSIDC and Planetary Computer: {identical:.2f} %")

diff_da = (nsidc_da.astype("int16") - mirror_da.astype("int16")).where(
    (mirror_da <= 100) & (nsidc_da <= 100)
)
diff_da.attrs.update(
    {"long_name": "NDSI, NSIDC minus Planetary Computer", "units": "%"}
)
ax = esd.plotting.map(
    diff_da.odc.reproject(grid, resampling="nearest"),
    cmap="RdBu",
    vmin=-10,
    vmax=10,
    title=f"NSIDC minus Planetary Computer, {clear_day}",
    cbar_label="NDSI difference [%]",
)
ax.figure.tight_layout()
print(f"largest absolute difference: {float(np.nanmax(np.abs(diff_da.values))):.0f} %")

/home/runner/work/easysnowdata/easysnowdata/.pixi/envs/docs/lib/python3.13/site-packages/earthaccess/results.py:348: FutureWarning: As of version 1.0, `DataGranule.size` will be accessed as an attribute; e.g. use `DataCollection.size` **not** `DataCollection.size()`
  self["size"] = self.size()
/home/runner/work/easysnowdata/easysnowdata/.pixi/envs/docs/lib/python3.13/site-packages/earthaccess/store.py:838: FutureWarning: As of version 1.0, `DataGranule.size` will be accessed as an attribute; e.g. use `DataCollection.size` **not** `DataCollection.size()`
  total_size = round(sum(granule.size() for granule in granules) / 1024, 2)
/home/runner/work/easysnowdata/easysnowdata/.pixi/envs/docs/lib/python3.13/site-packages/rasterio/__init__.py:367: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dataset = DatasetReader(path, driver=driver, sharing=sharing, thread_safe=thread_safe, **kwargs)
bytes identical between NSIDC and Planetary Computer: 100.00 %
largest absolute difference: 0 %