Python project to download, transpose and save HOSTRADA data

.gitignore

+.venv/
+__pycache__/
+*.pyc
+raw_cache/
+*.zarr/
+.DS_Store

README.md

+# hostrada-transpose
+
+One-time ETL that turns HOSTRADA's national, map-oriented NetCDF files into a
+small, time-major archive suited to point queries -- built to feed a future
+SimStadt HOSTRADA integration (as an alternative to Meteonorm) without
+making every workflow run download tens of GB from DWD.
+
+Depends on [`hostrada4py`](https://github.com/UdK-VPT/hostrada4py) for all
+download/caching/provider logic; this repo only adds the aggregation and
+rechunking step hostrada4py doesn't do itself.
+
+## Why this exists
+
+HOSTRADA's monthly NetCDF files are chunked as `(1, 938, 720)` -- one
+**entire national grid** per hour, gzip-compressed (verified via `h5py` on a
+real cached file). There is no byte range that gives you "just one point":
+extracting a single pixel's time series still requires decompressing every
+hourly chunk, i.e. close to 100% of the file. HTTP-range subsetting
+(`HOSTRADA_NETCDF_SUBSET_MODE=http_range` in hostrada4py) does not help here
+-- it only reduces what you ingest if the source is chunked in a way that
+matches your query, which this isn't.
+
+For a handful of SimStadt cities, downloading full national grids per
+request is wasteful. This repo inverts the layout **once**: aggregate the
+1 km grid to a coarser cell size and store the full multi-decade time series
+contiguously per cell, so a later point query is a cheap local array read
+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).
+
+Run it with:
+
+```bash
+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
+```
+
+Reading one point's full time series back out looks like:
+
+```python
+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")
+```
+
+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.
+
+## 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. |
+
+**Precipitation is deliberately out of scope here.** It's required by
+GreenWater but absent from both hostrada4py providers: not in DWD
+HOSTRADA's `BASE_URLS`, and not mapped in CERRA's `CERRA_VARIABLES` even
+though CERRA's upstream CDS catalogue has `total_precipitation`. Adding it
+belongs in `hostrada4py` itself (other consumers would want it too), as a
+small addition to `providers/cerra.py`. The agreed fallback resolution is
+CERRA's native ~5.5 km grid, not aggregated HOSTRADA (HOSTRADA has nothing
+to aggregate).
+
+## DHI: modelled vs. measured
+
+`HostradaDiffuse` only ever **estimates** DHI from GHI via the
+Erbs/Erbs-Driesse decomposition (solar geometry + an empirical cloudiness
+correlation -- no actual diffuse-radiation measurement involved). DWD
+separately publishes real pyranometer measurements:
+
+> "10-minute station observations of solar and sunshine for Germany"
+> `https://opendata.dwd.de/climate_environment/CDC/observations_germany/climate/10_minutes/solar/`
+> columns `STATIONS_ID;MESS_DATUM;QN;DS_10;GS_10;SD_10;LS_10;eor`
+> (`DS_10` = diffuse, `GS_10` = global, J/cm² per 10-min sum, -999 = missing)
+> **339 stations** (verified by parsing the live station list, not the ~500+
+> a first automated skim suggested), 1989/1991-present, CC BY 4.0.
+
+`hostrada_transpose.dhi_validation` fetches this and joins it against
+hostrada4py's modelled DHI at the same coordinates, as **ground truth to
+quantify the Erbs-Driesse model's error**, not as a gridded replacement
+(339 unevenly-spaced points is too sparse for that).
+
+```python
+from pathlib import Path
+import pandas as pd
+from hostrada_transpose.dhi_validation import fetch_station_list, compare_station_to_model
+
+stations = fetch_station_list()
+station = next(s for s in stations if s.station_id == "00867")  # Lautertal-Oberlauter
+result = compare_station_to_model(
+    station, pd.Timestamp("2024-06-01", tz="UTC"), pd.Timestamp("2024-06-07", tz="UTC"),
+    cache_dir=Path("./raw_cache"),
+)
+print(result["dhi_bias_wm2"].mean(), (result["dhi_bias_wm2"] ** 2).mean() ** 0.5)
+```
+
+**Verified working end-to-end** against real data (station 00867,
+Lautertal-Oberlauter, Bayern, one week of June 2024): mean bias +3.7 W/m²,
+RMSE 48.5 W/m² (59 W/m² restricted to daylight hours). **This is one
+station, one week -- illustrative that the pipeline works, not a validated
+error figure for the dataset.** A real error characterization needs many
+stations across multiple years and seasons; that aggregation is not
+implemented yet.
+
+**Known gotcha, found while testing:** `HostradaDiffuse.estimate()` does
+`ghi.fillna(0.0)` internally. A station whose nearest HOSTRADA cell is NaN
+(coastline/border -- roughly half the national bounding box has no land
+data) silently produces `dhi=0` all day instead of a visible error, which
+would masquerade as a huge bogus "bias" rather than a missing-data problem.
+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)
+
+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.
+
+## 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.

hostrada_transpose/__init__.py

+"""hostrada-transpose: build a compact, pre-aggregated HOSTRADA archive for SimStadt.
+
+See README.md for the design rationale (why this exists, what it does not
+try to be, and the open questions still worth revisiting).
+"""

hostrada_transpose/aggregate.py

+"""Spatial block-aggregation of one HOSTRADA monthly grid file."""
+from __future__ import annotations
+
+import xarray as xr
+from hostrada4py.hostrada import find_variable
+
+from .config import AGGREGATION_FACTOR
+
+
+def aggregate_month(
+    ds: xr.Dataset,
+    var: str,
+    factor: int = AGGREGATION_FACTOR,
+) -> xr.DataArray:
+    """Block-average one month of a national HOSTRADA grid to a coarser grid.
+
+    Uses a plain arithmetic mean: verified that HOSTRADA cells are a uniform
+    1 km x 1 km grid in EPSG:3034, so no area weighting is needed.
+    ``boundary="trim"`` drops the partial row/column at the domain edge
+    rather than padding with NaN. A block's mean is NaN only if every
+    1 km cell inside it is NaN (xarray's default skipna=True) -- this is
+    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).
+    """
+    name = find_variable(var, ds)
+    da = ds[name]
+    coarse = da.coarsen(X=factor, Y=factor, boundary="trim").mean()
+    coarse.name = var
+    coarse.attrs = dict(da.attrs)
+    # HOSTRADA attaches 2-D lon/lat auxiliary coords (dims Y, X). coarsen()
+    # block-averages them along with the data, which is both numerically
+    # sloppy for lon/lat and a mismatch for the archive template (which only
+    # has time/Y/X). Coarse-cell centre lon/lat can be computed once from
+    # X/Y via pyproj if ever needed; they aren't carried through here.
+    extra_coords = [c for c in coarse.coords if c not in ("time", "Y", "X")]
+    if extra_coords:
+        coarse = coarse.drop_vars(extra_coords)
+    return coarse

hostrada_transpose/archive.py

+"""Create and incrementally 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.
+
+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.
+"""
+from __future__ import annotations
+
+from pathlib import Path
+from typing import Mapping, Sequence
+
+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
+
+ARRAY_NAME = "weather"
+
+
+def build_time_index(start: str, end: str) -> pd.DatetimeIndex:
+    """Full hourly time axis for the archive, inclusive of the whole end month."""
+    end_ts = pd.Timestamp(end) + pd.offsets.MonthEnd(0) + pd.Timedelta(hours=23)
+    return pd.date_range(start=pd.Timestamp(start), end=end_ts, freq="h")
+
+
+def create_archive_template(
+    path: Path,
+    variables: Sequence[str],
+    time_index: pd.DatetimeIndex,
+    coarse_coords: Mapping[str, np.ndarray],
+) -> None:
+    """Write an empty archive: only metadata and chunk layout touch disk.
+
+    ``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`.
+    """
+    n_time = len(time_index)
+    y_vals = coarse_coords["Y"]
+    x_vals = coarse_coords["X"]
+    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),
+    )
+    template = xr.Dataset(
+        {ARRAY_NAME: (("time", "Y", "X", "variable"), placeholder)},
+        coords={"time": time_index, "Y": y_vals, "X": x_vals, "variable": list(variables)},
+    )
+    template.to_zarr(path, mode="w", compute=False)
+
+
+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.
+
+    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),
+        },
+    )

hostrada_transpose/build_archive.py

+"""CLI entry point: build the coarse, time-major HOSTRADA archive.
+
+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
+
+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.
+"""
+from __future__ import annotations
+
+import argparse
+import logging
+from pathlib import Path
+
+import hostrada4py.hostrada as hs
+
+from .aggregate import aggregate_month
+from .archive import build_time_index, create_archive_template, write_month
+from .config import (
+    ALL_VARIABLES,
+    DEFAULT_ARCHIVE_PATH,
+    DEFAULT_END,
+    DEFAULT_RAW_CACHE_DIR,
+    DEFAULT_START,
+)
+
+logger = logging.getLogger("hostrada_transpose")
+
+
+def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    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(
+        "--variables",
+        default=",".join(ALL_VARIABLES),
+        help="Comma-separated canonical HOSTRADA variable names",
+    )
+    parser.add_argument("--raw-cache-dir", type=Path, default=DEFAULT_RAW_CACHE_DIR)
+    parser.add_argument("--archive-path", type=Path, default=DEFAULT_ARCHIVE_PATH)
+    parser.add_argument(
+        "--keep-raw",
+        action="store_true",
+        help="Keep each month's full national-grid download after aggregating it "
+        "(default: delete it, since only the coarse archive is meant to persist).",
+    )
+    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.",
+    )
+    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
+
+    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 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()]
+
+    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)
+
+    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.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)
+                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)
+
+
+if __name__ == "__main__":
+    main()

hostrada_transpose/config.py

+from __future__ import annotations
+
+from pathlib import Path
+
+# --- Variables -----------------------------------------------------------
+# 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
+
+# Precipitation is intentionally NOT handled here. Neither the DWD HOSTRADA
+# provider nor hostrada4py's CERRA provider currently expose it (confirmed by
+# grepping providers/dwd_hostrada.py and providers/cerra.py -- CERRA's own
+# upstream CDS catalogue has `total_precipitation`, it's just unmapped in
+# CERRA_VARIABLES). Adding it belongs in hostrada4py itself, not here, since
+# it's provider-level plumbing other hostrada4py consumers would want too.
+
+# --- 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
+# arithmetic block mean (no area weighting needed: cells are equal-area in
+# the source grid). boundary="trim" drops the partial last row/column
+# (<1 coarse cell) rather than padding with NaN.
+AGGREGATION_FACTOR = 10
+
+# --- Time range -------------------------------------------------------
+# 1995-01 is HOSTRADA's own archive start. DEFAULT_END is kept close to the
+# live edge of DWD's publication (verified: tas/rsds for August 2026 were
+# already published, with roughly a one-month lag behind the current date)
+# rather than an arbitrary round year -- update it periodically if this is
+# re-run later to pick up newer months.
+DEFAULT_START = "1995-01"
+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.
+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
+# at the default of 2), not the whole country.
+SPATIAL_CHUNK_CELLS = 2
+
+# --- Paths --------------------------------------------------------------
+DEFAULT_RAW_CACHE_DIR = Path("raw_cache")
+DEFAULT_ARCHIVE_PATH = Path("hostrada_coarse.zarr")

hostrada_transpose/dhi_validation.py

+"""Validate modelled DHI against DWD's measured diffuse-radiation stations.
+
+hostrada4py's ``HostradaDiffuse`` only ever *estimates* DHI from GHI via the
+Erbs/Erbs-Driesse decomposition (solar-geometry + empirical cloudiness
+correlation, no actual diffuse measurement involved). DWD separately
+publishes real pyranometer measurements of diffuse sky radiation:
+
+    "10-minute station observations of solar and sunshine for Germany"
+    https://opendata.dwd.de/climate_environment/CDC/observations_germany/climate/10_minutes/solar/
+    columns: STATIONS_ID; MESS_DATUM; QN; DS_10 (diffuse); GS_10 (global);
+             SD_10 (sunshine); LS_10 (longwave); eor
+    units: J/cm^2 per 10-minute sum; -999 = missing. CC BY 4.0.
+
+Verified directly against a downloaded sample
+(10minutenwerte_SOLAR_00183_20200101_20251231_hist.zip) rather than assumed
+from the dataset's PDF description alone, which only lists the physical
+parameters, not the exact file layout.
+
+This module compares that measured DHI against the HOSTRADA-derived
+estimate at the same station locations, to get an actual documented error
+figure instead of trusting the decomposition model blindly.
+"""
+from __future__ import annotations
+
+import io
+import re
+import zipfile
+from dataclasses import dataclass
+from pathlib import Path
+
+import pandas as pd
+import requests
+
+STATION_LIST_URL = (
+    "https://opendata.dwd.de/climate_environment/CDC/observations_germany/"
+    "climate/10_minutes/solar/historical/zehn_min_sd_Beschreibung_Stationen.txt"
+)
+HISTORICAL_INDEX_URL = (
+    "https://opendata.dwd.de/climate_environment/CDC/observations_germany/"
+    "climate/10_minutes/solar/historical/"
+)
+ZIP_NAME_RE = re.compile(r"10minutenwerte_SOLAR_(\d{5})_(\d{8})_(\d{8})_hist\.zip")
+
+# J/cm^2 per 10-minute sum -> average W/m^2 over that interval:
+# 1 J/cm^2 = 1e4 J/m^2; divide by 600 s -> W/m^2.
+J_PER_CM2_TO_W_PER_M2_OVER_10MIN = 1e4 / 600
+
+
+@dataclass(frozen=True)
+class Station:
+    station_id: str  # zero-padded to 5 digits, matches the zip filename
+    lat: float
+    lon: float
+    height_m: float
+    name: str
+    state: str
+    start: pd.Timestamp
+    end: pd.Timestamp
+
+
+# All 16 German states -- used as an anchor to split the free-text station
+# description line, since it is NOT reliably fixed-width despite the dashed
+# header (verified against the live file: numeric field widths vary and
+# don't line up with the dashes). Station names never contain whitespace in
+# the current list (verified against all 339 stations), so
+# ``tokens[6:-2]`` joined is always the full station name.
+_BUNDESLAENDER = {
+    "Baden-Württemberg", "Bayern", "Berlin", "Brandenburg", "Bremen", "Hamburg",
+    "Hessen", "Mecklenburg-Vorpommern", "Niedersachsen", "Nordrhein-Westfalen",
+    "Rheinland-Pfalz", "Saarland", "Sachsen", "Sachsen-Anhalt",
+    "Schleswig-Holstein", "Thüringen",
+}
+
+
+def fetch_station_list(timeout: float = 30.0) -> list[Station]:
+    """Download and parse the station description file (339 stations as of writing).
+
+    Encoding is ISO-8859-1 (verified: umlauts are mis-decoded as UTF-8
+    otherwise). Parsed by token position anchored on the trailing
+    Bundesland/Abgabe columns rather than the dashed header widths, which do
+    not actually match the data's column boundaries.
+    """
+    resp = requests.get(STATION_LIST_URL, timeout=timeout)
+    resp.raise_for_status()
+    text = resp.content.decode("latin-1")
+    stations = []
+    for line in text.splitlines()[2:]:  # skip header + dashed separator
+        tokens = line.split()
+        if len(tokens) < 8 or tokens[-2] not in _BUNDESLAENDER:
+            continue
+        station_id, von_datum, bis_datum, height, lat, lon, *rest = tokens
+        name = " ".join(rest[:-2])
+        bundesland, abgabe = rest[-2], rest[-1]
+        stations.append(
+            Station(
+                station_id=f"{int(station_id):05d}",
+                lat=float(lat),
+                lon=float(lon),
+                height_m=float(height),
+                name=name,
+                state=bundesland,
+                start=pd.Timestamp(von_datum, tz="UTC"),
+                end=pd.Timestamp(bis_datum, tz="UTC"),
+            )
+        )
+    return stations
+
+
+def _list_station_zip_urls(station_id: str, timeout: float = 30.0) -> list[tuple[pd.Timestamp, pd.Timestamp, str]]:
+    """Find every historical zip URL for one station via the directory index.
+
+    There is no listing API -- DWD serves a plain Apache-style index, so this
+    greps hrefs out of the HTML. Each station typically has several
+    non-overlapping decade-ish zips, not one file per station.
+    """
+    resp = requests.get(HISTORICAL_INDEX_URL, timeout=timeout)
+    resp.raise_for_status()
+    matches = []
+    for m in ZIP_NAME_RE.finditer(resp.text):
+        filename_station, begin, end = m.groups()
+        if filename_station != station_id:
+            continue
+        url = HISTORICAL_INDEX_URL + m.group(0)
+        matches.append((pd.Timestamp(begin, tz="UTC"), pd.Timestamp(end, tz="UTC"), url))
+    return matches
+
+
+def fetch_station_timeseries(
+    station_id: str,
+    start: pd.Timestamp,
+    end: pd.Timestamp,
+    timeout: float = 120.0,
+) -> pd.DataFrame:
+    """Download, extract and parse one station's measured DS_10/GS_10.
+
+    Returns an hourly-mean DataFrame (not the raw 10-minute values) with
+    columns ``dhi_measured_wm2`` and ``ghi_measured_wm2``, to match the
+    hourly cadence hostrada4py/HostradaDiffuse operates at.
+    """
+    frames = []
+    for file_start, file_end, url in _list_station_zip_urls(station_id, timeout=timeout):
+        if file_end < start or file_start > end:
+            continue
+        resp = requests.get(url, timeout=timeout)
+        resp.raise_for_status()
+        with zipfile.ZipFile(io.BytesIO(resp.content)) as zf:
+            names = [n for n in zf.namelist() if n.startswith("produkt")]
+            if not names:
+                raise ValueError(f"No 'produkt*' data file found inside {url}")
+            with zf.open(names[0]) as fh:
+                df = pd.read_csv(fh, sep=";")
+        df.columns = [c.strip() for c in df.columns]
+        df["time"] = pd.to_datetime(df["MESS_DATUM"], format="%Y%m%d%H%M", utc=True)
+        for col in ("DS_10", "GS_10"):
+            df[col] = pd.to_numeric(df[col], errors="coerce")
+            df.loc[df[col] <= -999, col] = pd.NA
+        frames.append(df[["time", "DS_10", "GS_10"]])
+
+    if not frames:
+        return pd.DataFrame(columns=["time", "dhi_measured_wm2", "ghi_measured_wm2"]).set_index("time")
+
+    combined = pd.concat(frames, ignore_index=True).drop_duplicates("time").sort_values("time")
+    combined = combined.set_index("time").loc[start:end]
+    hourly = combined.resample("1h").mean() * J_PER_CM2_TO_W_PER_M2_OVER_10MIN
+    hourly = hourly.rename(columns={"DS_10": "dhi_measured_wm2", "GS_10": "ghi_measured_wm2"})
+    return hourly
+
+
+def compare_station_to_model(
+    station: Station,
+    start: pd.Timestamp,
+    end: pd.Timestamp,
+    cache_dir: Path,
+    apply_weather_correction: bool = False,
+) -> pd.DataFrame:
+    """Join measured DHI/GHI against hostrada4py's Erbs-Driesse DHI estimate.
+
+    Returns an hourly DataFrame with both series plus ``dhi_bias_wm2``
+    (modelled - measured), so callers can aggregate bias/RMSE by month,
+    season or sky condition as needed. Requires network access to both DWD
+    endpoints and hostrada4py's configured provider.
+    """
+    from hostrada4py.hostradaPoint import extract_diffuse_radiation_for_point
+
+    measured = fetch_station_timeseries(station.station_id, start, end)
+    modelled = extract_diffuse_radiation_for_point(
+        lon=station.lon,
+        lat=station.lat,
+        start=start,
+        end=end,
+        apply_weather_correction=apply_weather_correction,
+        cache_dir=cache_dir,
+    )
+    # HostradaDiffuse.estimate() does ghi.fillna(0.0) internally -- a station
+    # whose nearest HOSTRADA cell is NaN (coastline/border, where roughly
+    # half the national grid has no land data -- verified) silently yields
+    # dhi=0 all day, which would masquerade as a huge, bogus "bias" rather
+    # than a missing-data problem. Confirmed with station 00183 (Arkona, a
+    # lighthouse on the Baltic coast): rsds came back all-NaN at that point.
+    nan_fraction = modelled["rsds"].isna().mean()
+    if nan_fraction > 0:
+        raise ValueError(
+            f"{station.name} ({station.lat}, {station.lon}): HOSTRADA rsds is NaN for "
+            f"{nan_fraction:.0%} of the requested period at the nearest grid cell -- "
+            "likely too close to the coastline/border. Pick a more interior station."
+        )
+    modelled = modelled.set_index(pd.DatetimeIndex(pd.to_datetime(modelled["time"], utc=True)))
+    joined = measured.join(modelled[["dhi"]].rename(columns={"dhi": "dhi_modelled_wm2"}), how="inner")
+    joined["dhi_bias_wm2"] = joined["dhi_modelled_wm2"] - joined["dhi_measured_wm2"]
+    return joined

pyproject.toml

+[build-system]
+requires = ["setuptools>=68"]
+build-backend = "setuptools.build_meta"
+
+[project]
+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"
+license = {text = "MIT"}
+dependencies = [
+    "hostrada4py @ git+https://github.com/UdK-VPT/hostrada4py.git",
+    "xarray",
+    "zarr",
+    "dask",
+    "pandas",
+    "numpy",
+    "requests",
+    "pvlib",
+    "netcdf4",
+]
+
+[project.optional-dependencies]
+test = ["pytest>=8"]
+
+[tool.setuptools.packages.find]
+include = ["hostrada_transpose*"]