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

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 docs/architecture/testing.md ("Network-gated tests").
    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(timezone.utc).replace(tzinfo=None))

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

    # Cast to xr.Dataset to satisfy the type checker, as indexing with a list is misidentified as
    # returning a DataArray.
    ds = cast(xr.Dataset, 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")

    # Check for empty coordinates before computing bounds to fail gracefully.
    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.keys()
        ]
        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

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 docs/architecture/testing.md ("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.UInt8),
                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(NWP_MODEL_ID_DTYPE),
            init_time=pl.lit(ds["init_time"].values).cast(UTC_DATETIME_DTYPE),
            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)