Dynamical Data API

Dynamical Data Package

Download & process numerical weather predictions from Dynamical.org.

We convert the ECMWF ENS 0.25 degree data to these H3 resolution 5 hexagons:

Note: The generic geospatial logic for mapping latitude/longitude grids to H3 hexagons has been extracted to the packages/geo package. This package (dynamical_data) focuses specifically on the ingestion, processing, and storage of time-varying NWP datasets like ECMWF. The H3 grid weights are provided as a Dagster asset from the geo package, eliminating the need for precomputed static files.

Data storage experiments

The storage format itself lives in delta_store.nwp (writer properties, sort order, precision); this section records the measurements behind it. Full before/after detail is in PR #271; earlier experiments (UInt8/Int16 affine quantisation, codec and sort-order sweeps) are in this file's git history.

Current scheme: physical-unit Float32, every continuous variable rounded to a 13-bit significand (max relative error 2⁻¹³ ≈ 1.2×10⁻⁴ — measured ≤ 0.004 °C for temperature, ≤ 8 Pa for MSL pressure), rows sorted init_time → ensemble_member → valid_time → h3_index, plain ZSTD level 3.

How much space does GB-wide ECMWF ENS take? One daily run (1,671 H3 cells × 51 members × 85 lead times, up to ~7.24M rows) averages ~113 MB, so a year is ~41 GB. The full local development table — 810 daily runs (Apr 2024 → Jun 2026, 1.57 billion rows) — is 86 GB.

Storage (9 real partitions spread across every season):

Config	avg MB/partition	extrapolated GB/yr
Previous: Int16 12-bit quantisation, ZSTD-14	115.3	42.1
Float32 + 13-bit significand, ZSTD-3	110.6	40.4
Adopted: same + member-early sort	112.9	41.2
Same + BYTE_STREAM_SPLIT	133.8	48.8

BYTE_STREAM_SPLIT makes this table worse (unlike power_forecasts, where it wins): significand rounding collapses NWP values into repeats that parquet's default dictionary+RLE encoding captures directly, and BYTE_STREAM_SPLIT scatters that repetition across four byte planes. Writer properties are data-dependent — measure per table.

Read path — the member-early sort means each approx 1M-row parquet row group spans only a few ensemble members, so a single-member read (every training run reads just the control member) skips most row groups via min/max stats. Measured on a real 29-day, 9-cell, control-member collect: ~5× faster, ~5× less peak memory (0.15 s / ~1 GB → 0.02–0.04 s / approx 205 MB), for a ~2% storage cost.

Per-variable keep_bits: considered and rejected (2026-07). Since Dynamical's upstream precision caps the real information at 7–12 significand bits per variable, budgets matched to upstream would compress better than the uniform 13. Measured on 4 seasonal partitions through the exact production write path (error columns = error added on top of today's stored values; wind speed's power impact is ~3× its relative error because the speed-to-power curve is roughly cubic):

Config	Size vs today	GB/yr	wind speed rel err	→ power (×3)	temp	MSL
Uniform 13 (today)	—	41.1	0	0	0	0
Upstream-matched (temp 8, wind 7, pressure 12, flux 8)	−19.4%	33.1	0.76%	2.3%	0.06 °C	16 Pa
Uniform 10	−15.6%	34.7	0.10%	0.29%	0.016 °C	64 Pa
Wind-protected (wind speed stays 13, rest squeezed)	−13.0%	35.7	0	0	0.06 °C	16 Pa

The full squeeze adds ~2.3% power-equivalent wind error — not tolerable. The wind-safe ceiling is −13% ≈ 5.4 GB/yr, and the NWP table is GB-wide so it does not grow with the V2 scale-up to ~2,500 time series: a fixed ~5 GB/yr saving doesn't justify maintaining a dict of per-variable precision budgets. If disk ever becomes a real constraint, the wind-protected config is the one to reach for.

Why round at all, when Dynamical.org already rounds? Dynamical stores ECMWF ENS with 6–11 mantissa bits per variable (their binary_rounding.py). But trailing zeros do not survive arithmetic: our H3 aggregation is a weighted mean over grid points, and wind speed/direction are derived from their u/v via sqrt/arctan2, so by the time values reach our writer their mantissas are full entropy again (measured: 100% of Dynamical-style rounded values have zeroed low bits; after a weighted mean, 0.07% do). Our 13-significand-bit rounding restores compressibility while being 1–6 bits finer than the upstream precision, so it discards almost nothing beyond what Dynamical already dropped.

dynamical_data.ecmwf_ens.download

Opening and downloading one ECMWF ENS run from the Dynamical.org catalog.

Attributes

ECMWF_ENS_INSTANTANEOUS_VARS = frozenset(_ECMWF_ENS_VARS_TO_DOWNLOAD) - Nwp.deaccumulated_var_names - Nwp.categorical_var_names module-attribute

The downloaded variables describing conditions at one instant, under their download names.

None of these is ever legitimately null, anywhere in a run, which is why :func:dynamical_data.ecmwf_ens.upstream_nulls.assess_upstream_grid_point_nulls counts them separately from the de-accumulated ones and the ecmwf_ens asset gates a check on that count being zero.

Derived from the download list rather than from Nwp's fields, because the two namespaces differ: we download wind_u_10m/wind_v_10m (and the 100 m pair), and :func:dynamical_data.ecmwf_ens.convert_to_polars.convert_nwp_xarray_dataset_to_polars_dataframe derives wind_speed_*/wind_direction_* from them. A set taken from the contract would name four variables the downloaded dataset does not carry, and indexing it would raise KeyError.

Classes

NwpRunNotYetAvailable

Bases: Exception

Raised when nwp_init_time isn't in the catalog yet (Dynamical hasn't published it).

Source code in packages/dynamical_data/src/dynamical_data/ecmwf_ens/download.py

class NwpRunNotYetAvailable(Exception):
    """Raised when ``nwp_init_time`` isn't in the catalog yet (Dynamical hasn't published it)."""

Functions:

open_ecmwf_ens_run(nwp_init_time, h3_grid)

Lazily open the ECMWF ENS Icechunk store and slice it to the requested run and H3 grid.

No data is downloaded: the returned dataset is still backed by lazy Dask/Zarr arrays. Call :func:download_ecmwf_ens_data to actually fetch the data.

Parameters:

Name	Type	Description	Default
nwp_init_time	datetime	The initialization time to open. Must be timezone aware.	required
h3_grid	DataFrame[H3GridWeights]	The H3 grid to use for spatial bounds.	required

Source code in packages/dynamical_data/src/dynamical_data/ecmwf_ens/download.py

def open_ecmwf_ens_run(
    nwp_init_time: datetime,
    h3_grid: pt.DataFrame[H3GridWeights],
) -> xr.Dataset:
    """Lazily open the ECMWF ENS Icechunk store and slice it to the requested run and H3 grid.

    No data is downloaded: the returned dataset is still backed by lazy Dask/Zarr arrays.
    Call :func:`download_ecmwf_ens_data` to actually fetch the data.

    Args:
        nwp_init_time: The initialization time to open. Must be timezone aware.
        h3_grid: The H3 grid to use for spatial bounds.
    """
    # Convention-sensitive to the *real* Dynamical.org catalog: this function bakes in assumptions
    # about its shape (longitude in [-180, 180], descending latitude, coordinate/dimension names).
    # The offline tests share those assumptions and cannot catch a mismatch with the live
    # catalog, so after changing this function run the network-gated test manually:
    #     uv run pytest --run-network -m network
    # See
    # <https://openclimatefix.github.io/nged-substation-forecast/architecture/testing/#network-gated-tests>.

    # Reusable-package input validation, not a reachable production state: the `ecmwf_ens` asset
    # always sources `h3_grid` from `h3_grid_weights`, which raises on an empty cell list before
    # writing anything, so the file it reads can never hold zero rows.
    if h3_grid.is_empty():
        raise ValueError("h3_grid is empty. Cannot download ECMWF data for an empty grid.")

    if nwp_init_time.utcoffset() is None:
        raise ValueError(f"nwp_init_time must be timezone aware. {nwp_init_time.tzinfo=}")

    # We need to make nwp_init_time tz-naive for the xarray selection.
    utc_nwp_init_time = np.datetime64(nwp_init_time.astimezone(UTC).replace(tzinfo=None))

    ds = dynamical_catalog.open("ecmwf-ifs-ens-forecast-15-day-0-25-degree", chunks=None)

    ds = ds[list(_ECMWF_ENS_VARS_TO_DOWNLOAD)]

    if utc_nwp_init_time not in ds.init_time.values:
        raise NwpRunNotYetAvailable(f"{utc_nwp_init_time} is not in ds.init_time.values")

    # This guards the Dynamical.org catalog itself, which is an external substrate we neither
    # control nor version-pin, so its shape can change under us between runs.
    if ds.longitude.size == 0 or ds.latitude.size == 0:
        raise ValueError("Dataset has empty longitude or latitude coordinates.")

    # Validate longitude range.
    # NOTE: Dynamical.org converts the longitude range to [-180, 180].
    if ds.longitude.min() < -180 or ds.longitude.max() > 180:
        raise ValueError("Dataset longitude must be in the range [-180, 180]")

    min_lat, max_lat, min_lon, max_lon = h3_grid.select(
        min_lat=pl.col("nwp_lat").min(),
        max_lat=pl.col("nwp_lat").max(),
        min_lon=pl.col("nwp_lon").min(),
        max_lon=pl.col("nwp_lon").max(),
    ).row(0)

    lat_slice = _calc_slice_for_lat_or_lng("latitude", ds, min_lat, max_lat)
    lon_slice = _calc_slice_for_lat_or_lng("longitude", ds, min_lon, max_lon)

    # NOTE: This will fail if the region crosses the anti-meridian. But we do not anticipate
    # forecasting near the anti-meridian.
    ds_sliced = ds.sel(latitude=lat_slice, longitude=lon_slice, init_time=utc_nwp_init_time)

    # Explicitly check for an empty spatial intersection after slicing.
    # This prevents downstream KeyErrors during DataFrame conversion.
    if ds_sliced.longitude.size == 0 or ds_sliced.latitude.size == 0:
        raise ValueError("No spatial overlap found between H3 grid and NWP dataset.")

    return ds_sliced

download_ecmwf_ens_data(ds_sliced)

Download (compute) a lazily-opened, already-sliced ECMWF ENS dataset.

Parameters:

Name	Type	Description	Default
ds_sliced	Dataset	A lazy dataset as returned by :func:open_ecmwf_ens_run.	required

Source code in packages/dynamical_data/src/dynamical_data/ecmwf_ens/download.py

def download_ecmwf_ens_data(ds_sliced: xr.Dataset) -> xr.Dataset:
    """Download (compute) a lazily-opened, already-sliced ECMWF ENS dataset.

    Args:
        ds_sliced: A lazy dataset as returned by :func:`open_ecmwf_ens_run`.
    """

    def download_array(var_name: str) -> dict[str, xr.DataArray]:
        return {var_name: ds_sliced[var_name].compute()}

    # The download is I/O bound (S3 network requests). We use a ThreadPoolExecutor to parallelize
    # network latency across multiple variables. A ProcessPoolExecutor would be less efficient here
    # due to the high serialization overhead of Xarray objects between processes.
    #
    # max_workers is capped rather than left at the default (one thread per variable, i.e. 13).
    # Investigation of issue #276 found that 13 concurrent chunked-zarr fetches self-contend badly
    # (S3 rate limiting or connection-pool starvation): most variables finish in 5-20s, but a few
    # straggle for minutes, making the whole download 600s+. Capping at 4 removed the stragglers
    # entirely and cut a real download from 645s to 22.5s.
    #
    # This is a recent regression, not a pre-existing property of the download: Dagster's run
    # history shows per-partition downloads holding a steady ~48-54s right up to 2026-06-30
    # 12:26 UTC, then every run afterwards (2026-07-01 onwards) taking 3-12 min. That boundary
    # lines up exactly with an `icechunk` 2.0.6 -> 2.1.0 bump in the same `uv.lock` update
    # (commit b46d145, 2026-06-30 12:26:50 UTC) — the leading theory is a change in icechunk's
    # underlying S3 client (connection pooling/concurrency handling) between those two versions,
    # though this hasn't been confirmed by pinning back to 2.0.6 and re-testing.
    data_arrays: dict[str, xr.DataArray] = {}
    with concurrent.futures.ThreadPoolExecutor(max_workers=4) as executor:
        futures = [executor.submit(download_array, str(name)) for name in ds_sliced.data_vars]
        for future in concurrent.futures.as_completed(futures):
            data_arrays.update(future.result())

    return xr.Dataset(data_arrays)

dynamical_data.ecmwf_ens.convert_to_polars

Converting a downloaded ECMWF ENS xarray dataset into the Nwp Polars contract.

The regular lat/lon grid is aggregated onto H3 cells on the way.

Attributes

Classes

Functions:

convert_nwp_xarray_dataset_to_polars_dataframe(ds, h3_grid)

Vectorized processing of ECMWF dataset to H3 grid.

Source code in packages/dynamical_data/src/dynamical_data/ecmwf_ens/convert_to_polars.py

def convert_nwp_xarray_dataset_to_polars_dataframe(
    ds: xr.Dataset,
    h3_grid: pt.DataFrame[H3GridWeights],
) -> pt.DataFrame[Nwp]:
    """Vectorized processing of ECMWF dataset to H3 grid."""
    # Convention-sensitive to real ECMWF ENS data: the dim/coord order feeds the ravel + value-join,
    # and the physical units feed Nwp.validate. The offline tests share those assumptions, so after
    # changing this function run the network-gated test manually:
    #     uv run pytest --run-network -m network
    # See
    # <https://openclimatefix.github.io/nged-substation-forecast/architecture/testing/#network-gated-tests>.
    # Precompute latitude and longitude grids
    lat_grid, lon_grid = np.meshgrid(
        ds.latitude.values.astype(np.float32),
        ds.longitude.values.astype(np.float32),
        indexing="ij",
    )
    lat_grid_raveled = lat_grid.ravel()
    lon_grid_raveled = lon_grid.ravel()

    # Iterate over lead_time and ensemble_member to process in chunks.
    dfs: list[pl.DataFrame] = []
    for lead_time in ds.lead_time.values:
        for ensemble_member in ds.ensemble_member.values:
            ds_chunk = ds.sel(lead_time=lead_time, ensemble_member=ensemble_member)
            df = _process_chunk_for_1_lead_time_and_1_ens_member(
                ds_chunk,
                h3_grid,
                lat_grid=lat_grid_raveled,
                lon_grid=lon_grid_raveled,
            ).with_columns(
                ensemble_member=pl.lit(ensemble_member).cast(pl.Int8),
                valid_time=pl.lit(ds_chunk["valid_time"].values).cast(UTC_DATETIME_DTYPE),
            )
            dfs.append(df)

    df = (
        pl.concat(dfs)
        .with_columns(
            nwp_model_id=pl.lit(NwpModelId.ECMWF_ENS_0_25_degree.name).cast(pl.String),
            init_time=pl.lit(ds["init_time"].values).cast(UTC_DATETIME_DTYPE),
            # h3_grid.h3_index (H3GridWeights, joined in above) is UInt64 — that contract is
            # Parquet-backed, so it never had Delta's no-unsigned-integer constraint — but
            # Nwp.h3_index is Int64, so it needs an explicit cast here rather than at
            # H3GridWeights's boundary. The narrowing loses nothing at any resolution: H3 reserves
            # bit 63 of every index as zero, so no H3 value reaches the bit a signed Int64 gives
            # up. See Nwp.h3_index in contracts.weather_schemas for the full argument.
            h3_index=pl.col("h3_index").cast(pl.Int64),
            wind_speed_10m=_calc_wind_speed(height="10m"),
            wind_speed_100m=_calc_wind_speed(height="100m"),
            wind_direction_10m=_calc_wind_direction(height="10m"),
            wind_direction_100m=_calc_wind_direction(height="100m"),
        )
        .drop(cs.matches("^wind_u_.*") | cs.matches("^wind_v_.*"))
    )

    # No sort here: physical row order is delta_store.nwp's job (the single source of truth for
    # on-disk layout), and validation is order-independent.
    return Nwp.validate(df)

dynamical_data.ecmwf_ens.upstream_nulls

Measuring upstream corruption on the raw NWP grid, before the H3 aggregation sees it.

The H3 aggregation renormalises each cell over the grid points that supplied a value, so a corrupt grid point costs only its own share of its cell. That is what makes the stored cells robust, and it is also why counting null cells is a poor proxy for how corrupt the feed was. This module counts the nulls where they arrive.

Classes

UpstreamNullRate dataclass

How much of one ingested NWP run arrived null on the raw grid.

This is the provider channel: the number to quote to Dynamical.org when asking whether their feed is degrading. It counts grid points on the 0.25° lat/lon box we downloaded, before any H3 aggregation.

Read it alongside, never instead of, :class:contracts.weather_schemas.NwpQualityReport, which counts null H3 cells and answers the different question of how much the model lost. The two are not comparable as rates: different units over different populations.

See https://openclimatefix.github.io/nged-substation-forecast/architecture/ecmwf-ens-known-issues/.

Source code in packages/dynamical_data/src/dynamical_data/ecmwf_ens/upstream_nulls.py

@dataclass(frozen=True)
class UpstreamNullRate:
    """How much of one ingested NWP run arrived null on the **raw grid**.

    This is the provider channel: the number to quote to Dynamical.org when asking whether their
    feed is degrading. It counts grid points on the 0.25° lat/lon box we downloaded, before any H3
    aggregation.

    Read it alongside, never instead of, :class:`contracts.weather_schemas.NwpQualityReport`, which
    counts null H3 *cells* and answers the different question of how much the model lost. The two
    are not comparable as rates: different units over different populations.

    See
    <https://openclimatefix.github.io/nged-substation-forecast/architecture/ecmwf-ens-known-issues/>.
    """

    per_variable: pl.DataFrame
    """One row per counted variable, with its ``n_null``, ``n_affected_slices`` and ``n_total``
    grid-point counts, sorted by variable name. Every scalar below is derived from it, so a
    breakdown and a total cannot disagree."""

    @property
    def n_null_nwp_grid_points(self) -> int:
        """Null grid points in the counted variables and steps."""
        return int(self.per_variable["n_null"].sum())

    @property
    def n_total_nwp_grid_points(self) -> int:
        """The denominator: counted variables × ensemble members × counted steps × grid points."""
        return int(self.per_variable["n_total"].sum())

    @property
    def n_affected_nwp_slices(self) -> int:
        """``(variable, ensemble_member, lead_time)`` slices carrying at least one null grid point.

        Separates "one bad slice" from "a hundred" at the same overall rate, which the fraction
        alone cannot.
        """
        return int(self.per_variable["n_affected_slices"].sum())

    @property
    def affected_nwp_variables(self) -> tuple[str, ...]:
        """The counted variables carrying at least one null grid point.

        In ``per_variable``'s row order, which is sorted by variable name because
        :func:`assess_upstream_grid_point_nulls` builds it that way — ``filter`` preserves row
        order rather than imposing one.
        """
        return tuple(self.per_variable.filter(pl.col("n_null") > 0)["variable"])

    @property
    def null_nwp_grid_point_fraction(self) -> float:
        """Null grid points as a fraction of those counted; ``0.0`` when none were counted.

        A run with no step left to count has nothing to measure, and a warning path must not
        raise
        ([rule 7](https://openclimatefix.github.io/nged-substation-forecast/design-philosophy/inherent-stability/#the-rules)).
        """
        if self.n_total_nwp_grid_points == 0:
            return 0.0
        return self.n_null_nwp_grid_points / self.n_total_nwp_grid_points

    @property
    def is_healthy(self) -> bool:
        """True when no counted grid point arrived null."""
        return self.n_null_nwp_grid_points == 0

Attributes

per_variable instance-attribute

One row per counted variable, with its n_null, n_affected_slices and n_total grid-point counts, sorted by variable name. Every scalar below is derived from it, so a breakdown and a total cannot disagree.

n_null_nwp_grid_points property

Null grid points in the counted variables and steps.

n_total_nwp_grid_points property

The denominator: counted variables × ensemble members × counted steps × grid points.

n_affected_nwp_slices property

(variable, ensemble_member, lead_time) slices carrying at least one null grid point.

Separates "one bad slice" from "a hundred" at the same overall rate, which the fraction alone cannot.

affected_nwp_variables property

The counted variables carrying at least one null grid point.

In per_variable's row order, which is sorted by variable name because :func:assess_upstream_grid_point_nulls builds it that way — filter preserves row order rather than imposing one.

null_nwp_grid_point_fraction property

Null grid points as a fraction of those counted; 0.0 when none were counted.

A run with no step left to count has nothing to measure, and a warning path must not raise (rule 7).

is_healthy property

True when no counted grid point arrived null.

Methods:

__init__(per_variable)

Functions:

assess_upstream_grid_point_nulls(ds, variables, exclude_lead_0)

Count nulls on the raw NWP grid of one downloaded run.

Pure and Dagster-free (unit-testable in isolation); the ecmwf_ens asset calls it twice, once per null population, and publishes each result on its own WARN check.

Parameters:

Name	Type	Description	Default
ds	Dataset	One downloaded ECMWF ENS run, as returned by :func:dynamical_data.ecmwf_ens.download.download_ecmwf_ens_data — dimensions (lead_time, ensemble_member, latitude, longitude), with init_time already reduced to a scalar coordinate.	required
variables	Collection[str]	The variables to count over, named as ds names them rather than as the Nwp contract does — the two differ on wind, so :data:dynamical_data.ecmwf_ens.download.ECMWF_ENS_INSTANTANEOUS_VARS exists to be passed here. Their nulls must share one meaning, because a rate pooled over variables with opposite null semantics measures nothing: the asset passes the de-accumulated variables, whose nulls are known upstream corruption, and the instantaneous ones, whose nulls are anomalous, on separate calls.	required
exclude_lead_0	bool	Skip the lead-0 step. True for the de-accumulated variables, which are null there by design, so counting it would report every healthy run as corrupt. False for the instantaneous ones, where lead-0 is an ordinary step and a null in it means what a null in any other step means.	required

Source code in packages/dynamical_data/src/dynamical_data/ecmwf_ens/upstream_nulls.py

def assess_upstream_grid_point_nulls(
    ds: xr.Dataset, variables: Collection[str], exclude_lead_0: bool
) -> UpstreamNullRate:
    """Count nulls on the raw NWP grid of one downloaded run.

    Pure and Dagster-free (unit-testable in isolation); the ``ecmwf_ens`` asset calls it twice, once
    per null population, and publishes each result on its own WARN check.

    Args:
        ds: One downloaded ECMWF ENS run, as returned by
            :func:`dynamical_data.ecmwf_ens.download.download_ecmwf_ens_data` — dimensions
            ``(lead_time, ensemble_member, latitude, longitude)``, with ``init_time`` already
            reduced to a scalar coordinate.
        variables: The variables to count over, named as ``ds`` names them rather than as the
            ``Nwp`` contract does — the two differ on wind, so
            :data:`dynamical_data.ecmwf_ens.download.ECMWF_ENS_INSTANTANEOUS_VARS` exists to be
            passed here. Their nulls must share one meaning, because a rate pooled over variables
            with opposite null semantics measures nothing: the asset passes the de-accumulated
            variables, whose nulls are known upstream corruption, and the instantaneous ones, whose
            nulls are anomalous, on separate calls.
        exclude_lead_0: Skip the lead-0 step. True for the de-accumulated variables, which are null
            there by design, so counting it would report every healthy run as corrupt. False for
            the instantaneous ones, where lead-0 is an ordinary step and a null in it means what a
            null in any other step means.
    """
    # Selected per variable rather than once on `ds`: this runs inside `ecmwf_ens`, whose ECMWF
    # concurrency pool exists because the download is memory-intensive, and slicing the whole
    # dataset would copy all thirteen downloaded variables to read three.
    beyond_lead_0 = ds.lead_time > _LEAD_0
    rows = []
    for name in sorted(variables):
        values = ds[name].isel(lead_time=beyond_lead_0) if exclude_lead_0 else ds[name]
        nulls_per_slice = values.isnull().sum(dim=["latitude", "longitude"])
        rows.append(
            {
                "variable": name,
                "n_null": int(nulls_per_slice.sum()),
                "n_affected_slices": int((nulls_per_slice > 0).sum()),
                "n_total": values.size,
            }
        )
    return UpstreamNullRate(per_variable=pl.DataFrame(rows, schema=_PER_VARIABLE_SCHEMA))