Ingest all 11 HOSTRADA variables in a compact int16 archive

- Store values as int16 (scale 0.1), rounded to DWD's published
  resolution; transpose + Blosc/zstd5 codecs (2.0 MB vs 6.0 MB per chunk)
- Circular mean for wind direction; per-variable units, time convention
  and resolution as coordinates (fixes DWD's Pa/hPa label, documents
  hour-ending rsds)
- Stage months on disk and write one year at a time; track completed
  years in the store; refuse to overwrite without --resume
- --seed-from copies variables from an older archive instead of
  re-downloading them
- Set zarr fill_value to match _FillValue so unwritten hours are NaN
- Require zarr>=3 / Python>=3.11

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

.gitignore

@@ -2,6 +2,7 @@
 __pycache__/
 *.pyc
 raw_cache/
+staging/
 *.zarr/
 .DS_Store
 *.egg-info

Makefile

@@ -4,21 +4,23 @@ END_MONTH=2026-08
 help:     ## Show this help.
 	@egrep -h '(\s##\s|^##\s)' $(MAKEFILE_LIST) | egrep -v '^--' | awk 'BEGIN {FS = ":.*?## "}; {printf "\033[32m  %-35s\033[0m %s\n", $$1, $$2}'
 
-transpose: ## Download and transpose
-	uv run python -m hostrada_transpose.build_archive \
+ARCHIVE=./hostrada_coarse.zarr
+# Optional: older archive (same grid/time axis) whose variables are copied
+# instead of downloaded again, e.g. make transpose SEED=./hostrada_coarse_tas_rsds.zarr
+SEED=
+BUILD=uv run python -m hostrada_transpose.build_archive \
 		--start ${START_MONTH} --end ${END_MONTH} \
-		--variables tas,rsds \
-		--raw-cache-dir ./raw_cache --archive-path ./hostrada_coarse.zarr
+		--raw-cache-dir ./raw_cache --staging-dir ./staging --archive-path ${ARCHIVE} \
+		$(if ${SEED},--seed-from ${SEED})
+
+transpose: ## Download and transpose all variables
+	${BUILD}
 
 # Example with screen:
 # screen -S "Hostrada transpose" -L -Logfile "hostrada_transpose.log" make transpose
 
 resume: ## Resume already started transpose
-	uv run python -m hostrada_transpose.build_archive \
-		--start ${START_MONTH} --end ${END_MONTH} \
-		--variables tas,rsds \
-		--raw-cache-dir ./raw_cache --archive-path ./hostrada_coarse.zarr \
-		--resume
+	${BUILD} --resume
 
 # TODO: Show available data
 # TODO: Show data for given location / date range

README.md

@@ -29,9 +29,23 @@ with no network access at all.
 ## What it does
 
 1. **Downloads** each monthly national grid via `hostrada4py.hostrada.ensure_month_file` (full-grid mode, not point mode -- aggregation needs every cell).
-2. **Aggregates** it to a coarser grid (`hostrada_transpose.aggregate`): a plain block mean, `AGGREGATION_FACTOR` cells per edge (default 10, i.e. ~10 km from HOSTRADA's verified 1 km x 1 km grid).
-3. **Writes** the aggregated month into its region of a pre-allocated Zarr store (`hostrada_transpose.archive`): a single array shaped `(time, Y, X, variable)`, chunked `(whole time range, 2 coarse cells, 2 coarse cells, all variables)`. Variables are stacked into one array (not one array per variable) so that a point query for every variable at one location is a single chunk read, matching how SimStadt will actually consume this. The time axis is one chunk spanning the entire configured range by default (`TIME_CHUNK_HOURS = None` in `config.py`) -- a point query reads its whole multi-decade series in one chunk read. Measured the write-amplification this causes during incremental ingestion (see `config.py`'s comment): real, but adds only ~15-20 minutes across a full 30-year run, negligible next to the ~20-30s/month network download that otherwise dominates.
-4. **Discards** the raw national-grid download once aggregated (`--keep-raw` to keep it).
+2. **Aggregates** it to a coarser grid (`hostrada_transpose.aggregate`): a plain block mean, `AGGREGATION_FACTOR` cells per edge (default 10, i.e. ~10 km from HOSTRADA's verified 1 km x 1 km grid). Wind direction is averaged as unit vectors (an arithmetic mean of 350° and 10° would be 180°).
+3. **Rounds** the result to the precision DWD publishes (0.1 or 1 unit, see below) -- the archive never claims more detail than the source -- and **stages** it as a small `.npy` file in `--staging-dir`, then discards the raw download (`--keep-raw` to keep it).
+4. **Writes** one calendar year, all variables at once, into a pre-allocated Zarr v3 store (`hostrada_transpose.archive`), then records the year in the store's `completed_years` attribute. The store holds a single array shaped `(time, Y, X, variable)`, chunked `(whole time range, 2 coarse cells, 2 coarse cells, all variables)`: a point query for every variable at one location is a single chunk read of the whole multi-decade series. Writing into such a chunk re-encodes it, which is why writes are batched per year (~32 rewrites of the store, ~5 min each, instead of one per variable and month).
+
+### Storage format
+
+| Aspect | Choice | Why |
+|---|---|---|
+| Values | `int16`, CF `scale_factor=0.1`, `_FillValue=-32768` | DWD's own files are integers with 0.1 or 1 resolution; every range fits (largest: psl, ~10,600 steps). xarray decodes to float with NaN outside Germany. |
+| Layout on disk | `transpose` codec, so each cell's time series is contiguous | Consecutive hours compress far better than neighbouring cells/variables |
+| Compression | Blosc, zstd level 5, byte shuffle | Core Zarr v3 codecs only (readable without numcodecs extras) |
+| Metadata | `var_units`, `var_long_name`, `var_cell_methods`, `var_time_convention`, `var_resolution` coordinates along `variable` | DWD's units/time semantics survive aggregation |
+
+Measured on a real 2x2-cell chunk (tas+rsds, 1995-2026): **2.0 MB** vs. 6.0 MB
+for the former float32/zstd-0 encoding. Blosc level 9 (1.9 MB, 13x slower
+writes) and a numcodecs `Delta` filter (1.9 MB, non-core codec) were not worth
+it.
 
 Run it with:
 
@@ -39,32 +53,56 @@ Run it with:
 uv venv && uv pip install -e .
 python -m hostrada_transpose.build_archive \
     --start 1995-01 --end 2026-08 \
-    --variables tas,rsds,hurs,sfcWind,clt,uhi \
-    --raw-cache-dir ./raw_cache --archive-path ./hostrada_coarse.zarr
+    --raw-cache-dir ./raw_cache --staging-dir ./staging --archive-path ./hostrada_coarse.zarr
+```
+
+or `make transpose`. All variables are ingested by default (`--variables` to
+restrict). An existing archive is never overwritten: `--resume` (`make resume`)
+continues it, skipping completed years and already-staged months.
+`--seed-from OLD.zarr` copies variables an older archive already has (same
+grid and time axis) instead of downloading them again, e.g. to upgrade the
+former tas/rsds-only archive:
+
+```bash
+mv hostrada_coarse.zarr hostrada_coarse_tas_rsds.zarr
+make transpose SEED=./hostrada_coarse_tas_rsds.zarr
 ```
 
 Reading one point's full time series back out looks like:
 
 ```python
+import pyproj
 import xarray as xr
+
 ds = xr.open_zarr("hostrada_coarse.zarr")
-point = ds["weather"].sel(Y=..., X=..., method="nearest")  # all variables, full time range, one chunk read
-tas = point.sel(variable="tas")
+x, y = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3034", always_xy=True).transform(9.1737, 48.7802)
+df = ds["weather"].sel(X=x, Y=y, method="nearest").to_pandas()  # all variables, full time range, one chunk read
+ds["var_units"].to_pandas()                                     # units per variable
 ```
 
-This is a long-running (hours), high-bandwidth (see estimate below), one-time
-batch job. Run it unattended on a machine with a stable connection, not from
-an interactive SimStadt workflow step.
+This is a long-running (about a day), high-bandwidth (see estimate below),
+one-time batch job. Run it unattended on a machine with a stable connection,
+not from an interactive SimStadt workflow step.
 
 ## Variables
 
-| Variable | Why |
-|---|---|
-| `tas` | Ambient temperature -- confirmed required by SimStadt's `GreenWaterProcessor`/`EvapoTranspirationCalculator` |
-| `rsds` | GHI -- same, plus the base input for the derived DHI estimate |
-| `hurs` | Relative humidity -- same |
-| `sfcWind` | Wind speed -- same |
-| `clt`, `uhi` | Optional: only needed to enable `HostradaDiffuse.estimate(apply_weather_correction=True)`'s refinement of the DHI estimate. Drop them from `config.ALL_VARIABLES` if that refinement isn't worth the extra ingestion volume. |
+All 11 HOSTRADA variables. Resolution = precision of DWD's published values,
+to which the 10 km means are rounded. All values are UTC. All are
+instantaneous at the timestamp, except `rsds`.
+
+| Variable | Units | Resolution | Note |
+|---|---|---|---|
+| `tas` | °C | 0.1 | Ambient temperature, required by SimStadt's GreenWater |
+| `rsds` | W/m² | 0.1 | GHI. **Mean over the hour *ending* at the timestamp.** DWD's NetCDF `time_bnds` claim the hour starting at it; a clear-sky fit over 1995-2026 confirms hour-ending, as in DWD's PDF description. Use `t - 30 min` for solar geometry. |
+| `hurs` | % | 0.1 | Relative humidity, required by GreenWater |
+| `sfcWind` | m/s | 0.1 | Wind speed at 10 m, required by GreenWater |
+| `sfcWind_direction` | ° | 1 | Circular (unit-vector) mean |
+| `clt` | octa | 1 | Cloud cover; used by `HostradaDiffuse(apply_weather_correction=True)` |
+| `uhi` | °C | 0.1 | Urban heat island intensity; same |
+| `tdew` | °C | 0.1 | Dew point (DWD's long_name wrongly says "Daily Mean") |
+| `mixr` | g/kg | 0.1 | Water vapour mixing ratio |
+| `ps` | hPa | 1 | Surface pressure. DWD's files say `Pa`, but the values are hPa. |
+| `psl` | hPa | 1 | Sea-level pressure, same unit issue |
 
 **Precipitation is deliberately out of scope here.** It's required by
 GreenWater but absent from both hostrada4py providers: not in DWD
@@ -125,36 +163,39 @@ Confirmed with station 00183 (Arkona, a Baltic-coast lighthouse).
 `compare_station_to_model` now raises instead of silently returning garbage
 for such stations -- keep that guard if this logic is ever copied elsewhere.
 
-## Data volume (30 years, verified against real cached files)
+## Data volume
 
 HOSTRADA's grid is a verified uniform 1000 m x 1000 m EPSG:3034 raster,
-938 x 720 cells. Per-variable annual download size, measured from a real
-2025 cache:
-
-| Variable | GiB/year |
-|---|---|
-| `tas` | 1.84 |
-| `rsds` | 1.81 |
-| `hurs` | 2.61 |
-| `sfcWind` | 2.18 |
-| `clt` | 0.54 |
-| `uhi` | not sampled -- assumed similar to the others (~1.8) |
-
-**~10.8 GiB/year for the 6 configured variables -> ~324 GB one-time
-ingestion for 1995-2024.** This is the unavoidable cost of reading the
-source files once (chunking forces reading ~100% of the bytes regardless of
-how little of the data is ultimately kept); it is paid once, not per city,
-and not repeated on every run.
-
-Resulting archive after 10x spatial aggregation (93 x 72 coarse cells):
-uncompressed float32 estimate is on the order of several GiB for all 6
-variables over 30 years; real post-compression size is unmeasured --
-run a 1-2 year ingestion first and check before committing to the full
-30-year pull.
+938 x 720 cells. Download volume for 1995-01..2026-08 (380 monthly files per
+variable), summed from DWD's directory listings on 2026-10-07:
+
+| Variable | GiB | Variable | GiB |
+|---|---|---|---|
+| `tas` | 57.8 | `clt` | 17.3 |
+| `rsds` | 59.0 | `uhi` | 20.4 |
+| `hurs` | 83.2 | `tdew` | 55.1 |
+| `sfcWind` | 69.2 | `mixr` | 36.0 |
+| `sfcWind_direction` | 60.6 | `ps` | 62.5 |
+| | | `psl` | 13.9 |
+
+**535 GiB in total, 418 GiB when tas/rsds are seeded from the older archive.**
+This is the unavoidable cost of reading the source files once (chunking
+forces reading ~100% of the bytes regardless of how little of the data is
+ultimately kept); it is paid once, not per city. Measured throughput on a
+test month: ~25 s per file including aggregation.
+
+Resulting archive (93 x 72 coarse cells, ~46% of them outside Germany):
+estimated ~10 GB for all 11 variables, i.e. **~10 MB per location** for the
+full 1995-2026 series over HTTP (one chunk, 2x2 cells). This extrapolates
+from tas/rsds chunks (~1 MB per variable); ps/psl/clt should compress better,
+wind direction worse. Measure on the finished archive.
+
+Disk and memory during the build: one year in memory (~1.3 GB int16), staged
+months for one year (~1.2 GB), at most one raw month file (~160 MB).
 
 ## Open questions / not yet done
 
 - Real, multi-station, multi-year DHI bias/RMSE characterization (see above).
 - Precipitation: add `total_precipitation` to `hostrada4py`'s `CERRA_VARIABLES`, then decide how/whether to fold a 5.5 km CERRA series into this archive's output alongside the 10 km HOSTRADA-derived one.
-- `build_archive.py`'s `--resume` flag is a cheap single-cell NaN check, not real ingestion bookkeeping -- fine for restarting after a crash, not a substitute for a proper manifest if this is run incrementally over time.
-- Final archive compression ratio is unmeasured.
+- The time axis is fixed when the archive is created. Appending newer months requires a rebuild (seeding everything from the old archive makes that download-free except for the new months -- not implemented yet).
+- `sfcWind_direction`: whether DWD uses 0° for calm (which would bias the circular mean) is unverified.

hostrada_transpose/aggregate.py

 """Spatial block-aggregation of one HOSTRADA monthly grid file."""
 from __future__ import annotations
 
+import numpy as np
 import xarray as xr
 from hostrada4py.hostrada import find_variable
 
-from .config import AGGREGATION_FACTOR
+from .config import AGGREGATION_FACTOR, STORAGE_FILL, STORAGE_SCALE, VARIABLES, VariableSpec
 
 
 def aggregate_month(
@@ -22,10 +23,21 @@ def aggregate_month(
     expected along the coastline/border, where roughly half the national
     bounding box has no land data at all (verified: ~50% NaN at 1 km,
     ~46% NaN after 10x aggregation).
+
+    Angles (``aggregation="circular_mean"``) are averaged as unit vectors and
+    returned in [0, 360).
     """
     name = find_variable(var, ds)
     da = ds[name]
-    coarse = da.coarsen(X=factor, Y=factor, boundary="trim").mean()
+
+    def block_mean(a: xr.DataArray) -> xr.DataArray:
+        return a.coarsen(X=factor, Y=factor, boundary="trim").mean()
+
+    if VARIABLES[var].aggregation == "circular_mean":
+        rad = np.deg2rad(da)
+        coarse = np.rad2deg(np.arctan2(block_mean(np.sin(rad)), block_mean(np.cos(rad)))) % 360
+    else:
+        coarse = block_mean(da)
     coarse.name = var
     coarse.attrs = dict(da.attrs)
     # HOSTRADA attaches 2-D lon/lat auxiliary coords (dims Y, X). coarsen()
@@ -37,3 +49,13 @@ def aggregate_month(
     if extra_coords:
         coarse = coarse.drop_vars(extra_coords)
     return coarse
+
+
+def to_storage(values: np.ndarray, spec: VariableSpec) -> np.ndarray:
+    """Round to the source resolution and encode as int16 steps of STORAGE_SCALE."""
+    values = np.asarray(values, dtype="float64")
+    rounded = np.round(values / spec.resolution) * spec.resolution
+    if spec.aggregation == "circular_mean":
+        rounded = rounded % 360  # 359.6 rounds to 360, i.e. 0
+    steps = np.round(rounded / STORAGE_SCALE)
+    return np.where(np.isnan(steps), STORAGE_FILL, steps).astype("int16")

hostrada_transpose/archive.py

-"""Create and incrementally fill the time-major coarse-grid Zarr archive.
+"""Create and fill the time-major coarse-grid Zarr archive.
 
-The archive is created once as an empty, NaN-filled template whose full
-shape and chunk layout are known up front (we know the final time range and
-variable list before ingestion starts). Each monthly ingestion step then
-writes only its own region -- one variable's slice of one time range --
-without ever reading or rewriting unrelated data. This is what makes
-building a 30-year archive tractable without holding 30 years of data in
-memory at once.
+The archive is created once as an empty template whose full shape and chunk
+layout are known up front (we know the final time range and variable list
+before ingestion starts). Each ingestion step then writes one calendar year
+for all variables, without ever reading or rewriting unrelated data. This is
+what makes building a 30-year archive tractable without holding 30 years of
+data in memory at once.
 
 Variables are stacked into a single array along a ``variable`` dimension
 (rather than one array per variable) so that a point query for "every
 variable at this location" is one chunk read, matching how SimStadt will
 actually consume this archive (temperature+GHI+humidity+wind together, per
-point). Measured impact of this plus a time axis chunked as one single
-multi-decade block: see config.py's comment on TIME_CHUNK_HOURS.
+point).
+
+Storage: int16 steps of 0.1 (CF ``scale_factor``, decoded by xarray), a
+``transpose`` codec so that each cell's time series is contiguous on disk,
+and Blosc/zstd with byte shuffle. All three are core Zarr v3 codecs.
+Measured on a real 2x2-cell chunk with tas+rsds over 1995-2026: 2.0 MB vs.
+6.0 MB for float32/zstd0, writing in 0.05 s. Blosc clevel 9 (1.9 MB, 13x
+slower writes) and numcodecs Delta (1.9 MB, not a core codec) weren't worth
+it.
+
+Per-variable metadata (units, long_name, cell_methods, time convention,
+source resolution) is stored as coordinates along ``variable``, e.g.
+``ds.var_units.sel(variable="rsds")``.
 """
 from __future__ import annotations
 
+import warnings
 from pathlib import Path
 from typing import Mapping, Sequence
 
@@ -24,10 +35,21 @@ import dask.array as da
 import numpy as np
 import pandas as pd
 import xarray as xr
-
-from .config import SPATIAL_CHUNK_CELLS, TIME_CHUNK_HOURS
+import zarr
+from zarr.codecs import BloscCodec, TransposeCodec
+
+from .config import (
+    AGGREGATION_FACTOR,
+    SPATIAL_CHUNK_CELLS,
+    STORAGE_DTYPE,
+    STORAGE_FILL,
+    STORAGE_SCALE,
+    TIME_CHUNK_HOURS,
+    VARIABLES,
+)
 
 ARRAY_NAME = "weather"
+COMPLETED_YEARS_ATTR = "completed_years"
 
 
 def build_time_index(start: str, end: str) -> pd.DatetimeIndex:
@@ -46,8 +68,8 @@ def create_archive_template(
 
     ``compute=False`` means the dask-backed NaN placeholder array is never
     actually materialised -- Zarr only records the store's shape, dtype and
-    chunking. Real values are filled in later, one (variable, month) slice
-    at a time, by :func:`write_month`.
+    chunking. Real values are filled in later, one year at a time, by
+    :func:`write_block`.
     """
     n_time = len(time_index)
     y_vals = coarse_coords["Y"]
@@ -55,55 +77,68 @@ def create_archive_template(
     ny, nx = len(y_vals), len(x_vals)
     n_var = len(variables)
     time_chunk = n_time if TIME_CHUNK_HOURS is None else min(TIME_CHUNK_HOURS, n_time)
-    y_chunk = min(SPATIAL_CHUNK_CELLS, ny)
-    x_chunk = min(SPATIAL_CHUNK_CELLS, nx)
-
-    placeholder = da.full(
-        (n_time, ny, nx, n_var),
-        np.nan,
-        dtype="float32",
-        chunks=(time_chunk, y_chunk, x_chunk, n_var),
-    )
+    chunks = (time_chunk, min(SPATIAL_CHUNK_CELLS, ny), min(SPATIAL_CHUNK_CELLS, nx), n_var)
+
+    placeholder = da.full((n_time, ny, nx, n_var), np.nan, dtype="float32", chunks=chunks)
+    specs = [VARIABLES[v] for v in variables]
     template = xr.Dataset(
         {ARRAY_NAME: (("time", "Y", "X", "variable"), placeholder)},
-        coords={"time": time_index, "Y": y_vals, "X": x_vals, "variable": list(variables)},
+        coords={
+            "time": time_index,
+            "Y": ("Y", y_vals, {"units": "m", "standard_name": "projection_y_coordinate"}),
+            "X": ("X", x_vals, {"units": "m", "standard_name": "projection_x_coordinate"}),
+            "variable": list(variables),
+            "var_units": ("variable", [s.units for s in specs]),
+            "var_long_name": ("variable", [s.long_name for s in specs]),
+            "var_cell_methods": ("variable", [s.cell_methods for s in specs]),
+            "var_time_convention": ("variable", [s.time_convention for s in specs]),
+            "var_resolution": ("variable", [s.resolution for s in specs]),
+        },
+        attrs={
+            "title": "HOSTRADA v1.0, aggregated to ~10 km, time-major",
+            "source": "DWD HOSTRADA - High-resolution grids of hourly variables for Germany, v1.0",
+            "crs": "EPSG:3034",
+            "aggregation": f"{AGGREGATION_FACTOR}x{AGGREGATION_FACTOR} block mean of the 1 km grid "
+            "(circular mean for wind direction), rounded to the source resolution",
+            "time_zone": "UTC",
+            COMPLETED_YEARS_ATTR: [],
+        },
     )
-    template.to_zarr(path, mode="w", compute=False)
+    encoding = {
+        ARRAY_NAME: {
+            "dtype": STORAGE_DTYPE,
+            "scale_factor": STORAGE_SCALE,
+            "_FillValue": STORAGE_FILL,
+            # Zarr's own fill value, used for never-written elements. Must
+            # match _FillValue, or unwritten hours decode as 0.0, not NaN.
+            "fill_value": STORAGE_FILL,
+            "chunks": chunks,
+            "filters": [TransposeCodec(order=(1, 2, 3, 0))],
+            "compressors": [BloscCodec(cname="zstd", clevel=5, shuffle="shuffle")],
+        }
+    }
+    template.to_zarr(path, mode="w", compute=False, encoding=encoding, zarr_format=3)
+
+
+def write_block(path: Path, time_slice: slice, block: np.ndarray) -> None:
+    """Write already-encoded int16 data for all variables over ``time_slice``.
+
+    Bypasses xarray on purpose: ``block`` is already in storage units (see
+    :func:`hostrada_transpose.aggregate.to_storage`), which avoids a float
+    copy of a year of data.
+    """
+    zarr.open_array(Path(path) / ARRAY_NAME, mode="r+")[time_slice] = block
 
 
-def write_month(
-    path: Path,
-    var: str,
-    variables: Sequence[str],
-    month_data: xr.DataArray,
-    time_index: pd.DatetimeIndex,
-) -> None:
-    """Write one (variable, month) slice into its region of the archive.
+def completed_years(path: Path) -> set[int]:
+    return set(zarr.open_group(path, mode="r").attrs.get(COMPLETED_YEARS_ATTR, []))
 
-    Only ``var``'s data is written, as a width-1 slice along the ``variable``
-    axis -- the other variables already/later occupying that axis in the
-    same chunk are left untouched by this call (Zarr handles the
-    read-modify-write of the rest of the chunk internally).
-    """
-    first_time = month_data["time"].values[0]
-    start_idx = time_index.get_loc(first_time)
-    if not isinstance(start_idx, (int, np.integer)):
-        raise ValueError(
-            f"Month start {first_time} is not aligned to the archive's hourly "
-            "time index -- check for gaps or a mismatched time zone."
-        )
-    end_idx = start_idx + month_data.sizes["time"]
-    var_idx = list(variables).index(var)
-
-    expanded = month_data.astype("float32").expand_dims(variable=[var])
-    ds = expanded.to_dataset(name=ARRAY_NAME)
-    ds = ds.drop_vars(["time", "Y", "X", "variable"], errors="ignore")
-    ds.to_zarr(
-        path,
-        region={
-            "time": slice(start_idx, end_idx),
-            "Y": slice(None),
-            "X": slice(None),
-            "variable": slice(var_idx, var_idx + 1),
-        },
-    )
+
+def mark_year_completed(path: Path, year: int) -> None:
+    group = zarr.open_group(path, mode="r+")
+    done = sorted(set(group.attrs.get(COMPLETED_YEARS_ATTR, [])) | {year})
+    group.attrs[COMPLETED_YEARS_ATTR] = done
+    with warnings.catch_warnings():
+        # Consolidated metadata isn't in the v3 spec yet, but xarray uses it.
+        warnings.filterwarnings("ignore", message="Consolidated metadata")
+        zarr.consolidate_metadata(path)

hostrada_transpose/build_archive.py

@@ -2,48 +2,79 @@
 
 Usage:
     python -m hostrada_transpose.build_archive \\
-        --start 1995-01 --end 2024-12 \\
-        --variables tas,rsds,hurs,sfcWind,clt,uhi \\
-        --raw-cache-dir ./raw_cache --archive-path ./hostrada_coarse.zarr
+        --start 1995-01 --end 2026-08 \\
+        --raw-cache-dir ./raw_cache --archive-path ./hostrada_coarse.zarr \\
+        [--seed-from ./old_archive.zarr] [--resume]
 
 This is a long-running, one-time (or rarely-repeated) batch job -- expect it
-to run for hours and to transfer on the order of 150-250 GB from DWD for the
-default 30-year, 6-variable configuration (see README.md for the estimate
-and its assumptions). It is meant to be run unattended, not from an
-interactive SimStadt workflow step.
+to run for a day or more and to transfer ~535 GiB from DWD for all 11
+variables over 1995-01..2026-08 (sum of the published file sizes, checked
+2026-10). It is meant to be run unattended, not from an interactive
+SimStadt workflow step.
+
+Work is done one calendar year at a time. Each (variable, month) is
+downloaded, aggregated, encoded and staged as a small .npy file in
+``--staging-dir`` (the raw download is then deleted). Once all months of a
+year are staged, the year is written to the archive for all variables at
+once, recorded in the archive's ``completed_years`` attribute, and its
+staged files are removed. An interrupted run restarted with ``--resume``
+therefore skips completed years and already-staged months.
 """
 from __future__ import annotations
 
 import argparse
 import logging
+import os
 from pathlib import Path
 
+import numpy as np
+import pandas as pd
+import xarray as xr
+
 import hostrada4py.hostrada as hs
 
-from .aggregate import aggregate_month
-from .archive import build_time_index, create_archive_template, write_month
+from .aggregate import aggregate_month, to_storage
+from .archive import (
+    ARRAY_NAME,
+    build_time_index,
+    completed_years,
+    create_archive_template,
+    mark_year_completed,
+    write_block,
+)
 from .config import (
     ALL_VARIABLES,
     DEFAULT_ARCHIVE_PATH,
     DEFAULT_END,
     DEFAULT_RAW_CACHE_DIR,
+    DEFAULT_STAGING_DIR,
     DEFAULT_START,
+    STORAGE_FILL,
+    VARIABLES,
 )
 
 logger = logging.getLogger("hostrada_transpose")
 
 
 def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
-    parser = argparse.ArgumentParser(description=__doc__)
+    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
     parser.add_argument("--start", default=DEFAULT_START, help="First month, e.g. 1995-01")
-    parser.add_argument("--end", default=DEFAULT_END, help="Last month (inclusive), e.g. 2024-12")
+    parser.add_argument("--end", default=DEFAULT_END, help="Last month (inclusive), e.g. 2026-08")
     parser.add_argument(
         "--variables",
         default=",".join(ALL_VARIABLES),
-        help="Comma-separated canonical HOSTRADA variable names",
+        help="Comma-separated canonical HOSTRADA variable names (default: all)",
     )
     parser.add_argument("--raw-cache-dir", type=Path, default=DEFAULT_RAW_CACHE_DIR)
+    parser.add_argument("--staging-dir", type=Path, default=DEFAULT_STAGING_DIR)
     parser.add_argument("--archive-path", type=Path, default=DEFAULT_ARCHIVE_PATH)
+    parser.add_argument(
+        "--seed-from",
+        type=Path,
+        help="Existing archive (same grid and time axis, any encoding) to copy "
+        "variables from instead of downloading them again, e.g. an older "
+        "tas/rsds-only archive.",
+    )
     parser.add_argument(
         "--keep-raw",
         action="store_true",
@@ -53,70 +84,127 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
     parser.add_argument(
         "--resume",
         action="store_true",
-        help="Skip (var, year, month) combinations whose region appears already "
-        "non-NaN in the archive. Cheap safety net for restarting an interrupted run "
-        "-- NOT a substitute for keeping real ingestion-progress bookkeeping if this "
-        "job is run repeatedly.",
+        help="Continue an existing archive: skip years recorded as completed and "
+        "months already staged. Without it, an existing archive is never touched.",
     )
     return parser.parse_args(argv)
 
 
-def _first_coarse_cell_is_written(archive_path: Path, var: str, variables: list[str], start_idx: int) -> bool:
-    import numpy as np
-    import xarray as xr
+def _staged_month(var: str, year: int, month: int, expected_time: pd.DatetimeIndex, args) -> np.ndarray:
+    """Encoded coarse grid for one (variable, month), downloading it if not staged yet."""
+    staged = args.staging_dir / f"{var}_{year:04d}{month:02d}.npy"
+    if staged.exists():
+        return np.load(staged)
+
+    logger.info("Fetch %s %04d-%02d", var, year, month)
+    raw_path = hs.ensure_month_file(var, year, month, args.raw_cache_dir)
+    with hs.read_month_file(raw_path) as ds:
+        coarse = aggregate_month(ds, var)
+        month_time = pd.DatetimeIndex(coarse["time"].values)
+        if not month_time.equals(expected_time):
+            raise ValueError(
+                f"{var} {year:04d}-{month:02d}: time axis does not match the archive's hourly "
+                f"index ({len(month_time)} vs {len(expected_time)} steps, first "
+                f"{month_time[0] if len(month_time) else None} vs {expected_time[0]})."
+            )
+        encoded = to_storage(coarse.transpose("time", "Y", "X").values, VARIABLES[var])
+
+    args.staging_dir.mkdir(parents=True, exist_ok=True)
+    tmp = staged.with_suffix(".tmp.npy")
+    np.save(tmp, encoded)
+    os.replace(tmp, staged)
+    if not args.keep_raw:
+        raw_path.unlink(missing_ok=True)
+    return encoded
+
+
+def _probe_coarse_coords(var: str, time_index: pd.DatetimeIndex, args) -> dict[str, np.ndarray]:
+    """Coarse Y/X values from one real month (all variables share one grid)."""
+    year, month = time_index[0].year, time_index[0].month
+    raw_path = hs.ensure_month_file(var, year, month, args.raw_cache_dir)
+    with hs.read_month_file(raw_path) as ds:
+        coarse = aggregate_month(ds, var)
+    # The raw file is kept: the first year's ingestion picks it up again.
+    return {"Y": coarse["Y"].values, "X": coarse["X"].values}
 
-    from .archive import ARRAY_NAME
 
-    var_idx = variables.index(var)
-    with xr.open_zarr(archive_path) as ds:
-        sample = (
-            ds[ARRAY_NAME]
-            .isel(time=start_idx, Y=slice(0, 5), X=slice(0, 5), variable=var_idx)
-            .values
-        )
-    return bool(np.isfinite(sample).any())
+def _check_seed(seed: xr.Dataset, coords: dict[str, np.ndarray], time_index: pd.DatetimeIndex) -> None:
+    for dim in ("Y", "X"):
+        if not np.allclose(seed[dim].values, coords[dim]):
+            raise ValueError(f"--seed-from archive has a different {dim} grid.")
+    missing = time_index.difference(pd.DatetimeIndex(seed["time"].values))
+    if len(missing):
+        raise ValueError(f"--seed-from archive lacks {len(missing)} hours, e.g. {missing[0]}.")
 
 
 def main(argv: list[str] | None = None) -> None:
     logging.basicConfig(level=logging.INFO, format="%(asctime)s %(message)s")
     args = parse_args(argv)
     variables = [v.strip() for v in args.variables.split(",") if v.strip()]
+    unknown = sorted(set(variables) - set(VARIABLES))
+    if unknown:
+        raise SystemExit(f"Unknown variables {unknown}; known: {', '.join(VARIABLES)}")
 
     time_index = build_time_index(args.start, args.end)
     logger.info("Archive time range: %s - %s (%d hours)", time_index[0], time_index[-1], len(time_index))
 
-    if not args.archive_path.exists():
-        # Grab one real month's grid up front purely to read the coarse
-        # Y/X coordinate values the template needs -- the first variable
-        # in the list stands in for all of them, since they share one grid.
-        probe_year, probe_month = time_index[0].year, time_index[0].month
-        probe_path = hs.ensure_month_file(variables[0], probe_year, probe_month, args.raw_cache_dir)
-        with hs.read_month_file(probe_path) as probe_ds:
-            probe_coarse = aggregate_month(probe_ds, variables[0])
-        coarse_coords = {"Y": probe_coarse["Y"].values, "X": probe_coarse["X"].values}
-        create_archive_template(args.archive_path, variables, time_index, coarse_coords)
-        logger.info("Created archive template at %s", args.archive_path)
-        if not args.keep_raw:
-            probe_path.unlink(missing_ok=True)
+    seed = xr.open_zarr(args.seed_from) if args.seed_from else None
+    seeded = [v for v in variables if seed is not None and v in seed["variable"].values]
+    if seed is not None:
+        logger.info("Copying %s from %s instead of downloading", seeded or "nothing", args.seed_from)
 
-    for var in variables:
-        for year, month in hs.month_range(time_index[0], time_index[-1]):
-            start_idx = time_index.get_loc(
-                time_index[(time_index.year == year) & (time_index.month == month)][0]
+    if args.archive_path.exists():
+        if not args.resume:
+            raise SystemExit(
+                f"{args.archive_path} already exists. Pass --resume to continue it, "
+                "or move it away (e.g. to use it with --seed-from)."
             )
-            if args.resume and _first_coarse_cell_is_written(args.archive_path, var, variables, start_idx):
-                logger.info("Skip %s %04d-%02d (already written)", var, year, month)
+        with xr.open_zarr(args.archive_path) as existing:
+            if list(existing["variable"].values) != variables or not pd.DatetimeIndex(
+                existing["time"].values
+            ).equals(time_index):
+                raise SystemExit(f"{args.archive_path} was created with other variables or another time range.")
+            coords = {"Y": existing["Y"].values, "X": existing["X"].values}
+    else:
+        if seed is not None:
+            coords = {"Y": seed["Y"].values, "X": seed["X"].values}
+        else:
+            coords = _probe_coarse_coords(variables[0], time_index, args)
+        create_archive_template(args.archive_path, variables, time_index, coords)
+        logger.info("Created archive template at %s", args.archive_path)
+    if seed is not None:
+        _check_seed(seed, coords, time_index)
+
+    done = completed_years(args.archive_path)
+    for year in sorted(set(time_index.year)):
+        if year in done:
+            logger.info("Skip %d (already completed)", year)
+            continue
+        rows = np.flatnonzero(time_index.year == year)
+        year_slice = slice(int(rows[0]), int(rows[-1]) + 1)
+        year_time = time_index[year_slice]
+        block = np.full((len(year_time), len(coords["Y"]), len(coords["X"]), len(variables)), STORAGE_FILL, dtype="int16")
+
+        if seeded:
+            values = seed[ARRAY_NAME].sel(time=year_time, variable=seeded).transpose("time", "Y", "X", "variable").values
+            for k, var in enumerate(seeded):
+                block[..., variables.index(var)] = to_storage(values[..., k], VARIABLES[var])
+            del values
+
+        months = sorted(set(year_time.month))
+        for i, var in enumerate(variables):
+            if var in seeded:
                 continue
-
-            logger.info("Fetch %s %04d-%02d", var, year, month)
-            raw_path = hs.ensure_month_file(var, year, month, args.raw_cache_dir)
-            with hs.read_month_file(raw_path) as ds:
-                coarse = aggregate_month(ds, var)
-            write_month(args.archive_path, var, variables, coarse, time_index)
-            logger.info("Wrote %s %04d-%02d", var, year, month)
-
-            if not args.keep_raw:
-                raw_path.unlink(missing_ok=True)
+            for month in months:
+                month_rows = np.flatnonzero(year_time.month == month)
+                block[month_rows, :, :, i] = _staged_month(var, year, month, year_time[month_rows], args)
+
+        logger.info("Write %d (%d variables)", year, len(variables))
+        write_block(args.archive_path, year_slice, block)
+        mark_year_completed(args.archive_path, year)
+        for staged in args.staging_dir.glob(f"*_{year:04d}??.npy"):
+            staged.unlink()
+        logger.info("Completed %d", year)
 
 
 if __name__ == "__main__":

hostrada_transpose/config.py

 from __future__ import annotations
 
+from dataclasses import dataclass
 from pathlib import Path
 
 # --- Variables -----------------------------------------------------------
+
+
+@dataclass(frozen=True)
+class VariableSpec:
+    """How one HOSTRADA variable is aggregated, quantised and described.
+
+    ``resolution`` is the precision DWD itself publishes (int32 values, with
+    ``scale_factor=0.1`` for some variables, verified on the June 2025 files).
+    The 10 km block mean is rounded back to it, so the archive never claims
+    more precision than the source.
+    """
+
+    units: str
+    long_name: str
+    resolution: float
+    cell_methods: str = "time: point"
+    time_convention: str = "instantaneous value at the timestamp (UTC)"
+    aggregation: str = "mean"  # "mean" or "circular_mean" (angles in degrees)
+
+
+_HOUR_ENDING = (
+    "mean over the hour ending at the timestamp (UTC). DWD's description says "
+    "'sum over the last hour'; the NetCDF time_bnds claim the hour starting at "
+    "the timestamp, but a clear-sky fit (1995-2026) confirms hour-ending."
+)
+
+VARIABLES: dict[str, VariableSpec] = {
+    "tas": VariableSpec("degC", "Near-Surface Air Temperature", 0.1),
+    "rsds": VariableSpec(
+        "W m-2", "Surface Downwelling Shortwave Radiation (GHI)", 0.1,
+        cell_methods="time: mean", time_convention=_HOUR_ENDING,
+    ),
+    "hurs": VariableSpec("%", "Near-Surface Relative Humidity", 0.1),
+    "sfcWind": VariableSpec("m s-1", "Near-Surface Wind Speed (10 m)", 0.1),
+    # Block-averaging angles arithmetically is wrong (mean of 350 and 10 would
+    # be 180), so directions are averaged as unit vectors.
+    "sfcWind_direction": VariableSpec(
+        "degree", "Near-Surface Wind Direction (10 m)", 1.0, aggregation="circular_mean",
+    ),
+    "clt": VariableSpec("octa", "Total Cloud Fraction", 1.0),
+    "uhi": VariableSpec("degC", "Canopy Urban Heat Island Intensity", 0.1),
+    # DWD's long_name says "Daily Mean Dew-point Temperature", but the data
+    # are hourly like every other HOSTRADA variable.
+    "tdew": VariableSpec("degC", "Near-Surface Dew-point Temperature", 0.1),
+    "mixr": VariableSpec("g kg-1", "Near-Surface Water Vapor Mixing Ratio", 0.1),
+    # DWD's files say units="Pa", but the stored values are 750-1050: hPa.
+    "ps": VariableSpec("hPa", "Surface Air Pressure", 1.0),
+    "psl": VariableSpec("hPa", "Sea Level Pressure", 1.0),
+}
+
 # tas, rsds, hurs, sfcWind: confirmed against SimStadt's GreenWater module
 # (de.hft.stuttgart.simstadt2.greenwater) -- ambient temperature, GHI,
 # relative humidity, wind speed. DHI is NOT a separate HOSTRADA variable; it
 # is derived from rsds via pvlib's Erbs/Erbs-Driesse decomposition
-# (hostrada4py.hostradaDiffuse.HostradaDiffuse).
-CORE_VARIABLES: tuple[str, ...] = ("tas", "rsds", "hurs", "sfcWind")
-
-# clt, uhi: not required by SimStadt directly, only pulled in to enable
-# HostradaDiffuse's apply_weather_correction=True path, which refines the
-# Erbs-Driesse DHI estimate using cloud cover and urban-heat-island signals.
-# Drop these from ALL_VARIABLES if that refinement isn't worth the extra
-# ~2 GiB/year of ingestion.
-DHI_CORRECTION_VARIABLES: tuple[str, ...] = ("clt", "uhi")
-
-ALL_VARIABLES: tuple[str, ...] = CORE_VARIABLES + DHI_CORRECTION_VARIABLES
+# (hostrada4py.hostradaDiffuse.HostradaDiffuse). clt, uhi enable its
+# apply_weather_correction=True path. The rest completes the HOSTRADA set
+# (e.g. for TMY3/EPW exports).
+ALL_VARIABLES: tuple[str, ...] = tuple(VARIABLES)
 
 # Precipitation is intentionally NOT handled here. Neither the DWD HOSTRADA
 # provider nor hostrada4py's CERRA provider currently expose it (confirmed by
@@ -26,6 +70,15 @@ ALL_VARIABLES: tuple[str, ...] = CORE_VARIABLES + DHI_CORRECTION_VARIABLES
 # CERRA_VARIABLES). Adding it belongs in hostrada4py itself, not here, since
 # it's provider-level plumbing other hostrada4py consumers would want too.
 
+# --- Storage encoding ---------------------------------------------------
+# Every variable is stored as int16 in steps of 0.1 (CF scale_factor, decoded
+# transparently by xarray). Variables with a coarser source resolution (1 hPa,
+# 1 octa, 1 degree) are simply multiples of 10. All ranges fit: the largest is
+# psl, ~10,600 steps. -32768 marks missing data (outside Germany).
+STORAGE_DTYPE = "int16"
+STORAGE_SCALE = 0.1
+STORAGE_FILL = -32768
+
 # --- Spatial aggregation ---------------------------------------------------
 # HOSTRADA's native grid is a verified uniform 1000 m x 1000 m EPSG:3034
 # raster (938 x 720 cells). AGGREGATION_FACTOR=10 -> ~10 km cells, a plain
@@ -46,15 +99,10 @@ DEFAULT_END = "2026-08"
 # --- Archive chunking -------------------------------------------------
 # TIME_CHUNK_HOURS = None means "the whole configured time range in one
 # chunk" -- a single point query then reads its full multi-decade series in
-# one chunk read. Measured the write-amplification this causes during
-# incremental monthly ingestion (read-modify-write of an increasingly-full
-# chunk): ~0.9s/write when a 10-year chunk is nearly empty, rising to
-# ~1.9s/write once it's full (120-step benchmark). That adds on the order of
-# 15-20 minutes across a full 30-year ingestion -- negligible next to the
-# ~20-30s/month network download that otherwise dominates ingestion time
-# at this archive's few-GB total size. Set this to an hour count (e.g. one
-# decade) instead if this is ever re-run against a much larger archive
-# where that overhead would no longer be negligible.
+# one chunk read. Writing into such a chunk is a read-modify-write of the
+# whole chunk, so build_archive stages months on disk and writes one full
+# calendar year (all variables) at a time: ~32 rewrites of the store instead
+# of one per (variable, month).
 TIME_CHUNK_HOURS: int | None = None
 # Coarse cells per spatial chunk edge: a point query only ever decompresses
 # a SPATIAL_CHUNK_CELLS x SPATIAL_CHUNK_CELLS neighbourhood (20 km x 20 km
@@ -63,4 +111,5 @@ SPATIAL_CHUNK_CELLS = 2
 
 # --- Paths --------------------------------------------------------------
 DEFAULT_RAW_CACHE_DIR = Path("raw_cache")
+DEFAULT_STAGING_DIR = Path("staging")
 DEFAULT_ARCHIVE_PATH = Path("hostrada_coarse.zarr")

pyproject.toml

@@ -7,12 +7,12 @@ name = "hostrada-transpose"
 version = "0.1.0"
 description = "One-time ETL: rechunk HOSTRADA into a compact, time-major, spatially-aggregated archive for SimStadt."
 readme = "README.md"
-requires-python = ">=3.10,<3.14"
+requires-python = ">=3.11,<3.14"
 license = {text = "MIT"}
 dependencies = [
     "hostrada4py @ git+https://github.com/UdK-VPT/hostrada4py.git",
     "xarray",
-    "zarr",
+    "zarr>=3",
     "dask",
     "pandas",
     "numpy",