Aggregation API #

The temporal aggregator: reduce a downloaded stack into windowed composites (daily mean, monthly sum, …). For the guide with worked examples, see Temporal aggregation.

from earthlens.core import AggregationConfig, aggregate_netcdf

aggregate= is forwarded only when both conditions hold: the backend's OUTPUT_KIND is raster or mixed — the shapes a gridded reduction is defined for — and the backend declares SUPPORTS_AGGREGATE. A vector / tabular backend is refused because the aggregator has no meaning on GeoDataFrame / DataFrame rows; a raster backend that has not wired the reducer is refused for that reason instead. Either way the refusal is a NotImplementedError raised before the backend's download runs. See Base contracts.

`AggregationConfig`#

`earthlens.core.AggregationConfig` #

Bases: BaseModel

Frozen request shape consumed by :func:aggregate_netcdf.

Carries the windowing frequency, reduction operator, and output location. Frozen + extra="forbid" so a typo in a field name (e.g. freqency=) fails loud at construction time rather than silently using the default.

Attributes:

Name	Type	Description
`freq`	`str`	Pandas offset alias defining the window. Examples: `"1D"` (daily), `"7D"` (weekly), `"1MS"` (month-start), `"QS-DEC"` (DJF/MAM/JJA/SON climatological seasons), `"AS"` (annual). Any string accepted by `pandas.Grouper(freq=...)` is valid.
`op`	`OperationLiteral`	Reduction applied within each window. `"auto"` reads `Variable.is_flux` (state→`"mean"`, flux→`"sum"`); the other values are forwarded as-is to the dispatcher.
`out_dir`	`Path \| None`	Directory the per-window GeoTIFFs are written to. Created (with parents) if absent. `None` skips the write step entirely and returns arrays in memory only.
`cell_size`	`float`	Pixel size in degrees, informational only. `0.125` for ERA5 native, `0.1` for ERA5-Land. The geotransform written to each GeoTIFF is read off the NetCDF, not from this value, and it is not encoded in the filename.
`level`	`int \| float \| None`	When the NetCDF has a `pressure_level` dimension, pin this level via :meth:`pyramids.netcdf.NetCDF.sel`. `None` (default) requires a 3-D NetCDF; pass an explicit level (e.g. `1000`) to aggregate a single 4-D layer.
`skipna`	`bool`	When `True`, the reduction is NaN-aware (`np.nanmean` etc.). `False` propagates any NaN in a window to the output.
`min_count`	`int \| None`	Minimum non-NaN samples required for a window to produce a non-NaN value. Windows with fewer samples emit NaN. `None` (default) means no minimum.

Examples:

Daily-mean defaults — only freq is required, the rest stays at sensible CDS-shaped defaults:

>>> from earthlens.aggregate import AggregationConfig
>>> cfg = AggregationConfig(freq="1D")
>>> cfg.op
'auto'
>>> cfg.skipna
True
>>> cfg.cell_size
0.125
>>> cfg.out_dir is None
True

- Monthly sum into an explicit output directory:

>>> from pathlib import Path
>>> from earthlens.aggregate import AggregationConfig
>>> cfg = AggregationConfig(
...     freq="1MS",
...     op="sum",
...     out_dir=Path("out") / "monthly",
... )
>>> cfg.freq, cfg.op
('1MS', 'sum')
>>> cfg.out_dir.name
'monthly'

- Pin a pressure level for 4-D inputs and require a minimum sample count per window:

>>> from earthlens.aggregate import AggregationConfig
>>> cfg = AggregationConfig(
...     freq="7D", op="mean", level=1000, min_count=20,
... )
>>> cfg.level, cfg.min_count
(1000, 20)

Source code in libs/core/src/earthlens/aggregate.py

class AggregationConfig(BaseModel):
    """Frozen request shape consumed by :func:`aggregate_netcdf`.

    Carries the windowing frequency, reduction operator, and output
    location. Frozen + `extra="forbid"` so a typo in a field name
    (e.g. `freqency=`) fails loud at construction time rather than
    silently using the default.

    Attributes:
        freq: Pandas offset alias defining the window. Examples:
            `"1D"` (daily), `"7D"` (weekly), `"1MS"` (month-start),
            `"QS-DEC"` (DJF/MAM/JJA/SON climatological seasons),
            `"AS"` (annual). Any string accepted by
            `pandas.Grouper(freq=...)` is valid.
        op: Reduction applied within each window. `"auto"` reads
            `Variable.is_flux` (state→`"mean"`, flux→`"sum"`); the
            other values are forwarded as-is to the dispatcher.
        out_dir: Directory the per-window GeoTIFFs are written to.
            Created (with parents) if absent. `None` skips the write
            step entirely and returns arrays in memory only.
        cell_size: Pixel size in degrees, informational only. `0.125`
            for ERA5 native, `0.1` for ERA5-Land. The geotransform
            written to each GeoTIFF is read off the NetCDF, not from
            this value, and it is not encoded in the filename.
        level: When the NetCDF has a `pressure_level` dimension, pin
            this level via :meth:`pyramids.netcdf.NetCDF.sel`. `None`
            (default) requires a 3-D NetCDF; pass an explicit level
            (e.g. `1000`) to aggregate a single 4-D layer.
        skipna: When `True`, the reduction is NaN-aware
            (`np.nanmean` etc.). `False` propagates any NaN in a
            window to the output.
        min_count: Minimum non-NaN samples required for a window to
            produce a non-NaN value. Windows with fewer samples emit
            NaN. `None` (default) means no minimum.

    Examples:
        - Daily-mean defaults — only `freq` is required, the rest
          stays at sensible CDS-shaped defaults:

            ```python
            >>> from earthlens.aggregate import AggregationConfig
            >>> cfg = AggregationConfig(freq="1D")
            >>> cfg.op
            'auto'
            >>> cfg.skipna
            True
            >>> cfg.cell_size
            0.125
            >>> cfg.out_dir is None
            True

            ```
        - Monthly sum into an explicit output directory:

            ```python
            >>> from pathlib import Path
            >>> from earthlens.aggregate import AggregationConfig
            >>> cfg = AggregationConfig(
            ...     freq="1MS",
            ...     op="sum",
            ...     out_dir=Path("out") / "monthly",
            ... )
            >>> cfg.freq, cfg.op
            ('1MS', 'sum')
            >>> cfg.out_dir.name
            'monthly'

            ```
        - Pin a pressure level for 4-D inputs and require a minimum
          sample count per window:

            ```python
            >>> from earthlens.aggregate import AggregationConfig
            >>> cfg = AggregationConfig(
            ...     freq="7D", op="mean", level=1000, min_count=20,
            ... )
            >>> cfg.level, cfg.min_count
            (1000, 20)

            ```
    """

    model_config = ConfigDict(frozen=True, extra="forbid")

    freq: str
    op: OperationLiteral = "auto"
    out_dir: Path | None = None
    cell_size: float = 0.125
    level: int | float | None = None
    skipna: bool = True
    min_count: int | None = None
    keep_arrays: bool = True

`aggregate_netcdf`#

`earthlens.core.aggregate_netcdf(nc_path, var_info, config, date_range=None)` #

Slice a CDS-shaped NetCDF into per-window aggregated outputs.

Reads the NetCDF, groups its time axis by config.freq, reduces each group with config.op, and (when config.out_dir is set) writes one GeoTIFF per window. Returns the per-window arrays alongside their timestamps and output paths so callers can chain further processing without re-opening the files.

Parameters:

Name	Type	Description	Default
`nc_path`	`Path \| str`	Path to the NetCDF on disk.	required
`var_info`	`Variable`	Catalog row for the variable being aggregated. Used to pick the variable from the NetCDF (`var_info.nc_variable`), seed the output filename (`var_info.cds_variable`, plus `dataset_id` for a curated override — see :func:`_output_stem`), and resolve `op="auto"` (`var_info.is_flux`).	required
`config`	`AggregationConfig`	Frozen :class:`AggregationConfig` describing the window, reduction, and output location.	required
`date_range`	`tuple[Any, Any] \| None`	Optional inclusive `(start, end)` bounds. When set, samples outside the range are dropped before windowing, so the output stays faithful to the requested span even when the source cube over-covers it (a CDS `year`/`month`/`day` cross-product pulls stray dates for a window spanning month boundaries). `None` keeps every sample.	`None`

Returns:

Type	Description
`list[tuple[Timestamp, ndarray \| None, Path \| None]]`	list[tuple[pd.Timestamp, np.ndarray \| None, Path \| None]]: One
`list[tuple[Timestamp, ndarray \| None, Path \| None]]`	entry per window. The first item is the window's left-edge
`list[tuple[Timestamp, ndarray \| None, Path \| None]]`	timestamp; the second is the reduced 2-D array — `None` only when
`list[tuple[Timestamp, ndarray \| None, Path \| None]]`	the request set `keep_arrays=False` and the window was written to
`list[tuple[Timestamp, ndarray \| None, Path \| None]]`	disk; the third is the GeoTIFF path (or `None` when
`list[tuple[Timestamp, ndarray \| None, Path \| None]]`	`config.out_dir` was `None`).

Raises:

Type	Description
`KeyError`	If the NetCDF has no recognised time variable (`valid_time` / `time`); see :func:`_read_time_axis`.
`ValueError`	If `config.level` is set but the NetCDF has no pressure-level dimension, or vice versa; see :func:`_resolve_pressure_level`. Also raised by pandas when `config.freq` is not a recognised offset alias.

`iter_aggregate_netcdf`#

Streams one reduced window at a time instead of materialising the whole cube — this is what keeps memory bounded on a long time series.

`earthlens.core.iter_aggregate_netcdf(nc_path, var_info, config, date_range=None)` #

Yield one :class:AggregatedWindow per time window, streaming.

The streaming counterpart to :func:aggregate_netcdf, and the implementation it is built on. Two properties make it usable on cubes that do not fit in memory:

Only one window is resident at a time. The time steps for the current window are read band by band and stacked; the whole (time, y, x) cube is never materialised. A ten-year hourly ERA5 request over Europe at 0.25° is ~33.6 GB as one array but ~9 MB per daily window.
Windows are not accumulated. Each is yielded and then dropped, so a caller that writes and discards holds nothing. Pair it with keep_arrays=False to drop the reduced array too once it is on disk.

Every handle opened — the container, the variable subset, and the level-pinned view — is closed when the generator finishes or is abandoned, so the file can be deleted or overwritten straight after. Closing the container alone is not sufficient: the variable subset holds its own handle, and one left open keeps a Windows lock on the file.

Parameters:

Name	Type	Description	Default
`nc_path`	`Path \| str`	Path to the NetCDF on disk.	required
`var_info`	`Variable`	Catalog row for the variable being aggregated. Used to pick the variable from the NetCDF (`var_info.nc_variable`), seed the output filename (`_output_stem` — `var_info.cds_variable`, plus `dataset_id` for a curated override), and resolve `op="auto"` (`var_info.is_flux`).	required
`config`	`AggregationConfig`	Frozen :class:`AggregationConfig` describing the window, reduction, output location, and whether to retain arrays.	required
`date_range`	`tuple[Any, Any] \| None`	Optional inclusive `(start, end)` bounds. Samples outside the range are dropped before windowing, so the streamed windows match the requested span even when the source cube over-covers it; a window left empty by the filter is skipped. `None` keeps every sample. See :func:`_date_range_mask`.	`None`

Yields:

Name	Type	Description
`AggregatedWindow`	`AggregatedWindow`	One per window, in time order.

Raises:

Type	Description
`KeyError`	If the NetCDF has no recognised time variable (`valid_time` / `time`); see :func:`_read_time_axis`.
`ValueError`	If `config.level` is set but the NetCDF has no pressure-level dimension, or vice versa; see :func:`_resolve_pressure_level`. Also raised by pandas when `config.freq` is not a recognised offset alias.

`AggregatedWindow`#

`earthlens.core.AggregatedWindow` `dataclass` #

One reduced time window: its label, its array, and where it was written.

Yielded by :func:iter_aggregate_netcdf. array is None when the request set keep_arrays=False and the window was written to disk — the point of that mode is that a long run does not accumulate every window in memory alongside the GeoTIFFs it has already produced.

Comparison is left as identity (eq=False). A generated __eq__ would compare the array field with ==, which for a numpy array yields an elementwise array and then raises ValueError: truth value ... ambiguous; the matching __hash__ would raise TypeError because an ndarray is unhashable. Compare the fields you actually care about instead.

Attributes:

Name	Type	Description
`label`	`Timestamp`	The window's left-edge timestamp.
`array`	`ndarray \| None`	The reduced 2-D array, or `None` when the request asked not to retain it.
`path`	`Path \| None`	The written GeoTIFF, or `None` when `out_dir` was `None`.

Examples:

A window keeps both its array and its path when it was written:

>>> import numpy as np
>>> import pandas as pd
>>> from pathlib import Path
>>> from earthlens.aggregate import AggregatedWindow
>>> window = AggregatedWindow(
...     label=pd.Timestamp("2020-01-01"),
...     array=np.array([[1.0, 2.0]]),
...     path=Path("out/t2m_reanalysis-era5-single-levels_1D_20200101.tif"),
... )
>>> window.label.strftime("%Y-%m-%d")
'2020-01-01'
>>> float(window.array.mean())
1.5
>>> window.path.name
't2m_reanalysis-era5-single-levels_1D_20200101.tif'

A discarded array leaves only the label and the path to read back:

>>> import pandas as pd
>>> from pathlib import Path
>>> from earthlens.aggregate import AggregatedWindow
>>> window = AggregatedWindow(
...     label=pd.Timestamp("2020-02-01"),
...     array=None,
...     path=Path("out/t2m_reanalysis-era5-single-levels_1D_20200201.tif"),
... )
>>> window.array is None
True
>>> window.path.name
't2m_reanalysis-era5-single-levels_1D_20200201.tif'

Source code in libs/core/src/earthlens/aggregate.py

@dataclass(frozen=True, eq=False)
class AggregatedWindow:
    """One reduced time window: its label, its array, and where it was written.

    Yielded by :func:`iter_aggregate_netcdf`. `array` is `None` when the
    request set `keep_arrays=False` and the window was written to disk — the
    point of that mode is that a long run does not accumulate every window in
    memory alongside the GeoTIFFs it has already produced.

    Comparison is left as identity (`eq=False`). A generated `__eq__` would
    compare the `array` field with `==`, which for a numpy array yields an
    elementwise array and then raises `ValueError: truth value ... ambiguous`;
    the matching `__hash__` would raise `TypeError` because an ndarray is
    unhashable. Compare the fields you actually care about instead.

    Attributes:
        label: The window's left-edge timestamp.
        array: The reduced 2-D array, or `None` when the request asked not to
            retain it.
        path: The written GeoTIFF, or `None` when `out_dir` was `None`.

    Examples:
        - A window keeps both its array and its path when it was written:
            ```python
            >>> import numpy as np
            >>> import pandas as pd
            >>> from pathlib import Path
            >>> from earthlens.aggregate import AggregatedWindow
            >>> window = AggregatedWindow(
            ...     label=pd.Timestamp("2020-01-01"),
            ...     array=np.array([[1.0, 2.0]]),
            ...     path=Path("out/t2m_reanalysis-era5-single-levels_1D_20200101.tif"),
            ... )
            >>> window.label.strftime("%Y-%m-%d")
            '2020-01-01'
            >>> float(window.array.mean())
            1.5
            >>> window.path.name
            't2m_reanalysis-era5-single-levels_1D_20200101.tif'

            ```
        - A discarded array leaves only the label and the path to read back:
            ```python
            >>> import pandas as pd
            >>> from pathlib import Path
            >>> from earthlens.aggregate import AggregatedWindow
            >>> window = AggregatedWindow(
            ...     label=pd.Timestamp("2020-02-01"),
            ...     array=None,
            ...     path=Path("out/t2m_reanalysis-era5-single-levels_1D_20200201.tif"),
            ... )
            >>> window.array is None
            True
            >>> window.path.name
            't2m_reanalysis-era5-single-levels_1D_20200201.tif'

            ```
    """

    label: pd.Timestamp
    array: np.ndarray | None
    path: Path | None

Reduction helpers #

Public since 0.12.0 (previously _reduce and _window_groups).

`earthlens.aggregate.reduce_time_axis(arr, op, skipna, min_count)` #

Reduce a (time, lat, lon) slice along axis 0 with the named op.

Dispatches op to the matching numpy reducer (np.nanmean etc. when skipna=True, plain np.mean etc. when skipna=False), then masks pixels whose non-NaN sample count falls below min_count.

op="auto" is not accepted here — aggregate_netcdf resolves auto to a concrete operator before calling this helper. Passing "auto" raises KeyError to surface the mistake at the call site.

Parameters:

Name	Type	Description	Default
`arr`	`ndarray`	Array to reduce. The first axis is collapsed; the remaining axes pass through unchanged. Typically `(N_in_window, lat, lon)`.	required
`op`	`str`	One of `"mean" / "sum" / "min" / "max" / "std"`. Resolved to a numpy reducer via the dispatch table.	required
`skipna`	`bool`	When `True`, the NaN-aware reducer is used (`np.nanmean` etc.); when `False`, the strict variant is used and any NaN in a window propagates to the output.	required
`min_count`	`int \| None`	When set, pixels with fewer than this many non-NaN samples along axis 0 emit NaN regardless of the reduction result. `None` disables the floor.	required

Returns:

Type	Description
`ndarray`	np.ndarray: Reduced array with axis 0 collapsed.

Raises:

Type	Description
`KeyError`	If `op` is not in the dispatch table (in particular, `"auto"` is rejected — resolve it to a concrete op first).

Examples:

NaN-aware mean over the time axis:

>>> import numpy as np
>>> from earthlens.aggregate import reduce_time_axis
>>> arr = np.array([[[1.0, 2.0]], [[3.0, np.nan]], [[5.0, 6.0]]])
>>> reduce_time_axis(arr, op="mean", skipna=True, min_count=None).tolist()
[[3.0, 4.0]]

- Strict mean propagates NaN when skipna=False:

>>> import numpy as np
>>> from earthlens.aggregate import reduce_time_axis
>>> arr = np.array([[[1.0, np.nan]], [[3.0, 4.0]]])
>>> result = reduce_time_axis(arr, op="mean", skipna=False, min_count=None)
>>> bool(np.isnan(result[0, 1])), float(result[0, 0])
(True, 2.0)

- min_count masks under-sampled pixels:

>>> import numpy as np
>>> from earthlens.aggregate import reduce_time_axis
>>> arr = np.array([[[1.0, np.nan]], [[2.0, np.nan]]])
>>> result = reduce_time_axis(arr, op="mean", skipna=True, min_count=2)
>>> float(result[0, 0]), bool(np.isnan(result[0, 1]))
(1.5, True)

Source code in libs/core/src/earthlens/aggregate.py

def reduce_time_axis(
    arr: np.ndarray,
    op: str,
    skipna: bool,
    min_count: int | None,
) -> np.ndarray:
    """Reduce a `(time, lat, lon)` slice along axis 0 with the named op.

    Dispatches `op` to the matching numpy reducer (`np.nanmean` etc.
    when `skipna=True`, plain `np.mean` etc. when `skipna=False`),
    then masks pixels whose non-NaN sample count falls below
    `min_count`.

    `op="auto"` is **not** accepted here — `aggregate_netcdf` resolves
    `auto` to a concrete operator before calling this helper. Passing
    `"auto"` raises `KeyError` to surface the mistake at the call site.

    Args:
        arr: Array to reduce. The first axis is collapsed; the
            remaining axes pass through unchanged. Typically
            `(N_in_window, lat, lon)`.
        op: One of `"mean" / "sum" / "min" / "max" / "std"`. Resolved
            to a numpy reducer via the dispatch table.
        skipna: When `True`, the NaN-aware reducer is used
            (`np.nanmean` etc.); when `False`, the strict variant is
            used and any NaN in a window propagates to the output.
        min_count: When set, pixels with fewer than this many non-NaN
            samples along axis 0 emit NaN regardless of the reduction
            result. `None` disables the floor.

    Returns:
        np.ndarray: Reduced array with axis 0 collapsed.

    Raises:
        KeyError: If `op` is not in the dispatch table (in particular,
            `"auto"` is rejected — resolve it to a concrete op first).

    Examples:
        - NaN-aware mean over the time axis:

            ```python
            >>> import numpy as np
            >>> from earthlens.aggregate import reduce_time_axis
            >>> arr = np.array([[[1.0, 2.0]], [[3.0, np.nan]], [[5.0, 6.0]]])
            >>> reduce_time_axis(arr, op="mean", skipna=True, min_count=None).tolist()
            [[3.0, 4.0]]

            ```
        - Strict mean propagates NaN when `skipna=False`:

            ```python
            >>> import numpy as np
            >>> from earthlens.aggregate import reduce_time_axis
            >>> arr = np.array([[[1.0, np.nan]], [[3.0, 4.0]]])
            >>> result = reduce_time_axis(arr, op="mean", skipna=False, min_count=None)
            >>> bool(np.isnan(result[0, 1])), float(result[0, 0])
            (True, 2.0)

            ```
        - `min_count` masks under-sampled pixels:

            ```python
            >>> import numpy as np
            >>> from earthlens.aggregate import reduce_time_axis
            >>> arr = np.array([[[1.0, np.nan]], [[2.0, np.nan]]])
            >>> result = reduce_time_axis(arr, op="mean", skipna=True, min_count=2)
            >>> float(result[0, 0]), bool(np.isnan(result[0, 1]))
            (1.5, True)

            ```
    """
    table = _REDUCERS_SKIPNA if skipna else _REDUCERS_STRICT
    if op not in table:
        raise KeyError(
            f"unknown reduction op {op!r}; expected one of "
            f"{sorted(table)!r} (resolve 'auto' before calling reduce_time_axis)"
        )
    reducer = table[op]
    result = reducer(arr, axis=0)
    if min_count is not None:
        non_nan_count = np.count_nonzero(~np.isnan(arr), axis=0)
        result = np.where(non_nan_count >= min_count, result, np.nan)
    return np.asarray(result)

`earthlens.aggregate.window_groups(time_axis, freq)` #

Yield (window_label, mask) pairs that bucket time_axis by freq.

Builds a pandas.Series indexed by time_axis and groups it with pandas.Grouper(freq=freq). Each group's index gives the timestamps belonging to that window; the boolean mask is built by membership against time_axis so callers can use it to slice a numpy array along its first axis.

Empty groups (windows with no samples) are silently skipped — aggregate_netcdf doesn't write a GeoTIFF for a window it has no data for.

Parameters:

Name	Type	Description	Default
`time_axis`	`DatetimeIndex`	Time coordinate as a :class:`pandas.DatetimeIndex`. Typically the result of :func:`_read_time_axis`.	required
`freq`	`str`	Pandas offset alias (`"1D"`, `"7D"`, `"1MS"`, `"QS-DEC"`, `"AS"`, ...). Anything :class:`pandas.Grouper` accepts.	required

Yields:

Type	Description
`Timestamp`	tuple[pd.Timestamp, np.ndarray]: For each non-empty window:
`ndarray`	the group key (window's left-edge timestamp) paired with a
`tuple[Timestamp, ndarray]`	boolean mask of length `len(time_axis)`.

Examples:

Group four 6-hourly slots into one daily window:

>>> import pandas as pd
>>> from earthlens.aggregate import window_groups
>>> idx = pd.date_range("2022-01-01", periods=4, freq="6h")
>>> windows = list(window_groups(idx, "1D"))
>>> len(windows)
1
>>> label, mask = windows[0]
>>> label
Timestamp('2022-01-01 00:00:00')
>>> mask.tolist()
[True, True, True, True]

- Group two days of 6-hourly samples into two daily windows:

>>> import pandas as pd
>>> from earthlens.aggregate import window_groups
>>> idx = pd.date_range("2022-01-01", periods=8, freq="6h")
>>> [label.strftime("%Y-%m-%d") for label, _ in window_groups(idx, "1D")]
['2022-01-01', '2022-01-02']

Source code in libs/core/src/earthlens/aggregate.py

def window_groups(
    time_axis: pd.DatetimeIndex,
    freq: str,
) -> Iterator[tuple[pd.Timestamp, np.ndarray]]:
    """Yield `(window_label, mask)` pairs that bucket `time_axis` by `freq`.

    Builds a `pandas.Series` indexed by `time_axis` and groups it
    with `pandas.Grouper(freq=freq)`. Each group's `index` gives the
    timestamps belonging to that window; the boolean mask is built
    by membership against `time_axis` so callers can use it to slice
    a numpy array along its first axis.

    Empty groups (windows with no samples) are silently skipped —
    `aggregate_netcdf` doesn't write a GeoTIFF for a window it has
    no data for.

    Args:
        time_axis: Time coordinate as a :class:`pandas.DatetimeIndex`.
            Typically the result of :func:`_read_time_axis`.
        freq: Pandas offset alias (`"1D"`, `"7D"`, `"1MS"`, `"QS-DEC"`,
            `"AS"`, ...). Anything :class:`pandas.Grouper` accepts.

    Yields:
        tuple[pd.Timestamp, np.ndarray]: For each non-empty window:
        the group key (window's left-edge timestamp) paired with a
        boolean mask of length `len(time_axis)`.

    Examples:
        - Group four 6-hourly slots into one daily window:

            ```python
            >>> import pandas as pd
            >>> from earthlens.aggregate import window_groups
            >>> idx = pd.date_range("2022-01-01", periods=4, freq="6h")
            >>> windows = list(window_groups(idx, "1D"))
            >>> len(windows)
            1
            >>> label, mask = windows[0]
            >>> label
            Timestamp('2022-01-01 00:00:00')
            >>> mask.tolist()
            [True, True, True, True]

            ```
        - Group two days of 6-hourly samples into two daily windows:

            ```python
            >>> import pandas as pd
            >>> from earthlens.aggregate import window_groups
            >>> idx = pd.date_range("2022-01-01", periods=8, freq="6h")
            >>> [label.strftime("%Y-%m-%d") for label, _ in window_groups(idx, "1D")]
            ['2022-01-01', '2022-01-02']

            ```
    """
    indexer = pd.Series(np.arange(len(time_axis)), index=time_axis)
    timestamps = pd.Index(time_axis)
    for window_label, group in indexer.groupby(pd.Grouper(freq=freq)):
        if group.empty:
            continue
        mask = np.asarray(timestamps.isin(group.index))
        yield window_label, mask