The state catalog

Summary¶

Decision. Model states are indexed by geocatalog, through a new geocatalog.states namespace over the existing backends, and scored by xrtoolz’s existing forecast metrics. There is no statecatalog package, no StateSlice and no second GeoParquet format. The overview is updated to match in the same change. No library starts importing another: pipekit-cycle, geocatalog and xrtoolz stay independent, and scripts join them.

• A. State columns in geocatalog. Today a forecast’s run time and member can live in extra columns, but nothing names, types or fills them. Reserve five columns in schema v1 and add a builder that reads them from CF or cfgrib metadata.

• B. Run-aware loading and typed queries. load_xarray concatenates on time and drops duplicates, so overlapping runs lose steps silently. Add load_states, which stacks files into an (init time, member, lead time) cube, or into an xr.DataTree with one leaf per producer when producers use different grids, and typed query functions that replace hand-written where strings.

• C. Forecast alignment in xrtoolz. SkillByLeadTime and CRPS score a forecast cube today, but the caller builds the verifying reference by hand, and pixel metrics cannot reduce over a dimension the reference lacks. Add align_reference, broadcast the reference in the pixel metrics, and let geo.Reduce weight by cell area, so multi-model area means need no user code.

• D. Ensemble histories in pipekit-cycle. Cycle and DACycle keep their trajectory in history; EnsembleCycle and EnsembleDACycle keep nothing. Add save_history to both. Writing files stays in scripts, as AGENTS.md requires.

End state. A script runs a cycle and writes one file per run and member. geocatalog indexes the files with typed state columns and loads them back as a cube without mixing runs. xrtoolz aligns the cube with its analyses and scores it by lead time. Runs from several models on several grids load as a tree, are summarised on their own grids and are regridded conservatively only for the comparison. pipekit-experiment’s Run.pipeline_hash names the producer. Every one of the scoping report’s queries is either a function in geocatalog.states or an existing call.

The problem, as the scoping report put it¶

An observation has one time. A model state has several, and the scoping report named them well. Let $r$ be the run time (the initialisation time of a forecast), $v$ the valid time (the time the state describes), and

the lead time ((1)). An analysis is the state a data-assimilation step produces at its own valid time, so $\tau = 0$. A forecast run started at $r$ produces states at $\tau = \Delta t, 2 \Delta t, \ldots$ for a model step $\Delta t$. An ensemble repeats each run $N_e$ times, with members $m = 0, \ldots, N_e - 1$. Each state also has a producer, the hash of the model and configuration that made it, and its parents, the states it was computed from: one for a forecast (the state its run started from), and the whole forecast ensemble for an ensemble-Kalman analysis, whose gain is estimated from all members at once. Table 1 lists each field with its CF name and where it lives today.

The questions the scoping report wanted answered are the right ones: which states are valid at a time, what the latest analysis before a time is, which states one run produced, which members share a run, what a state descends from, and which observations an analysis assimilated. The rest of this page answers each one with the code that exists, and adds what is missing.

Ownership rules¶

States are produced by pipekit-cycle, indexed by geocatalog and scored by xrtoolz, and only scripts join the three (Figure 1).

geocatalog owns the index of states (geocatalog.states): the reserved columns, the state builder, run-aware loading and the state queries. It never imports pipekit or xrtoolz.

The other three follow the stack rules: pipekit-cycle produces states in memory and writes no files, xrtoolz verifies them, pipekit-experiment names the producer (Run.pipeline_hash), none of them imports another, and scripts own every side effect.

The running example: a regional methane forecast archive¶

The archive. A regional methane system covers a 4 × 4 grid over the Permian Basin (31° to 34° N, 104° to 101° W). Each day it runs one data-assimilation step, which writes an analysis of the column-methane anomaly $a(v, x)$ in ppb at grid cell $x$. From each of the first three analyses it launches a two-member ensemble forecast with leads of 0, 1 and 2 days. The archive therefore holds 5 analysis files and $R \times N_e = 3 \times 2 = 6$ forecast files, one per run and member, each NetCDF with a time axis of valid times. Each file’s attributes record its kind, run time, member, producer hash and parent file.

The data. The truth is a random walk in time; the analyses add noise with a standard deviation of 0.05 ppb, and the forecast error grows with lead as $\sigma_\tau = 0.1 + 0.5 \tau$ (in days). The data are synthetic, so the scores in stage 4 can be checked against their expected values.

The notebook notebooks/reports/report13_state_catalog.py builds the archive from a fixed seed in a temporary directory, runs every stage, and asserts each number on this page; it runs with geocatalog 0.5.1 and xrtoolz 0.0.5 from references/. Stage 6 widens the example from one regional system to a multi-model comparison, with its own six-file archive in the same notebook. The blocks below are excerpts with their real outputs; names such as cat, t3 and five_days are defined in the notebook. Lines marked user (markers) are glue the workstreams remove.

Stage 1 · Index the archive (scripts → geocatalog)¶

TL;DR. Index the eleven files so that later questions never open a file.

Problem. Each file $f$ has a footprint, a valid-time interval $[v_{\min}(f), v_{\max}(f)]$ and a path, which geocatalog’s xarray kind already records. A state index also needs $r(f)$, $m(f)$, $h(f)$ and the parent $p(f)$ on every row.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# Shapes: R = 3 runs, N_e = 2 members, L = 3 leads, H = W = 4 cells; 11 files
from geocatalog.build import build_xarray_catalog
from geocatalog.storage import from_geoparquet, to_geoparquet

# one row per file: footprint, [v_min, v_max], path
cat = build_xarray_catalog(paths, crs="EPSG:4326")      # 11 files → 11 rows  today
# r(f), m(f), h(f), p(f): read from the file attributes
attrs = [xr.open_dataset(p).attrs for p in cat.gdf.filepath]           # user
cat.gdf["run_time"] = pd.to_datetime([a["run_time"] for a in attrs])   # user
cat.gdf["member"] = [int(a["member"]) for a in attrs]                  # user
cat.gdf["producer"] = [a["producer"] for a in attrs]                   # user
cat.gdf["parent_uri"] = [a["parent_uri"] for a in attrs]               # user
to_geoparquet(cat, "states.parquet")                     # extras kept          today
cat = from_geoparquet("states.parquet")                  # run_time datetime64  today

The extra columns survive the GeoParquet round trip with their types (run_time as datetime64[us], member as int64), because geocatalog keeps every non-reserved column in CatalogRow.extras. What does not exist is a builder that fills them: the xarray builder reads neither forecast_reference_time nor realization, and a file without a time coordinate gets the placeholder interval 1900 to 2100 with n_timesteps=0. Workstream A adds the builder and reserves the names.

Stage 2 · Ask the scoping report’s questions (geocatalog)¶

TL;DR. Find one run’s forecasts, the latest analysis before a time, and the members of a run, using only the index.

Problem. With $\mathcal{C}$ the set of rows, the three questions are the filters (2):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# chain(r₀ = 01-02, m₀ = 0)
cat.where("kind == 'forecast' and member == 0 and run_time == '2026-01-02'")
# → 0 rows: a string equality on a datetime column matches nothing    today
cat.where("kind == 'forecast' and member == 0"
          " and run_time >= '2026-01-02' and run_time < '2026-01-03'")
# → ['fc_20260102_m0.nc']                                               today

# latest(01-03 06Z): end times live in the IntervalIndex after a round trip
ana = cat.where("kind == 'analysis'").gdf                              # today
ana = ana[ana.index.right <= pd.Timestamp("2026-01-03 06:00")]         # user
ana.iloc[ana.index.right.argmax()].filepath   # → 'ana_20260103.nc'    user

# members(r₀ = 01-01)
cat.where("kind == 'forecast' and run_time < '2026-01-02'").gdf.member
# → [0, 1]                                                              today

Every answer is right, but each one is a string the user must get exactly right. The first query shows why that matters: on pandas 3.0.6, the version in pixi.lock, DataFrame.query("t == '2026-01-02'") on a datetime column returns no rows, while df[df.t == "2026-01-02"] returns one. The failure is silent. Typed functions that take a pd.Timestamp (workstream B2) remove the trap.

Stage 3 · Load a forecast cube (geocatalog → xarray)¶

TL;DR. Load the six forecast files as one array indexed by run, member and lead.

Problem. The cube xrtoolz expects is (3),

with a two-dimensional coordinate $v = r + \tau$. Three runs overlap in valid time: all three are valid on 3 January, at leads 2, 1 and 0.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# Shapes: R = 3, N_e = 2, L = 3, H = W = 4
from geocatalog import GeoSlice
from geocatalog.load import load_xarray

rows = cat.where("kind == 'forecast'").gdf                # 6 rows   today
jan3 = GeoSlice(bounds=(-104.5, 30.5, -100.5, 34.5),
                interval=pd.Interval(t3, t3, closed="both"),
                resolution=(1.0, 1.0), crs="EPSG:4326")          # today
m0 = cat.where("kind == 'forecast' and member == 0").query(jan3)  # 3 rows   today
load_xarray(m0, jan3).sizes["time"]                  # → 1: two runs dropped  today

# F[r, m, τ, x]: one file per (r, m), relabelled from v to τ = v − r
pieces = []
for (r, m), row in rows.set_index(["run_time", "member"]).iterrows():  # user
    ds = xr.open_dataset(row.filepath).load()                         # user
    ds = (ds.assign_coords(lead_time=("time", ds.time.values - r))    # user
            .swap_dims(time="lead_time").drop_vars("time"))          # user
    pieces.append(ds.expand_dims(init_time=[r], member=[m]))          # user
fc = xr.combine_by_coords(pieces, combine_attrs="drop")              # user
fc = fc.assign_coords(valid_time=fc.init_time + fc.lead_time)        # user
# → sizes {'init_time': 3, 'member': 2, 'lead_time': 3, 'lat': 4, 'lon': 4}

load_xarray is right for observations, where two files covering the same time are duplicates. For forecasts the same valid time from two runs is two different states, and drop_duplicates on time keeps one of them without a warning. The eight user lines are the cube builder that workstream B1 puts in load_states.

Stage 4 · Score the forecasts by lead time (xrtoolz)¶

TL;DR. Compare each forecast with the analysis valid at the same time, and report the error at each lead.

Problem. The verifying reference for run $r$ at lead $\tau$ is the analysis at its valid time, $A[r, \tau, x] = a(r + \tau, x)$. With $\bar{F}$ the ensemble mean over members and $N = R H W$ the number of scored values per lead, the two scores are (4),

with the CRPS of an ensemble as in the overview ((12)). For this archive the expected ensemble-mean RMSE is $(\sigma_\tau^2 / 2 + 0.05^2)^{1/2}$, which is 0.087, 0.427 and 0.779 ppb at leads 0, 1 and 2 days.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# Shapes: F (R, N_e, L, H, W); A (R, L, H, W)
from xrtoolz.metrics import CRPS, RMSE, SkillByLeadTime

an = load_xarray(cat.where("kind == 'analysis'"), five_days)  # (T = 5, H, W)  today
# A[r, τ, x] = a(r + τ, x)
ref = an.sel(time=fc.valid_time).drop_vars("time")       # → (R, L, H, W)      user
dims = ("init_time", "lat", "lon")
# RMSE(τ) of the ensemble mean
SkillByLeadTime(RMSE("xch4", dims=dims))(fc.mean("member"), ref)  # → (L,)   today
# → [0.076, 0.467, 0.746]
# CRPS(τ) over members
SkillByLeadTime(CRPS("xch4", ensemble_dim="member", dims=dims))(fc, ref)  # today
# → [0.048, 0.356, 0.605]
# per-member RMSE: reduce over the member dimension too
SkillByLeadTime(RMSE("xch4", dims=("member", *dims)))(fc, ref)           # today
# → ValueError: tuple.index(x): x not in tuple
SkillByLeadTime(RMSE("xch4", dims=("member", *dims)))(*xr.broadcast(fc, ref))
# → [0.109, 0.617, 1.132]                                               user

The measured RMSE (0.076, 0.467 and 0.746) sits within sampling error of the expected values, with $N = 48$ cells per lead. The per-member RMSE matches its expected $(\sigma_\tau^2 + 0.05^2)^{1/2}$, which is 0.112, 0.602 and 1.101. The scoring is shipped. What is missing is the alignment line (workstream C1) and the broadcast that lets pixel metrics reduce over a dimension only the forecast has (C2).

Stage 5 · Trace lineage and the assimilation window (geocatalog)¶

TL;DR. Follow a state back through the runs it came from, and find the observations an analysis used.

Problem. With $P(s)$ the set of parents of state $s$, the lineage of $s$ is every state reachable by following parents, $s$, $P(s)$, $P(P(s))$, and so on, until no state has parents left: a directed acyclic graph, which is a chain when every state has one parent. This archive is the simple case: each analysis lists one parent, the member-0 forecast. An ensemble-Kalman analysis would list every forecast member, because its gain is computed from all of them. For an analysis at time $t$ with an assimilation window of width $w$, the observations used are the rows $o$ of an observation catalog whose valid time falls in $[t - w/2, t + w/2]$.

1
2
3
4
5
6
7
8
9
10
11
12
13
# s, P(s), P(P(s)), … until P = ∅; one parent per state in this archive
by_uri = cat.gdf.set_index("filepath")                                  # user
uri, chain = "ana_20260104.nc", []                                      # user
while uri:                                                              # user
    chain.append(uri)                                                   # user
    uri = by_uri.loc[uri, "parent_uri"]                                 # user
# → ana_20260104 ← fc_20260103_m0 ← ana_20260103 ← fc_20260102_m0
#   ← ana_20260102 ← fc_20260101_m0 ← ana_20260101

# obs with v ∈ [t − w/2, t + w/2]                              # arguments sketched
window = GeoSlice(bounds=aoi, interval=pd.Interval(t - w / 2, t + w / 2,
                  closed="both"), resolution=res, crs="EPSG:4326")  # today
obs_used = obs_cat.query(window)                       # (M rows)          today

The lineage walk is five lines of user code when each state has one parent, and a graph search when states have several. The scoping report’s assimilation_window needs no new code: it is GeoCatalog.query with the window as the interval, or matchup with WithinWindow(start, end) when many analyses are matched at once. Workstream B2 adds lineage; the window stays an existing call.

Stage 6 · Many producers on many grids: a CMIP-style comparison (geocatalog → xrtoolz)¶

TL;DR. Compare models that run on different grids, either by computing each model’s statistics on its own grid or by regridding all of them to one grid first, and see when the two give different answers.

Problem. A multi-model ensemble such as CMIP6 Eyring et al., 2016 holds runs from many models $k \in \mathcal{K}$. In CMIP6 terms the model is the source_id, a run is a member_id (a variant label such as r1i1p1f1), and each model writes on its own grid (grid_label). Model $k$ has $N_k$ members and a grid of cells $i$ with areas $a_i$; the field is $x_k(t, i)$ at day $t$. The archive here is synthetic: model A on a 1° grid with 3 members, B on 2.5° with 1, C on 2° with 2, 90 days each, with independent cell-scale weather of standard deviation 2 K on a zonal climatology. The six files are indexed like stage 1, with producer and member columns.

The area mean of a field $y$ on a grid with cells $i$ is $\langle y \rangle = \sum_i a_i y_i / \sum_i a_i$. A conservative regrid $\mathcal{R}$ onto a common grid with cells $j$ and areas $A_j$ weights each source cell by its overlap $w_{ji} = |c_j \cap c_i|$, so that the area mean survives whenever the common grid covers the source grid:

A statistic $S$ over time and members, applied cell by cell, either commutes with $\mathcal{R}$ or does not. The time mean is linear, so $S(\mathcal{R} x) = \mathcal{R} S(x)$. The 99th percentile $q_{0.99}$ is not: averaging cells before taking the percentile removes the cell-scale extremes that the percentile measures. That gives the three pipelines of (6):

where (3) equals (1) by (5), and (2) equals (1) only when $S$ is linear. Pipeline (3) is the usual way to draw a multi-model map: each statistic is computed where the model resolved it, and only the result is regridded. Combining the models needs one more choice: one vote per model, $\bar{x} = |\mathcal{K}|^{-1} \sum_k \langle S(x_k) \rangle$, or one vote per member, $\sum_k N_k \langle S(x_k) \rangle / \sum_k N_k$.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# Shapes: K = 3 models; model k: (N_k, T = 90, H_k, W_k), H_k × W_k its own grid
from xrtoolz.interpolate import RegridConservative

# one leaf per producer, members stacked inside it; the grids differ
leaves = {}
for k, rows in cat.gdf.groupby("producer"):                             # user
    runs = [xr.open_dataset(p).load() for p in rows.filepath]           # user
    leaves[f"/{k}"] = xr.concat(runs, dim="member")                     # user
tree = xr.DataTree.from_dict(leaves)                                    # today
# → A (3, 90, 180, 360) · B (1, 90, 72, 144) · C (2, 90, 90, 180)

# 𝓡: conservative regrid of every leaf onto the coarsest grid (2.5°)
regrid = RegridConservative(grid_25, mode="mean", normalize="destarea")  # today
common = regrid(tree)                       # every leaf → (N_k, 90, 72, 144)  today

# ⟨y⟩ = Σᵢ aᵢ yᵢ / Σᵢ aᵢ, with aᵢ ∝ sin φᵢ₊ − sin φᵢ₋
def area_mean(y): ...                                                   # user

# S = time mean (linear): (1) and (2) agree
# → native    A 277.9999 · B 279.4971 · C 277.2991 K
# → regridded A 277.9999 · B 279.4971 · C 277.2991 K   (max difference 1.2e-11)

# S = q₀.₉₉ − time mean (nonlinear), per cell, then ⟨·⟩
# → (1) native        A 4.289 · B 4.279 · C 4.288 K
# → (2) regrid first  A 1.544 · B 4.279 · C 2.564 K
# → (3) hybrid        A 4.289 · B 4.279 · C 4.288 K    (equals (1) to 1.2e-11)

# one vote per model vs one vote per member, S = time mean
# → 278.2654 K vs 278.0159 K

The tree, the regrid and the map over models all ship today: xrcore’s Operator maps any xrtoolz operator over the leaves of an xr.DataTree, and RegridConservative keeps area means to 1.2e-11 K across a 1° → 2.5° and a non-nested 2° → 2.5° regrid. The numbers show why the order of operations matters. Regridding first leaves the time mean untouched but removes 64 % of model A’s warm extremes and 40 % of model C’s, so the fine models look less extreme than the coarse one only because they were averaged. Model B is unchanged because it is already on the common grid. Two things are user code: grouping the catalog into leaves, which workstream B1 adds as load_states(..., group_by="producer"), and the area mean, which geo.Reduce cannot weight (workstream C3). The vote is a choice the script must state: here it moves the ensemble mean by 0.25 K.

Comparing each model with observations runs the other way: the observations are regridded onto each model’s native grid, so that the statistic is computed where the model resolved it. RegridConservative and RegridLike take their target grid when they are built, so today that is one operator per leaf, built in a loop (regrid_conservative(obs, leaf) inside xr.map_over_datasets does the same). This is also the step where a pipekit Graph earns its place: two inputs, the model leaf and the observations, feed one bias operator.

The scoping report’s names, checked¶

Table 2 checks every stack name the scoping report’s examples use.

Current-state audit¶

Eleven pieces of the proposed package were checked: seven are healthy as they ship, three are split between a library that has them and a gap, and one is missing (Table 3); Figure 2 counts the page’s calls by marker. The audit read the revisions pinned in references/ on 13 October 2026.

Smaller findings:

• A plain Cycle does not advance time. _ForwardModelAdapter passes the carry state through unchanged, so DAState.t stays 0.0. Only the DA cycles call _advance_da_state. A history entry’s valid time must be computed from its index: entry $k$ with stride $s$ is valid at $r + (ks + 1) \Delta t$ (checked: a 4-step Cycle with stride 2 keeps the states after steps 1 and 3).

• EnsembleCycle threads one carry state through all members. Its loop calls step_op(member_j, state) member by member, so a step that counts calls cannot tell members apart. This is one reason to write from the history (workstream D) rather than from inside the step.

• DAState.t is a float with no epoch. There is no run time or calendar anywhere in pipekit-cycle; the script that writes files supplies $r$.

• CatalogDataset is documented but not built. pipekit-train’s docs/design/api/datasets.md already says its catalog can be “A geocatalog.GeoCatalog (or StateCatalog) reference”. With this redesign the two are the same object.

• xrreader has no forecast products. Its CDS presets are ERA5 reanalyses (product_type="reanalysis") and its CMEMS presets are reanalyses and NRT observations. _TIME_COORDS in cds/source.py includes forecast_reference_time only to trim on time. The download cache maps a request to a path; it keeps no index.

• pipekit-evaluate has the level but not the scoring. LEVEL_NAMES["L4"] is "reanalysis_forecast", and references come from Task.datasets(). No lead-time or CRPS scoring ships there; xrtoolz’s does.

Options considered¶

Table 4 weighs four ways to index model states.

Recommended: option 2. Alternative: option 3 for a single producer with one store, which needs only workstream C.

Workstream A · State columns in geocatalog¶

Decision. The state fields are reserved row columns of geocatalog’s schema, not fields of a new slice (Table 5).

A1 · Reserve the columns¶

1. Register the migration. Bump SCHEMA_VERSION_CURRENT to 1 and register a v0 → v1 function in parquet._MIGRATIONS that adds the four columns as nulls. A v0 catalog then reads as a catalog of states with no state fields.

2. Add the names to the reserved set. CatalogRow.extras stops carrying them, and a CatalogRow.state property returns them typed.

3. Document them in docs/catalog/schema-versions.md, the place the schema module points to.

Lead time is not stored: a row covers $[v_{\min}, v_{\max}]$, so its leads are $[v_{\min} - r, v_{\max} - r]$, derived when needed.

A2 · Build the columns from metadata¶

1
2
3
4
5
6
7
8
9
10
11
12
# Shapes: 11 files → 11 rows
# Today:
cat = build_xarray_catalog(paths, crs="EPSG:4326")                     # today
attrs = [xr.open_dataset(p).attrs for p in cat.gdf.filepath]           # user
cat.gdf["run_time"] = pd.to_datetime([a["run_time"] for a in attrs])   # user
# … member, producer, parents the same way                            # user

# After:
from geocatalog.build import build_state_catalog
from geocatalog.states import CF_NAMES
cat = build_state_catalog(paths, crs="EPSG:4326", names=CF_NAMES)  # proposed (A2)
cat.gdf[["kind", "run_time", "member", "producer", "parents"]]    # proposed (A1)

names maps each field to where it is found: a coordinate, a variable, a global attribute or a callable on the path (for archives that encode the run in the file name, such as fc_20260102_m0.nc). kind, producer and parents have no CF name, so they come from attributes the script writes, or, for kind, from a callable on the path. The member column is a string, not the integer of stage 1’s glue: a CF realization is a number, but a CMIP member is a variant label such as r1i1p1f1 whose four indices (realization, initialization, physics, forcing) all distinguish runs, so an integer column would lose them. Stage 1’s members become "0" and "1", and its analyses, written with member=-1, get a null.

Workstream B · Run-aware loading and typed queries¶

Decision. geocatalog.states holds the loader and the queries that know what a run is; the observation loader stays as it is (Table 6); Table 7 places every new name of workstreams A to C.

B1 · Load without mixing runs¶

load_states is stage 3’s eight user lines, with a check that each (run, member) pair appears once. Members missing from a run come back as NaN, not as a shorter cube. With group_by="producer" it is stage 6’s grouping loop instead: one leaf per producer, each on its own grid, so that every xrtoolz operator maps over the models unchanged. The return type follows the argument, never the data: a plain call on rows with different grids raises rather than returning a tree by surprise.

B2 · Typed queries¶

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Shapes: R = 3, N_e = 2, L = 3, H = W = 4
import geocatalog.states as gs

# Today:
cat.where("kind == 'forecast' and member == 0"
          " and run_time >= '2026-01-02' and run_time < '2026-01-03'")  # today

# After:
gs.forecast_chain(cat, run_time=pd.Timestamp("2026-01-02"), member="0")
# → 1 row, fc_20260102_m0.nc                                       proposed (B2)
gs.latest_analysis(cat, before=pd.Timestamp("2026-01-03 06:00"))
# → ana_20260103.nc                                               proposed (B2)
gs.ensemble_at(cat, run_time=pd.Timestamp("2026-01-01"))  # 2 rows  proposed (B2)
fc = gs.load_states(gs.forecasts(cat), aoi)    # → (3, 2, 3, 4, 4)  proposed (B1)
gs.lineage(cat, "ana_20260104.nc")             # → 7 rows           proposed (B2)

latest_analysis filters on the reserved kind column. Zero leads are not enough: in a one-row-per-step archive (a GRIB file per step, or open question 1’s explode="step") a lead-0 forecast row also has $v_{\min} = v_{\max} = r$, so a lead test would return it as an analysis. The member filter takes a string, as the column does (Section ??); an integer raises TypeError rather than silently matching nothing.

Workstream C · Forecast alignment in xrtoolz¶

Decision. xrtoolz builds the verifying reference and broadcasts it; the metrics stay as they are (Table 8).

1
2
3
4
5
6
7
8
9
10
11
12
# Shapes: F (R, N_e, L, H, W); a (T, H, W)
# Today:
ref = an.sel(time=fc.valid_time).drop_vars("time")    # (R, L, H, W)     user
SkillByLeadTime(RMSE("xch4", dims=("member", *dims)))(*xr.broadcast(fc, ref))  # user

# After:
from xrtoolz.metrics import align_reference
# A[r, τ, x] = a(r + τ, x)
ref = align_reference(fc, an, init_dim="init_time", lead_dim="lead_time",
                      time_dim="time")               # → (R, L, H, W) proposed (C1)
SkillByLeadTime(RMSE("xch4", dims=("member", *dims)))(fc, ref)  # (L,)  proposed (C2)
# → [0.109, 0.617, 1.132]

C3 · Area-weighted reductions¶

1
2
3
4
5
6
7
8
# Shapes: leaf k: (N_k, T, H_k, W_k) → (N_k, T)
# ⟨y⟩ = Σᵢ aᵢ yᵢ / Σᵢ aᵢ, aᵢ ∝ Δλ (sin φᵢ₊ − sin φᵢ₋)
# Today:
geo.Reduce("mean", dim=("lat", "lon"))(tree)   # unweighted: wrong on a sphere  today
def area_mean(y): ...                          # sin-latitude weights            user

# After:
geo.Reduce("mean", dim=("lat", "lon"), weights="cell_area")(tree)  # proposed (C3)

Because Reduce is an xrcore Operator, the weighted form maps over the tree like every other operator, and each leaf is weighted by its own cell areas.

Workstream D · Ensemble histories in pipekit-cycle¶

Decision. The cycle’s history is the write seam; files are written by the script, after the call (Table 9).

1
2
3
4
5
6
7
8
9
10
11
12
13
# Shapes: T steps, N_e members, H × W cells; Δt the model step
fcst = pc.EnsembleCycle(step_op=model, n_steps=T, n_members=N_e,
                        save_history=True)                       # proposed (D1)
members, _ = fcst([x_a] * N_e, pc.DAState())                     # today
# v_k = r + (k + 1) Δt for history entry k (stride 1)
for m in range(N_e):                                             # user (script)
    # history entry h = (members, state), so h[0][m] is member m's carrier
    ds = to_dataset([h[0][m] for h in fcst.history], r, dt)      # arguments sketched
    ds.attrs |= {"forecast_reference_time": str(r), "realization": m,
                 "producer": h_hash, "parents": [x_a_uri]}      # user (script)
    ds.to_zarr(root / f"fc_{r:%Y%m%dT%H}_m{m}.zarr")             # user (script)
cat = build_state_catalog(sorted(root.glob("fc_*.zarr")),
                          crs=crs, names=CF_NAMES)               # proposed (A2)

to_dataset is the model’s own carrier-to-xarray step, which differs per model; it is not a library function. h_hash is the producer: sha256_hex(stable_json(fcst.get_config())) today, or Report 9’s proposed pipeline_hash.

Testing gates¶

Each workstream merges only when its gates pass (Table 10). The fixtures are the running example’s eleven files.

Sequencing and risks¶

The three phase-1 pieces are independent, so they ship side by side; the builder, loader and queries follow because they use A1’s columns (Figure 3).

• The first schema migration. A1 is geocatalog’s first registered migration, and a mistake in it touches every saved catalog. Fallback: ship A2 and B with the four names as documented extra columns and no schema bump, and reserve them once the API has settled.

• geocatalog keeps growing. Report 6 found it had already outgrown its scope. Fallback: if geocatalog.states passes about 500 lines, move it to its own distribution in the geotoolz monorepo that depends on geocatalog, still with one index format.

• Metadata conventions vary. CF files, cfgrib-decoded GRIB, and archives that encode the run in the file name all store $r$ and $m$ differently. Fallback: names accepts a callable on the path, so any archive can be indexed without changing the builder.

• One row per file is coarse. A file per (run, member) answers “valid at $t$” at file granularity, and a per-step query must open the file. Fallback: when steps are separate files (common for GRIB), rows are per step without any change; a per-step explode option is open question 1.

• Histories held in memory. A 72-step, 40-member trajectory at 0.25° does not fit in memory. Fallback: run the cycle in chunks of steps from the last state and write after each chunk; a streaming hook on the cycles is open question 3.

• Pandas query semantics. Users who keep writing where strings still hit stage 2’s silent miss. Fallback: B2’s functions use boolean masks, and the where docstring gets a warning about datetime equality.

• Curvilinear model grids. regrid_conservative works along two rectilinear dimensions, and many CMIP ocean models use curvilinear or tripolar native grids (grid_label="gn"). Fallback: use the regridded products the modelling centres publish (grid_label="gr"), or regrid those models with an external tool such as xESMF in the script; open question 5.

Open questions¶

1. One row per file, or per step? Per file keeps catalogs small (40 rows for a 40-member run instead of the scoping report’s 2,880). Per step answers “valid at $t$” without opening files. Should build_state_catalog offer an explode="step" option?

2. Which kinds? kind is reserved with three values: analysis, forecast, reanalysis. Do hindcasts and reforecasts need their own values, or are they forecasts with a producer that says so?

3. A streaming hook on the cycles. Should Cycle and EnsembleCycle take an on_step(k, carriers, state) callable, so scripts can write while the cycle runs? It would also give SmootherCycle and WindowedCycle a history.

4. Recording assimilated observations. The scoping report’s states_that_used(obs_slice) needs each analysis to record which observation rows it used. A CatalogBundle matchup table between the analysis and observation catalogs would do it; is that the right home?

5. Curvilinear grids. Should xrtoolz grow a conservative regrid for curvilinear grids (cell polygons rather than two 1-D weight matrices), or should the stack defer to xESMF for them?

6. Targets that differ per leaf. RegridConservative and RegridLike fix their target grid when built, so regridding one observation field onto each model’s own grid needs one operator per leaf. Should they accept the target as a second input, so the step maps over a tree of targets?

Upcoming work¶

Table 11 lists the follow-ups.

References¶

Sources¶

• The scoping report, “Report 13 — statecatalog: catalog for model states” (supplied by the author).

• The pinned code in references/: geocatalog (geotoolz-catalog), xrtoolz and xrreader, pipekit-cycle, pipekit-train, pipekit-experiment and pipekit-evaluate.

• CF conventions, standard name table (forecast_reference_time, forecast_period, realization).

• The running example’s script and outputs, described in Section ??.