Concepts

Ten rules that every loader in the package follows. They are short on purpose: once you know them, a product you have never used behaves the way you expect, and adding one is a matter of filling in a catalog entry rather than re-deciding an interface.

One spatial input: aoi

Every loader’s first argument is an area of interest, and it accepts four formats:

import geopandas as gpd
import shapely

import easysnowdata as esd

esd.terrain.dem.load((-121.94, 46.72, -121.54, 46.99))  # (W, S, E, N) in EPSG:4326
esd.terrain.dem.load(shapely.box(-121.94, 46.72, -121.54, 46.99))
esd.terrain.dem.load(gpd.read_file("basin.geojson"))  # any CRS; reprojected for you
esd.terrain.dem.load(geobox)  # an odc.geo.GeoBox — exact grid

Internally all four become an easysnowdata.aoi.AOI: a target grid (GeoBox) plus a footprint (GeoDataFrame). Pass a GeoBox when you need output on a grid you already have — the loader will not resample to a grid of its own choosing.

aoi = esd.parse_aoi(gpd.read_file("basin.geojson"))
aoi.bounds, aoi.utm_crs, aoi.to_geobox(resolution=30)

Two knobs worth knowing:

clip

a flag on the AOI, not a loader keyword: esd.parse_aoi(aoi, clip=False). True by default — the result is cut to the AOI. clip=False returns the covering tiles or granules whole, which is what you want when the AOI is a point, when you are mosaicking yourself, or when you need the untouched source pixels.

Antimeridian

An AOI that crosses ±180° is split before it reaches a STAC API, which is the form those APIs expect. You do not need to do anything.

One temporal input: time

Anything pandas or STAC understands, and partial dates expand to the period they name:

esd.snow.snodas.load(aoi, "2024-03")  # the whole month
esd.snow.snodas.load(aoi, ("2024-03-01", "2024-03-15"))  # inclusive pair
esd.snow.snodas.load(aoi, "2023-10/2024-06")  # STAC syntax
esd.snow.snodas.load(aoi, "2024-03/")  # open end = now, at call time

“Now” is computed when you call, never when the module is imported, so a long-running kernel does not quietly freeze its own idea of today.

Water years are first class:

from easysnowdata.processing import add_water_year_coords

obs_ds = add_water_year_coords(obs_ds)  # adds `water_year` and `dowy` coords
obs_ds.resample(time="YS-OCT").max()  # water-year aggregation, the xarray way

dowy is day of water year, 1 on 1 October (northern hemisphere; pass hemisphere="southern" for 1 April). It is a plain integer coordinate, so it can be swapped in as a plotting axis or grouped on directly.

A southern water year starts 1 April and is named for the year it starts, so water year N is the snow season that peaks in calendar year N in either hemisphere. The station loaders (stations.load, stations.archive.load) default to hemisphere="auto": each station gets its own hemisphere’s water year from its latitude, recorded in a hemisphere coordinate. When one load spans both hemispheres, water_year and dowy become (station, time) coordinates. groupby("water_year") still works on them, while swap_dims needs one station selected first: ds.sel(station=code).swap_dims({"time": "dowy"}).

Lazy by default

Loaders return Dask-backed xarray objects and never call .compute(). A load that returns instantly has done a metadata read, not a data read; the bytes move when you compute, plot, or write.

s1_ds = esd.sar.sentinel1.load(aoi, "2024-03")  # seconds: STAC search only
s1_ds.mean("time").compute()  # now the COGs are read

chunks= means the same thing on every loader:

• left out: lazy and Dask-backed, with the source’s native chunking where it is known;

• chunks=None: no Dask. The data are read into memory before load returns, as chunks=None means in rioxarray, xarray, and odc-stac;

• anything else ({}, True, "auto", or a dict of sizes): passed to the library underneath as it is.

load(aoi, chunks=None) and load(aoi).compute() give the same result.

Vector products — basins, stations inventories, boundaries, mountain ranges, glacier outlines — are the exception: they are small, so they come back as an in-memory GeoDataFrame. What keeps them cheap is that each archive is downloaded once into the easysnowdata cache (EASYSNOWDATA_CACHE_DIR moves it) and read with the AOI pushed down, so only the intersecting features are parsed.

Subset wide stores before you chunk them

Opening a 273-variable hourly store like ARCO-ERA5 with a chunks= argument builds a Dask graph over the whole archive before you have selected anything — enough to exhaust memory. The ERA5 loader opens with chunks=None, selects, and only then chunks. Do the same with any wide cloud-native store you add.

Search and load are separate

search returns a GeoDataFrame of the items or granules that match, with the STAC properties as columns. Inspect it, filter it, and hand what survives to load. Nothing prints; nothing is downloaded.

items_gdf = esd.optical.sentinel2.search(aoi, "2024-03", cloud_cover=30)
items_gdf[["datetime", "eo:cloud_cover", "s2:mgrs_tile"]]
best_gdf = items_gdf.sort_values("eo:cloud_cover").head(3)
s2_ds = esd.optical.sentinel2.load(aoi, items=best_gdf)

What comes back

Dims

time, y, x for projected grids; time, latitude, longitude for geographic (EPSG:4326) ones. That is the odc-stac convention and the one ERA5-style data already uses.

CRS

always written with both the .rio and .odc accessors — they read the same metadata, and writing both means neither library has to guess.

Nodata

categorical products keep the source sentinel with rio.nodata set, so a land-cover class map still has integer classes; continuous products are NaN-masked with encoded_nodata preserved so the original value is recoverable. mask= overrides either way.

Attributes

only strings and numbers, never Python objects — that is what lets a result round-trip through Zarr and netCDF. Every product carries source, source_url, product_id, data_citation, license, and easysnowdata_version. ds.attrs["source"] is the source id, so it hands straight back to load(source=...); the human-readable name is in source_title.

Units

metric, always. Station SWE, snow depth, and snowfall are centimetres and precipitation millimetres, converted in the client whatever the network reports (AWDB reports inches); gridded SWE and snow depth are metres.

Vector products

a GeoDataFrame in EPSG:4326 (reprojected from the source’s CRS, 2-D even where the archive stores a constant Z) holding the features that intersect the AOI, whole rather than cut at its edge. The provenance attributes above are in gdf.attrs. A product family leads with the same columns whatever the source — name, iso3, admin_level for administrative units; rgi_id, name, area_km2, o1region for either RGI version — and keeps the source’s own columns after them. A returned GeoDataFrame is itself an AOI:

wa_gdf = esd.boundaries.admin.states(name="Washington", country="USA")
huc4_gdf = esd.hydro.basins.huc(wa_gdf, level=4)

Categorical products additionally carry CF flag attributes — flag_values, flag_meanings, flag_colors, long_name — which is what lets the plotting helpers draw a correct legend without any callable living in .attrs:

esd.plotting.categorical(esd.land.landcover.load(aoi))

Several routes to the same product

A product can be served by more than one provider, and that is treated as normal rather than exceptional. The first source listed is the default; source= picks another; the return contract does not change.

esd.terrain.dem.load(aoi)  # Planetary Computer
esd.terrain.dem.load(aoi, source="earth-search")  # AWS Open Data, same output
esd.terrain.dem.load(aoi, source="opentopography")  # OpenTopography, the 2023_1 release

Routes are not always byte-identical copies: here the three serve the 2021, 2021 and 2023_1 releases of GLO-30, and attrs["source"] and the catalog page say which one you got. OpenTopography’s catalog is a static STAC catalog rather than a STAC API — there is nothing to search — so easysnowdata.providers.opentopography downloads its item once, picks the tiles by their bounding boxes, and hands them to the same odc-stac load as every other STAC route.

This is the escape hatch when a provider has an outage, when you want to avoid an account, or when you are in a cloud region where one route is direct S3 and the other is not. Each product’s catalog page compares its routes side by side, with resolution, extent, temporal coverage, latency, and what differs.

esd.catalog.describe("copernicus-dem")

Credentials are lazy, uniform, and checked early

A product declares what it needs (requires=("earthdata",)); the provider is initialised on first use; a missing credential raises CredentialError before any network request, naming the setup steps and any credential-free alternative for that product. See Credentials.

esd.auth.status()  # a table: provider, configured?, how, needed by

Processing is pure, plotting is separate

Masking, metadata-driven scaling, baseline harmonization, dB conversion, water-year coordinates, and the local-incidence-angle computation are standalone functions in easysnowdata.processing that take and return xarray objects and do no I/O. Loaders call them for you through keyword options, but the raw product is always one call away.

What is not wrapped is band arithmetic. A normalized difference, a threshold on a snow-cover byte, or an RGB stretch is one line of xarray, and putting it behind a function name hid which bands were used and what happened to the sentinel values. The gallery writes these out every time:

s2_ds = esd.optical.sentinel2.load(aoi, "2024-03", mask="scl-default")  # convenience
raw_ds = esd.optical.sentinel2.load(aoi, "2024-03")  # nothing applied
ndsi_da = (raw_ds["green"] - raw_ds["swir16"]) / (raw_ds["green"] + raw_ds["swir16"])

easysnowdata.plotting is optional — every product plots fine with plain xarray and geopandas — but it is where the package’s map and time-series conventions live: map, categorical, and points give equal-aspect axes (a latitude-corrected aspect plus a GeographicAxesWarning for data in degrees), a colorbar matched to the map, a scale bar, a light latitude/longitude graticule, and an optional web basemap; timeseries puts calendar dates on the x axis; and label writes long name [units], square brackets always. Every piece of furniture is a keyword argument.

No import-time side effects, and logging instead of printing

Importing the package configures nothing global: GDAL and rasterio settings live in context managers around the reads that need them, and xr.set_options is never called on your behalf. Progress goes to logging.getLogger("easysnowdata") — INFO for progress, DEBUG for URLs, WARNING for fallbacks — so you choose the verbosity:

import logging

logging.getLogger("easysnowdata").setLevel(logging.DEBUG)

Two exceptions print: the one-line credential summary in interactive sessions, which makes no network request and is silenced with EASYSNOWDATA_QUIET=1, and the progress bar pooch shows while it downloads a file into the cache.

Every product has six artefacts

A product is not finished until it has a catalog entry, an offline test, a live smoke test, a health probe, a gallery example, and an upstream-watch entry. The catalog entry is what generates its docs page and its health row, so the set stays in step. Contributing turns that into a checklist.