Skip to content

climate-pipeline

The modules appear in pipeline order: normalize a raw source, fetch periods from it, ingest them one commit at a time, derive indices, and publish.

normalize

Normalizing source data into the service's own conventions.

Real sources disagree about everything: dimension names (lat/lon vs latitude/longitude), units (Kelvin vs Celsius, m vs mm), and axis direction (south-up vs north-up). open-climate-service resolves all of that at ingest so that every stored dataset looks identical to everything downstream, which is what makes its API uniform across sources.

This module is that step: source in, canonical (time, y, x) in degrees Celsius or millimetres out.

Functions:

rename_dims(ds)

Rename whatever the source calls its axes to (time, y, x).

Parameters:

Name Type Description Default
ds Dataset

A dataset using any of the recognized dimension spellings.

 required

Returns:

Type Description
Dataset

The dataset with its dimensions and coordinates renamed.

Raises:

Type Description
ValueError

If a spatial or time axis cannot be identified.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/normalize.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def rename_dims(ds: xr.Dataset) -> xr.Dataset:
    """Rename whatever the source calls its axes to ``(time, y, x)``.

    Args:
        ds: A dataset using any of the recognized dimension spellings.

    Returns:
        The dataset with its dimensions and coordinates renamed.

    Raises:
        ValueError: If a spatial or time axis cannot be identified.
    """
    mapping: dict[str, str] = {}
    for aliases, target in ((TIME_ALIASES, "time"), (Y_ALIASES, "y"), (X_ALIASES, "x")):
        found = next((name for name in aliases if name in ds.dims or name in ds.coords), None)
        if found is None:
            raise ValueError(f"no dimension matching {target!r} among {tuple(ds.dims)}")
        if found != target:
            mapping[found] = target
    return ds.rename(mapping) if mapping else ds

orient_north_up(ds)

Ensure the y axis descends, so row 0 is the northernmost.

GeoZarr places a raster with an affine whose y step is negative for a north-up grid. A south-up source silently renders upside down, so the orientation is fixed here rather than trusted.

Parameters:

Name Type Description Default
ds Dataset

A dataset with a y coordinate.

 required

Returns:

Type Description
Dataset

The dataset, reversed along y if it was ascending.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/normalize.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def orient_north_up(ds: xr.Dataset) -> xr.Dataset:
    """Ensure the y axis descends, so row 0 is the northernmost.

    GeoZarr places a raster with an affine whose y step is negative for a
    north-up grid. A south-up source silently renders upside down, so the
    orientation is fixed here rather than trusted.

    Args:
        ds: A dataset with a ``y`` coordinate.

    Returns:
        The dataset, reversed along y if it was ascending.
    """
    if "y" not in ds.coords or ds.sizes.get("y", 0) < 2:
        return ds
    if float(ds["y"][0]) < float(ds["y"][-1]):
        return ds.isel(y=slice(None, None, -1))
    return ds

convert_units(ds)

Convert known variables to the service's canonical units.

Kelvin becomes Celsius and metres of precipitation become millimetres. Attributes are inert in xarray, so the units attribute is rewritten by hand -- forgetting that is how a dataset ends up labelled K while holding Celsius.

Parameters:

Name Type Description Default
ds Dataset

A dataset whose variables carry a units attribute.

 required

Returns:

Type Description
Dataset

A new dataset with converted values and corrected unit attributes.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/normalize.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def convert_units(ds: xr.Dataset) -> xr.Dataset:
    """Convert known variables to the service's canonical units.

    Kelvin becomes Celsius and metres of precipitation become millimetres.
    Attributes are inert in xarray, so the ``units`` attribute is rewritten by
    hand -- forgetting that is how a dataset ends up labelled ``K`` while
    holding Celsius.

    Args:
        ds: A dataset whose variables carry a ``units`` attribute.

    Returns:
        A new dataset with converted values and corrected unit attributes.
    """
    out = ds.copy()
    for name, var in ds.data_vars.items():
        units = str(var.attrs.get("units", "")).strip()
        if units in ("K", "kelvin", "Kelvin"):
            out[name] = var - 273.15
            out[name].attrs = {**var.attrs, "units": "degC"}
        elif units in ("m", "metre", "meters") and name == "tp":
            out[name] = var * 1000.0
            out[name].attrs = {**var.attrs, "units": "mm"}
        else:
            out[name].attrs = dict(var.attrs)
    return out

sort_time(ds)

Sort along time and drop duplicate timestamps, keeping the last.

A re-fetched period arrives with timestamps the store may already hold. Appending it blindly produces a store with duplicate coordinates that cannot be indexed sanely, so duplicates are resolved here.

Parameters:

Name Type Description Default
ds Dataset

A dataset with a time dimension.

 required

Returns:

Type Description
Dataset

The dataset sorted by time with unique timestamps.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/normalize.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def sort_time(ds: xr.Dataset) -> xr.Dataset:
    """Sort along time and drop duplicate timestamps, keeping the last.

    A re-fetched period arrives with timestamps the store may already hold.
    Appending it blindly produces a store with duplicate coordinates that
    cannot be indexed sanely, so duplicates are resolved here.

    Args:
        ds: A dataset with a ``time`` dimension.

    Returns:
        The dataset sorted by time with unique timestamps.
    """
    if "time" not in ds.dims:
        return ds
    ds = ds.sortby("time")
    index = ds.indexes["time"]
    if index.has_duplicates:
        keep = ~index.duplicated(keep="last")
        ds = ds.isel(time=np.flatnonzero(keep))
    return ds

normalize(ds)

Run the full normalization pipeline on a source dataset.

Parameters:

Name Type Description Default
ds Dataset

Raw source data in whatever conventions it arrived with.

 required

Returns:

Type Description
Dataset

A dataset with dims (time, y, x), a north-up y axis, canonical
Dataset

units, and a sorted, duplicate-free time axis.

Raises:

Type Description
ValueError

If the dataset's axes cannot be identified.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/normalize.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def normalize(ds: xr.Dataset) -> xr.Dataset:
    """Run the full normalization pipeline on a source dataset.

    Args:
        ds: Raw source data in whatever conventions it arrived with.

    Returns:
        A dataset with dims ``(time, y, x)``, a north-up y axis, canonical
        units, and a sorted, duplicate-free time axis.

    Raises:
        ValueError: If the dataset's axes cannot be identified.
    """
    ds = rename_dims(ds)
    ds = orient_north_up(ds)
    ds = convert_units(ds)
    ds = sort_time(ds)
    return ds.transpose("time", "y", "x", ...)

sources

Synthetic data sources, deliberately messy.

Each source stands in for a real one and arrives in its own conventions, so the normalization step has something to actually do. Periods are enumerated the way open-climate-service's streaming ingest does it: the plugin lists the periods it can supply, and the framework fetches them one at a time.

Classes

Period dataclass

One ingestable period.

Attributes:

Name Type Description
period_id str

Stable identifier, such as "2024-01".
start str

First day of the period, as an ISO date string.
days int

Number of daily steps in the period.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/sources.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
@dataclass(frozen=True)
class Period:
    """One ingestable period.

    Attributes:
        period_id: Stable identifier, such as ``"2024-01"``.
        start: First day of the period, as an ISO date string.
        days: Number of daily steps in the period.
    """

    period_id: str
    start: str
    days: int

Functions:

enumerate_periods(year=2024, months=6)

List the monthly periods a source can supply.

Parameters:

Name Type Description Default
year int

Calendar year to enumerate.

 2024
months int

How many months from January to include; must be 1..12.

 6

Returns:

Name Type Description
One list[Period]

class:Period per month, in chronological order.

Raises:

Type Description
ValueError

If months is outside 1..12.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/sources.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def enumerate_periods(year: int = 2024, months: int = 6) -> list[Period]:
    """List the monthly periods a source can supply.

    Args:
        year: Calendar year to enumerate.
        months: How many months from January to include; must be 1..12.

    Returns:
        One :class:`Period` per month, in chronological order.

    Raises:
        ValueError: If months is outside 1..12.
    """
    if not 1 <= months <= 12:
        raise ValueError(f"months must be between 1 and 12, got {months}")
    periods: list[Period] = []
    for month in range(1, months + 1):
        start = pd.Timestamp(year=year, month=month, day=1)
        periods.append(
            Period(
                period_id=f"{year}-{month:02d}",
                start=start.strftime("%Y-%m-%d"),
                days=int(start.days_in_month),
            )
        )
    return periods

fetch_temperature(period, ny=24, nx=24, seed=0)

Fetch one period of temperature, in a deliberately awkward source format.

This source publishes Kelvin on a south-up lat/lon grid -- exactly the kind of thing normalization exists to fix.

Parameters:

Name Type Description Default
period Period

The period to fetch.

 required
ny int

Grid height.

 24
nx int

Grid width.

 24
seed int

Base seed; the period id is mixed in so periods differ.

 0

Returns:

Type Description
Dataset

A dataset with variable t2m in Kelvin, dims (time, lat, lon).
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/sources.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def fetch_temperature(period: Period, ny: int = 24, nx: int = 24, seed: int = 0) -> xr.Dataset:
    """Fetch one period of temperature, in a deliberately awkward source format.

    This source publishes Kelvin on a south-up ``lat``/``lon`` grid -- exactly
    the kind of thing normalization exists to fix.

    Args:
        period: The period to fetch.
        ny: Grid height.
        nx: Grid width.
        seed: Base seed; the period id is mixed in so periods differ.

    Returns:
        A dataset with variable ``t2m`` in Kelvin, dims ``(time, lat, lon)``.
    """
    rng = np.random.default_rng(seed + _period_seed(period))
    time = pd.date_range(period.start, periods=period.days, freq="D")
    lat, lon = _grid(ny, nx, ascending_y=True)

    day_of_year = time.dayofyear.to_numpy().reshape(-1, 1, 1)
    seasonal = 3.0 * np.sin(2 * np.pi * day_of_year / 365.25)
    gradient = np.linspace(-1.5, 1.5, ny).reshape(1, ny, 1)
    kelvin = 273.15 + 27.0 + seasonal + gradient + rng.normal(0.0, 0.7, size=(period.days, ny, nx))

    return xr.DataArray(
        kelvin,
        dims=("time", "lat", "lon"),
        coords={"time": time, "lat": lat, "lon": lon},
        name="t2m",
        attrs={"units": "K", "long_name": "2 metre temperature"},
    ).to_dataset()

fetch_precipitation(period, ny=24, nx=24, seed=1)

Fetch one period of precipitation, in metres on a south-up grid.

Parameters:

Name Type Description Default
period Period

The period to fetch.

 required
ny int

Grid height.

 24
nx int

Grid width.

 24
seed int

Base seed; the period id is mixed in so periods differ.

 1

Returns:

Type Description
Dataset

A dataset with variable tp in metres, dims (time, lat, lon).
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/sources.py 
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def fetch_precipitation(period: Period, ny: int = 24, nx: int = 24, seed: int = 1) -> xr.Dataset:
    """Fetch one period of precipitation, in metres on a south-up grid.

    Args:
        period: The period to fetch.
        ny: Grid height.
        nx: Grid width.
        seed: Base seed; the period id is mixed in so periods differ.

    Returns:
        A dataset with variable ``tp`` in metres, dims ``(time, lat, lon)``.
    """
    rng = np.random.default_rng(seed + _period_seed(period))
    time = pd.date_range(period.start, periods=period.days, freq="D")
    lat, lon = _grid(ny, nx, ascending_y=True)

    # West African rainfall peaks mid-year; zero-inflated, like the real thing.
    month = time.month.to_numpy().reshape(-1, 1, 1)
    wetness = np.clip(np.sin((month - 2) / 12 * np.pi), 0.05, None)
    wet = rng.random((period.days, ny, nx)) < wetness
    amounts_mm = rng.gamma(2.0, 5.0, size=(period.days, ny, nx))
    metres = np.where(wet, amounts_mm, 0.0) / 1000.0

    return xr.DataArray(
        metres,
        dims=("time", "lat", "lon"),
        coords={"time": time, "lat": lat, "lon": lon},
        name="tp",
        attrs={"units": "m", "long_name": "total precipitation"},
    ).to_dataset()

ingest

Streaming ingest: one period at a time, committed as it lands.

This is the heart of the open-climate-service ingestion contract. The source enumerates periods; the framework fetches each one, normalizes it, appends it to the store, and commits. Because every period is its own transaction, an interrupted ingest leaves a store that is complete up to the last commit -- never half a period -- and resuming is a matter of asking the store what it already holds.

Classes

IngestReport dataclass

What one ingest run did.

Attributes:

Name Type Description
ingested list[str]

Period ids written during this run.
skipped list[str]

Period ids already present and therefore skipped.
failed dict[str, str]

Period ids that raised, with the error message.
snapshots dict[str, str]

Snapshot id produced by each ingested period.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/ingest.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
@dataclass
class IngestReport:
    """What one ingest run did.

    Attributes:
        ingested: Period ids written during this run.
        skipped: Period ids already present and therefore skipped.
        failed: Period ids that raised, with the error message.
        snapshots: Snapshot id produced by each ingested period.
    """

    ingested: list[str] = field(default_factory=list)
    skipped: list[str] = field(default_factory=list)
    failed: dict[str, str] = field(default_factory=dict)
    snapshots: dict[str, str] = field(default_factory=dict)

    @property
    def total(self) -> int:
        """Number of periods considered in this run."""
        return len(self.ingested) + len(self.skipped) + len(self.failed)
Attributes
total property

Number of periods considered in this run.

Functions:

chunking_for(ds)

Choose chunk sizes for a dataset, capping the spatial dimensions.

Parameters:

Name Type Description Default
ds Dataset

The dataset about to be written.

 required

Returns:

Type Description
dict[str, int]

A chunk mapping suitable for Dataset.chunk.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/ingest.py 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def chunking_for(ds: xr.Dataset) -> dict[str, int]:
    """Choose chunk sizes for a dataset, capping the spatial dimensions.

    Args:
        ds: The dataset about to be written.

    Returns:
        A chunk mapping suitable for ``Dataset.chunk``.
    """
    chunks = {"time": min(TIME_CHUNK, int(ds.sizes.get("time", 1)))}
    for dim in ("y", "x"):
        if dim in ds.sizes:
            chunks[dim] = min(SPATIAL_CHUNK_CAP, int(ds.sizes[dim]))
    return chunks

committed_periods(repo, *, period_type='month')

Return the period ids already committed to a store.

The store's time coordinate is authoritative: whatever is committed is what exists, regardless of what any external bookkeeping claims. That is what makes resume safe after a crash.

Parameters:

Name Type Description Default
repo Any

An icechunk repository.

 required
period_type str

Granularity of the period ids; only "month" is currently supported.

 'month'

Returns:

Type Description
set[str]

The set of period ids present, empty if the store has no data yet.

Raises:

Type Description
ValueError

If period_type is not supported.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/ingest.py 
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def committed_periods(repo: Any, *, period_type: str = "month") -> set[str]:
    """Return the period ids already committed to a store.

    The store's time coordinate is authoritative: whatever is committed is
    what exists, regardless of what any external bookkeeping claims. That is
    what makes resume safe after a crash.

    Args:
        repo: An icechunk repository.
        period_type: Granularity of the period ids; only ``"month"`` is
            currently supported.

    Returns:
        The set of period ids present, empty if the store has no data yet.

    Raises:
        ValueError: If period_type is not supported.
    """
    if period_type != "month":
        raise ValueError(f"unsupported period_type: {period_type!r}")
    try:
        ds = xr.open_zarr(repo.readonly_session("main").store, consolidated=False)
    except Exception:
        # A repository with no committed data yet has no group to open.
        return set()
    if "time" not in ds.coords or ds.sizes.get("time", 0) == 0:
        return set()
    stamps = pd.DatetimeIndex(ds["time"].values)
    return {f"{ts.year}-{ts.month:02d}" for ts in stamps}

ingest_period(repo, period, fetch)

Fetch, normalize, and commit a single period.

Parameters:

Name Type Description Default
repo Any

An icechunk repository to write into.

 required
period Period

The period to ingest.

 required
fetch Fetcher

Callable returning the raw source dataset for the period.

 required

Returns:

Type Description
str

The snapshot id of the commit.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/ingest.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def ingest_period(repo: Any, period: Period, fetch: Fetcher) -> str:
    """Fetch, normalize, and commit a single period.

    Args:
        repo: An icechunk repository to write into.
        period: The period to ingest.
        fetch: Callable returning the raw source dataset for the period.

    Returns:
        The snapshot id of the commit.
    """
    ds = normalize(fetch(period))
    ds = ds.chunk(chunking_for(ds))

    session = repo.writable_session("main")
    existing = _has_data(session)
    if existing:
        # align_chunks=True is not optional here. Months are 28-31 days and the
        # store's time chunk is 30, so after a few appends the final zarr chunk
        # is partial and the incoming period straddles it. Without alignment
        # xarray refuses the write outright -- "would overlap multiple Dask
        # chunks" -- because a parallel write across a shared chunk can corrupt
        # it. Alignment rechunks the incoming data to the store's boundaries.
        ds.to_zarr(session.store, append_dim="time", consolidated=False, align_chunks=True)
    else:
        ds.to_zarr(session.store, mode="w", zarr_format=3, consolidated=False)
    return str(session.commit(f"ingest {period.period_id}"))

ingest(repo, periods, fetch, *, resume=True, stop_after=None)

Ingest a list of periods, one commit each.

Parameters:

Name Type Description Default
repo Any

An icechunk repository to write into.

 required
periods list[Period]

Periods to ingest, in chronological order.

 required
fetch Fetcher

Callable returning the raw source dataset for a period.

 required
resume bool

When True, skip periods the store already holds.

 True
stop_after int | None

Stop once this many periods have been ingested, simulating an interrupted run.

 None

Returns:

Name Type Description
An IngestReport

class:IngestReport describing what happened.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/ingest.py 
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def ingest(
    repo: Any,
    periods: list[Period],
    fetch: Fetcher,
    *,
    resume: bool = True,
    stop_after: int | None = None,
) -> IngestReport:
    """Ingest a list of periods, one commit each.

    Args:
        repo: An icechunk repository to write into.
        periods: Periods to ingest, in chronological order.
        fetch: Callable returning the raw source dataset for a period.
        resume: When True, skip periods the store already holds.
        stop_after: Stop once this many periods have been ingested, simulating
            an interrupted run.

    Returns:
        An :class:`IngestReport` describing what happened.
    """
    report = IngestReport()
    present = committed_periods(repo) if resume else set()

    for period in periods:
        if stop_after is not None and len(report.ingested) >= stop_after:
            break
        if period.period_id in present:
            report.skipped.append(period.period_id)
            continue
        try:
            snapshot = ingest_period(repo, period, fetch)
        except Exception as exc:
            report.failed[period.period_id] = f"{type(exc).__name__}: {exc}"
            continue
        report.ingested.append(period.period_id)
        report.snapshots[period.period_id] = snapshot

    return report

store_path(base, dataset_id)

Return the on-disk path for a dataset's store.

Mirrors the open-climate-service layout, {data_dir}/downloads/{dataset_id}.icechunk.

Parameters:

Name Type Description Default
base Path | str

The instance's data directory.

 required
dataset_id str

Public identifier of the dataset.

 required

Returns:

Type Description
Path

The store path.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/ingest.py 
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def store_path(base: Path | str, dataset_id: str) -> Path:
    """Return the on-disk path for a dataset's store.

    Mirrors the open-climate-service layout,
    ``{data_dir}/downloads/{dataset_id}.icechunk``.

    Args:
        base: The instance's data directory.
        dataset_id: Public identifier of the dataset.

    Returns:
        The store path.
    """
    return Path(base) / "downloads" / f"{dataset_id}.icechunk"

indices

Climate indices: the derived products a service actually publishes.

Raw temperature and rainfall are inputs, not answers. What a health ministry or planning office asks for is "how many hot days", "was this month unusually dry", "when does the rainy season start" -- indices computed from the stored series. These are the shape of the processes open-climate-service exposes over openEO, implemented directly here so the arithmetic is visible.

Functions:

climatological_normal(ds, variable='t2m')

Return the per-month mean over all years: the climatological normal.

Parameters:

Name Type Description Default
ds Dataset

A dataset with a daily time axis.

 required
variable str

Which variable to summarize.

 't2m'

Returns:

Type Description
DataArray

A (month, y, x) array of long-run monthly means.

Raises:

Type Description
KeyError

If the variable is absent.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/indices.py 
34
35
36
37
38
39
40
41
42
43
44
45
46
47
def climatological_normal(ds: xr.Dataset, variable: str = "t2m") -> xr.DataArray:
    """Return the per-month mean over all years: the climatological normal.

    Args:
        ds: A dataset with a daily ``time`` axis.
        variable: Which variable to summarize.

    Returns:
        A ``(month, y, x)`` array of long-run monthly means.

    Raises:
        KeyError: If the variable is absent.
    """
    return _require(ds, variable).groupby("time.month").mean()

monthly_anomaly(ds, variable='t2m')

Return each timestep's departure from its month's normal.

Anomalies, not absolute values, are what make two places or two years comparable -- which is why nearly every published climate product is one.

Parameters:

Name Type Description Default
ds Dataset

A dataset with a daily time axis.

 required
variable str

Which variable to anomalize.

 't2m'

Returns:

Type Description
DataArray

An array shaped like the input, in the same units.

Raises:

Type Description
KeyError

If the variable is absent.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/indices.py 
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
def monthly_anomaly(ds: xr.Dataset, variable: str = "t2m") -> xr.DataArray:
    """Return each timestep's departure from its month's normal.

    Anomalies, not absolute values, are what make two places or two years
    comparable -- which is why nearly every published climate product is one.

    Args:
        ds: A dataset with a daily ``time`` axis.
        variable: Which variable to anomalize.

    Returns:
        An array shaped like the input, in the same units.

    Raises:
        KeyError: If the variable is absent.
    """
    values = _require(ds, variable)
    normal = values.groupby("time.month").mean()
    anomaly = values.groupby("time.month") - normal
    anomaly.attrs = {**values.attrs, "long_name": f"{variable} anomaly"}
    return anomaly

hot_days(ds, threshold=30.0, variable='t2m')

Count days per month above a temperature threshold.

Parameters:

Name Type Description Default
ds Dataset

A dataset with daily temperature in degrees Celsius.

 required
threshold float

The temperature above which a day counts as hot.

 30.0
variable str

Which variable to threshold.

 't2m'

Returns:

Type Description
DataArray

A (time, y, x) array of monthly counts, resampled to month ends.

Raises:

Type Description
KeyError

If the variable is absent.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/indices.py 
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def hot_days(ds: xr.Dataset, threshold: float = 30.0, variable: str = "t2m") -> xr.DataArray:
    """Count days per month above a temperature threshold.

    Args:
        ds: A dataset with daily temperature in degrees Celsius.
        threshold: The temperature above which a day counts as hot.
        variable: Which variable to threshold.

    Returns:
        A ``(time, y, x)`` array of monthly counts, resampled to month ends.

    Raises:
        KeyError: If the variable is absent.
    """
    values = _require(ds, variable)
    counts = (values > threshold).resample(time="1ME").sum()
    counts.attrs = {"units": "days", "long_name": f"days above {threshold} degC"}
    return counts

wet_days(ds, threshold=1.0, variable='tp')

Count days per month with rainfall at or above a threshold.

One millimetre is the conventional cutoff for a "wet day": below it, the reading is indistinguishable from dew or gauge noise.

Parameters:

Name Type Description Default
ds Dataset

A dataset with daily precipitation in millimetres.

 required
threshold float

Millimetres at or above which a day counts as wet.

 1.0
variable str

Which variable to threshold.

 'tp'

Returns:

Type Description
DataArray

A (time, y, x) array of monthly counts.

Raises:

Type Description
KeyError

If the variable is absent.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/indices.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def wet_days(ds: xr.Dataset, threshold: float = 1.0, variable: str = "tp") -> xr.DataArray:
    """Count days per month with rainfall at or above a threshold.

    One millimetre is the conventional cutoff for a "wet day": below it, the
    reading is indistinguishable from dew or gauge noise.

    Args:
        ds: A dataset with daily precipitation in millimetres.
        threshold: Millimetres at or above which a day counts as wet.
        variable: Which variable to threshold.

    Returns:
        A ``(time, y, x)`` array of monthly counts.

    Raises:
        KeyError: If the variable is absent.
    """
    values = _require(ds, variable)
    counts = (values >= threshold).resample(time="1ME").sum()
    counts.attrs = {"units": "days", "long_name": f"days with at least {threshold} mm"}
    return counts

monthly_total(ds, variable='tp')

Sum a variable per month -- the right reduction for rainfall.

Temperature is intensive and gets averaged; rainfall is extensive and gets summed. Using the wrong one is a classic and silent error.

Parameters:

Name Type Description Default
ds Dataset

A dataset with a daily time axis.

 required
variable str

Which variable to total.

 'tp'

Returns:

Type Description
DataArray

A (time, y, x) array of monthly totals.

Raises:

Type Description
KeyError

If the variable is absent.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/indices.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def monthly_total(ds: xr.Dataset, variable: str = "tp") -> xr.DataArray:
    """Sum a variable per month -- the right reduction for rainfall.

    Temperature is intensive and gets averaged; rainfall is extensive and gets
    summed. Using the wrong one is a classic and silent error.

    Args:
        ds: A dataset with a daily ``time`` axis.
        variable: Which variable to total.

    Returns:
        A ``(time, y, x)`` array of monthly totals.

    Raises:
        KeyError: If the variable is absent.
    """
    values = _require(ds, variable)
    totals = values.resample(time="1ME").sum()
    totals.attrs = {**values.attrs, "long_name": f"monthly total {variable}"}
    return totals

spi_like(ds, variable='tp')

Standardize monthly rainfall totals against their own month's history.

A simplified standardized precipitation index: for each calendar month, subtract that month's long-run mean and divide by its standard deviation, so -2 means "far drier than this month usually is". The real SPI fits a gamma distribution first; the standardization idea is the same.

Parameters:

Name Type Description Default
ds Dataset

A dataset with daily precipitation.

 required
variable str

Which variable to standardize.

 'tp'

Returns:

Type Description
DataArray

A (time, y, x) array of dimensionless standardized anomalies.

Raises:

Type Description
KeyError

If the variable is absent.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/indices.py 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def spi_like(ds: xr.Dataset, variable: str = "tp") -> xr.DataArray:
    """Standardize monthly rainfall totals against their own month's history.

    A simplified standardized precipitation index: for each calendar month,
    subtract that month's long-run mean and divide by its standard deviation,
    so -2 means "far drier than this month usually is". The real SPI fits a
    gamma distribution first; the standardization idea is the same.

    Args:
        ds: A dataset with daily precipitation.
        variable: Which variable to standardize.

    Returns:
        A ``(time, y, x)`` array of dimensionless standardized anomalies.

    Raises:
        KeyError: If the variable is absent.
    """
    totals = monthly_total(ds, variable)
    grouped = totals.groupby("time.month")
    mean = grouped.mean()
    std = grouped.std()
    # Guard against a month with no variation, which would divide by zero.
    safe_std = std.where(std > 0, np.nan)
    index = (totals.groupby("time.month") - mean).groupby("time.month") / safe_std
    index.attrs = {"units": "1", "long_name": "standardized precipitation index (simplified)"}
    return index

pyramid_levels(ds, levels=3)

Build coarser resolutions by repeated 2x2 mean downsampling.

This is how open-climate-service builds the multiscale GeoZarr pyramid a map viewer needs: level 0 is full resolution, each subsequent level halves both spatial dimensions so a zoomed-out tile reads a small array instead of the whole grid.

Parameters:

Name Type Description Default
ds Dataset

A dataset with y and x dimensions.

 required
levels int

Total number of levels including level 0; must be at least 1.

 3

Returns:

Type Description
list[Dataset]

A list of datasets, coarsest last.

Raises:

Type Description
ValueError

If levels is less than 1.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/indices.py 
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
def pyramid_levels(ds: xr.Dataset, levels: int = 3) -> list[xr.Dataset]:
    """Build coarser resolutions by repeated 2x2 mean downsampling.

    This is how open-climate-service builds the multiscale GeoZarr pyramid a
    map viewer needs: level 0 is full resolution, each subsequent level halves
    both spatial dimensions so a zoomed-out tile reads a small array instead
    of the whole grid.

    Args:
        ds: A dataset with ``y`` and ``x`` dimensions.
        levels: Total number of levels including level 0; must be at least 1.

    Returns:
        A list of datasets, coarsest last.

    Raises:
        ValueError: If levels is less than 1.
    """
    if levels < 1:
        raise ValueError(f"levels must be at least 1, got {levels}")
    out = [ds]
    current = ds
    for _ in range(levels - 1):
        if current.sizes.get("y", 1) < 2 or current.sizes.get("x", 1) < 2:
            break
        # xarray injects the reduction methods onto Coarsen at runtime, so
        # neither type checker can see .mean(); go through Any deliberately.
        coarsened: Any = current.coarsen(y=2, x=2, boundary="trim")
        current = coarsened.mean()
        out.append(current)
    return out

publish

Publishing: GeoZarr attributes and STAC metadata.

Storing the data is not the same as publishing it. A client that finds this store needs to know where on Earth the grid sits, in which CRS, what time range it covers, and what the variables mean. GeoZarr answers the first two with root attributes; STAC answers the rest with a collection document that a catalogue can index.

Functions:

grid_transform(ds)

Return the affine transform placing a north-up grid on Earth.

The six values are [stepX, rotX, originX, rotY, stepY, originY] with the origin on the OUTER EDGE of the first cell -- pixel registration, not cell centres -- and a negative y step for a north-up grid. Getting the half-cell offset wrong shifts every rendered tile by half a pixel.

Parameters:

Name Type Description Default
ds Dataset

A dataset with y and x coordinates.

 required

Returns:

Type Description
list[float]

The six affine coefficients.

Raises:

Type Description
ValueError

If the grid has fewer than two cells on an axis.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/publish.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def grid_transform(ds: xr.Dataset) -> list[float]:
    """Return the affine transform placing a north-up grid on Earth.

    The six values are ``[stepX, rotX, originX, rotY, stepY, originY]`` with
    the origin on the OUTER EDGE of the first cell -- pixel registration, not
    cell centres -- and a negative y step for a north-up grid. Getting the
    half-cell offset wrong shifts every rendered tile by half a pixel.

    Args:
        ds: A dataset with ``y`` and ``x`` coordinates.

    Returns:
        The six affine coefficients.

    Raises:
        ValueError: If the grid has fewer than two cells on an axis.
    """
    if ds.sizes.get("x", 0) < 2 or ds.sizes.get("y", 0) < 2:
        raise ValueError("a transform needs at least two cells on each axis")

    x = ds["x"].values
    y = ds["y"].values
    step_x = float(x[1] - x[0])
    step_y = float(y[1] - y[0])
    origin_x = float(x[0]) - step_x / 2.0
    origin_y = float(y[0]) - step_y / 2.0
    return [step_x, 0.0, origin_x, 0.0, step_y, origin_y]

bounding_box(ds)

Return [west, south, east, north] covering the grid's outer edges.

Parameters:

Name Type Description Default
ds Dataset

A dataset with y and x coordinates.

 required

Returns:

Type Description
list[float]

The bounding box in the stored CRS.

Raises:

Type Description
ValueError

If the grid has fewer than two cells on an axis.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/publish.py 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def bounding_box(ds: xr.Dataset) -> list[float]:
    """Return ``[west, south, east, north]`` covering the grid's outer edges.

    Args:
        ds: A dataset with ``y`` and ``x`` coordinates.

    Returns:
        The bounding box in the stored CRS.

    Raises:
        ValueError: If the grid has fewer than two cells on an axis.
    """
    step_x, _, origin_x, _, step_y, origin_y = grid_transform(ds)
    far_x = origin_x + step_x * ds.sizes["x"]
    far_y = origin_y + step_y * ds.sizes["y"]
    return [min(origin_x, far_x), min(origin_y, far_y), max(origin_x, far_x), max(origin_y, far_y)]

geozarr_attrs(ds)

Build the GeoZarr root attributes for a dataset.

Parameters:

Name Type Description Default
ds Dataset

A normalized dataset with dims (time, y, x).

 required

Returns:

Type Description
dict[str, Any]

The attribute mapping to write at the store root.

Raises:

Type Description
ValueError

If the grid is too small to place.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/publish.py 
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
def geozarr_attrs(ds: xr.Dataset) -> dict[str, Any]:
    """Build the GeoZarr root attributes for a dataset.

    Args:
        ds: A normalized dataset with dims ``(time, y, x)``.

    Returns:
        The attribute mapping to write at the store root.

    Raises:
        ValueError: If the grid is too small to place.
    """
    return {
        "spatial:transform": grid_transform(ds),
        # Array order, y first: read positionally by clients. Naming these
        # x-first transposes every raster that reads the store.
        "spatial:dimensions": ["y", "x"],
        "spatial:shape": [int(ds.sizes["y"]), int(ds.sizes["x"])],
        "spatial:bbox": bounding_box(ds),
        "proj:code": CRS,
        "zarr_conventions": [{"name": "geozarr", "version": "0.4"}],
    }

temporal_extent(ds)

Return the ISO 8601 start and end of a dataset's time axis.

Parameters:

Name Type Description Default
ds Dataset

A dataset with a time coordinate.

 required

Returns:

Type Description
list[str]

A two-element list of ISO timestamps.

Raises:

Type Description
ValueError

If the dataset has no time values.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/publish.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
def temporal_extent(ds: xr.Dataset) -> list[str]:
    """Return the ISO 8601 start and end of a dataset's time axis.

    Args:
        ds: A dataset with a ``time`` coordinate.

    Returns:
        A two-element list of ISO timestamps.

    Raises:
        ValueError: If the dataset has no time values.
    """
    if "time" not in ds.coords or ds.sizes.get("time", 0) == 0:
        raise ValueError("dataset has no time coordinate to describe")
    stamps = pd.DatetimeIndex(np.asarray(ds["time"].values))
    return [stamps[0].isoformat() + "Z", stamps[-1].isoformat() + "Z"]

stac_collection(ds, dataset_id, *, title=None, description='', zarr_href=None, now=None)

Build a STAC Collection document describing a published dataset.

STAC is how a client discovers what an instance holds without knowing anything about its internals: one document per dataset, with spatial and temporal extent, variable summaries, and a link to the actual store.

Parameters:

Name Type Description Default
ds Dataset

The published dataset.

 required
dataset_id str

Stable public identifier, used as the collection id.

 required
title str | None

Human-readable title; defaults to the id.

 None
description str

Longer prose description.

 ''
zarr_href str | None

URL where the store is served, if it is served.

 None
now datetime | None

Timestamp to record as the publication time; defaults to now.

 None

Returns:

Type Description
dict[str, Any]

A STAC Collection as a plain dict, ready to serialize as JSON.

Raises:

Type Description
ValueError

If the dataset lacks the extents STAC requires.
Source code in climate-pipeline/src/ocs_stack_climate_pipeline/publish.py 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def stac_collection(
    ds: xr.Dataset,
    dataset_id: str,
    *,
    title: str | None = None,
    description: str = "",
    zarr_href: str | None = None,
    now: datetime | None = None,
) -> dict[str, Any]:
    """Build a STAC Collection document describing a published dataset.

    STAC is how a client discovers what an instance holds without knowing
    anything about its internals: one document per dataset, with spatial and
    temporal extent, variable summaries, and a link to the actual store.

    Args:
        ds: The published dataset.
        dataset_id: Stable public identifier, used as the collection id.
        title: Human-readable title; defaults to the id.
        description: Longer prose description.
        zarr_href: URL where the store is served, if it is served.
        now: Timestamp to record as the publication time; defaults to now.

    Returns:
        A STAC Collection as a plain dict, ready to serialize as JSON.

    Raises:
        ValueError: If the dataset lacks the extents STAC requires.
    """
    bbox = bounding_box(ds)
    interval = temporal_extent(ds)
    stamp = (now or datetime.now(UTC)).isoformat()

    variables: dict[str, Any] = {}
    for name, var in ds.data_vars.items():
        values = np.asarray(var.values, dtype="float64")
        variables[str(name)] = {
            "units": var.attrs.get("units", "unknown"),
            "long_name": var.attrs.get("long_name", str(name)),
            "min": round(float(np.nanmin(values)), 4),
            "max": round(float(np.nanmax(values)), 4),
        }

    collection: dict[str, Any] = {
        "type": "Collection",
        "stac_version": "1.0.0",
        "id": dataset_id,
        "title": title or dataset_id,
        "description": description,
        "license": "proprietary",
        "extent": {
            "spatial": {"bbox": [bbox]},
            "temporal": {"interval": [interval]},
        },
        "summaries": {
            "variables": variables,
            "proj:code": [CRS],
        },
        "properties": {"published": stamp},
        "links": [],
    }
    if zarr_href:
        collection["assets"] = {
            "zarr": {
                "href": zarr_href,
                "type": "application/vnd+zarr",
                "roles": ["data"],
                "title": "Zarr store",
            }
        }
    return collection