fix: partition regional land use retrieval

CHANGELOG.md

@@ -11,6 +11,7 @@
 
 - Generalized the proven Statbel population operator from a hardcoded Mol import to an approved geographic scope while keeping Mol as the backwards-compatible default.
 - Added one explicit regional synchronization command for official 2021-2025 population and 2013-2025 modern forest snapshots.
+- Added resumable municipality-partitioned WCS retrieval and native-resolution raster mosaicking after the upstream service rejected the complete regional response at its documented size limit.
 - Kept every fetch operator-triggered, idempotent and behind the canonical DatasetService upload path; no startup fetch, migration or API contract change was introduced.
 - Preserved separate Mol/regional series keys and honest partial-sector population and 10 m forest-area limitations.
 - Replaced internal provider identifiers with readable source labels in the primary map.

backend/README.md

@@ -1038,6 +1038,10 @@ existing immutable datasets are reused. Complete statistical sectors use exact
 published totals; a rectangle cutting a sector remains an area-weighted
 estimate. Forest area is measured within the official 10 m representation.
 Use `--fetch-only` to validate source artifacts without database mutation.
+The regional forest path partitions WCS requests by official municipality to
+stay within upstream response limits, then builds one retained 10 m mosaic and
+one normal regional vector Dataset. A failed source request leaves completed
+partition artifacts reusable and never lowers source resolution silently.
 
 ## Helpful repository scripts
 

backend/tests/test_sprint194_regional_timeseries.py

@@ -8,6 +8,10 @@ from pathlib import Path
 import sys
 import zipfile
 
+import numpy as np
+import rasterio
+from rasterio.transform import from_origin
+
 
 ROOT = Path(__file__).resolve().parents[2]
 SCRIPTS = ROOT / "scripts"
@@ -98,7 +102,8 @@ def test_regional_coordinator_builds_explicit_population_and_forest_commands(tmp
         max_landuse_features=500000,
     )
 
-    commands = dict(module.build_operator_commands(args, scope, tmp_path / "boundary.geojson"))
+    members_path = tmp_path / "municipalities.geojson"
+    commands = dict(module.build_operator_commands(args, scope, tmp_path / "boundary.geojson", members_path))
 
     assert set(commands) == {"population", "forest"}
     assert commands["population"][0] == sys.executable
@@ -107,6 +112,8 @@ def test_regional_coordinator_builds_explicit_population_and_forest_commands(tmp
     assert scope.project_name in commands["population"]
     assert commands["forest"][0] == sys.executable
     assert "--max-features" in commands["forest"]
+    assert "--partition-boundaries-path" in commands["forest"]
+    assert str(members_path) in commands["forest"]
     assert ",".join(scope.nis_codes) in commands["forest"]
     assert "--force" not in commands["population"]
     assert "--fetch-only" not in commands["forest"]
@@ -128,6 +135,39 @@ def test_regional_forest_provenance_does_not_claim_one_municipality() -> None:
     assert municipal["nis_code"] == "13025"
 
 
+def test_regional_forest_partition_rasters_merge_without_resolution_loss(tmp_path: Path) -> None:
+    module = load_script("provision_official_landuse_timeseries.py")
+    left = tmp_path / "left.tif"
+    right = tmp_path / "right.tif"
+    profile = {
+        "driver": "GTiff",
+        "height": 2,
+        "width": 2,
+        "count": 1,
+        "dtype": "uint8",
+        "crs": "EPSG:31370",
+        "transform": from_origin(100000, 200000, 10, 10),
+        "nodata": 0,
+    }
+    with rasterio.open(left, "w", **profile) as target:
+        target.write(np.full((1, 2, 2), 12, dtype="uint8"))
+    with rasterio.open(
+        right,
+        "w",
+        **{**profile, "transform": from_origin(100020, 200000, 10, 10)},
+    ) as target:
+        target.write(np.full((1, 2, 2), 17, dtype="uint8"))
+
+    destination = tmp_path / "regional.tif"
+    result = module.merge_partition_rasters([left, right], destination)
+
+    assert result["width"] == 4
+    assert result["height"] == 2
+    assert result["resolution_metres"] == 10.0
+    with rasterio.open(destination) as merged:
+        assert merged.read(1).tolist() == [[12, 12, 17, 17], [12, 12, 17, 17]]
+
+
 def test_regional_timeseries_operator_is_packaged_and_release_checked() -> None:
     dockerfile = (ROOT / "deploy/unraid/Dockerfile.all-in-one").read_text(encoding="utf-8")
     readiness = (ROOT / "scripts/run_readiness_check.sh").read_text(encoding="utf-8")

docs/CODEX_EXECUTION_LOG.md

@@ -7978,6 +7978,11 @@ Validation before deployment:
 Next:
 - Deploy the packaged operator, execute the live regional synchronization, verify persisted counts and temporal comparisons, then mark the regional time-series TODO complete only if PostGIS and browser evidence agree.
 
+Live source finding:
+- The regional Statbel synchronization imported all five requested snapshots successfully.
+- The first complete-region forest request was rejected by MercatorNet before download because the estimated 97.68 MB response exceeded its 78.12 MB service limit.
+- The operator now uses the 28 retained official municipality boundaries as resumable WCS partitions, merges them locally at the unchanged 10 m grid and applies the exact region-union clip before polygonization. It does not lower resolution or truncate the source.
+
 ## Sprint 190 Regional Kempen GRB buildings (2026-07-14)
 
 Implemented:

docs/DATA_SOURCES.md

@@ -163,6 +163,10 @@ The command is explicit and operator-triggered. No source fetch happens during
 startup or map interaction. Partial statistical sectors remain area-weighted
 population estimates; forest hectares remain measurements in the harmonized
 10 m source representation and are not cadastral forest boundaries.
+The WCS source limits large responses, so regional forest retrieval is
+partitioned by the same 28 official municipality boundaries. GeoIntel retains
+those source rasters, merges them on their native 10 m grid and applies the
+exact regional union clip before polygon persistence.
 
 Official catalogues:
 

scripts/README.md

@@ -1386,7 +1386,9 @@ coordinates Statbel 2021-2025 with Departement Omgeving
 canonical dataset upload API. It never runs at application startup. Prepare
 artifacts without persistence using `--fetch-only`; bound a run with
 `--skip-population`, `--skip-landuse`, `--population-years` or
-`--landuse-years`.
+`--landuse-years`. Regional WCS downloads use the 28 official municipality
+boundaries as resumable request partitions, preserve the native 10 m
+resolution and merge locally before exact clipping to the regional union.
 
 ## Tower deployment
 

scripts/provision_official_landuse_timeseries.py

@@ -28,6 +28,7 @@ import rasterio
 import requests
 from pyproj import Transformer
 from rasterio.features import geometry_mask, shapes
+from rasterio.merge import merge
 from requests.adapters import HTTPAdapter
 from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, mapping, shape
 from shapely.ops import transform, unary_union
@@ -118,6 +119,12 @@ def parse_args() -> argparse.Namespace:
         type=Path,
         default=Path(os.environ.get("OFFICIAL_LANDUSE_BOUNDARY_PATH", DEFAULT_BOUNDARY_PATH)),
     )
+    parser.add_argument(
+        "--partition-boundaries-path",
+        type=Path,
+        default=None,
+        help="Optional FeatureCollection of non-overlapping member boundaries used for bounded WCS requests.",
+    )
     parser.add_argument(
         "--output-dir",
         type=Path,
@@ -208,6 +215,30 @@ def load_boundary(path: Path):
     return boundary
 
 
+def load_partition_boundaries(path: Path | None) -> list[tuple[str, Any]]:
+    if path is None:
+        return []
+    if not path.is_file():
+        raise RuntimeError(f"Partition boundary artifact is missing at {path}")
+    payload = json.loads(path.read_text(encoding="utf-8"))
+    features = payload.get("features") or []
+    if not features:
+        raise RuntimeError("Partition boundary artifact contains no features")
+    partitions: list[tuple[str, Any]] = []
+    seen: set[str] = set()
+    for index, feature in enumerate(features):
+        properties = feature.get("properties") or {}
+        key = str(properties.get("nis_code") or properties.get("NISCODE") or feature.get("id") or index).strip()
+        if not key or key in seen:
+            raise RuntimeError(f"Partition boundary artifact contains an empty or duplicate key: {key!r}")
+        geometry = normalize_polygonal(shape(feature.get("geometry")))
+        if geometry is None:
+            raise RuntimeError(f"Partition boundary {key} is empty, invalid or non-polygonal")
+        seen.add(key)
+        partitions.append((key, geometry))
+    return partitions
+
+
 def normalize_polygonal(geometry):
     if geometry is None or geometry.is_empty:
         return None
@@ -325,6 +356,99 @@ def download_raster(
     return {**profile, "request_url": response.url, "retrieved_at": utc_now()}
 
 
+def merge_partition_rasters(paths: list[Path], destination: Path) -> dict[str, Any]:
+    if not paths:
+        raise RuntimeError("At least one partition raster is required for a regional mosaic")
+    sources = [rasterio.open(path) for path in paths]
+    temporary = destination.with_suffix(f"{destination.suffix}.partial")
+    try:
+        mosaic, transform_matrix = merge(sources)
+        profile = sources[0].profile.copy()
+        profile.pop("blockxsize", None)
+        profile.pop("blockysize", None)
+        tile_width = min(512, (mosaic.shape[2] // 16) * 16)
+        tile_height = min(512, (mosaic.shape[1] // 16) * 16)
+        profile.update(
+            driver="GTiff",
+            width=mosaic.shape[2],
+            height=mosaic.shape[1],
+            transform=transform_matrix,
+            count=mosaic.shape[0],
+            compress="deflate",
+            tiled=tile_width >= 16 and tile_height >= 16,
+            BIGTIFF="IF_SAFER",
+        )
+        if tile_width >= 16 and tile_height >= 16:
+            profile.update(blockxsize=tile_width, blockysize=tile_height)
+        with rasterio.open(temporary, "w", **profile) as target:
+            target.write(mosaic)
+        validated = validate_raster(temporary)
+        temporary.replace(destination)
+        return validated
+    except Exception:
+        temporary.unlink(missing_ok=True)
+        raise
+    finally:
+        for source in sources:
+            source.close()
+
+
+def download_partitioned_raster(
+    session: requests.Session,
+    *,
+    year: int,
+    partitions: list[tuple[str, Any]],
+    output_dir: Path,
+    destination: Path,
+    timeout: int,
+    force: bool,
+) -> dict[str, Any]:
+    partition_dir = output_dir / "source-partitions" / str(year)
+    partition_dir.mkdir(parents=True, exist_ok=True)
+    paths: list[Path] = []
+    details: list[dict[str, Any]] = []
+    for key, boundary in partitions:
+        path = partition_dir / f"land_use_{year}_{key}.tif"
+        boundary_metric = metric_boundary(boundary)
+        if force or not path.is_file():
+            profile = download_raster(
+                session,
+                year=year,
+                bounds=boundary_metric.bounds,
+                path=path,
+                timeout=timeout,
+            )
+            status = "downloaded"
+        else:
+            profile = validate_raster(path)
+            profile["request_url"] = requests.Request(
+                "GET", WCS_URL, params=build_wcs_params(year, boundary_metric.bounds)
+            ).prepare().url
+            profile["retrieved_at"] = None
+            status = "reused"
+        paths.append(path)
+        details.append(
+            {
+                "partition_key": key,
+                "path": str(path),
+                "sha256": sha256_file(path),
+                "status": status,
+                "request_url": profile.get("request_url"),
+                "width": profile["width"],
+                "height": profile["height"],
+            }
+        )
+
+    mosaic_profile = merge_partition_rasters(paths, destination)
+    return {
+        **mosaic_profile,
+        "request_mode": "partitioned_scope_members",
+        "partition_count": len(details),
+        "partitions": details,
+        "retrieved_at": utc_now(),
+    }
+
+
 def polygonize_snapshot(
     *,
     raster_path: Path,
@@ -492,6 +616,7 @@ def prepare_snapshot(
     args: argparse.Namespace,
     boundary,
     boundary_metric,
+    partition_boundaries: list[tuple[str, Any]],
     year: int,
     theme: ThemeDefinition,
 ) -> PreparedSnapshot:
@@ -512,13 +637,24 @@ def prepare_snapshot(
 
     raster_profile: dict[str, Any]
     if args.force or not raster_path.exists():
-        raster_profile = download_raster(
+        if partition_boundaries:
-            session,
+            raster_profile = download_partitioned_raster(
-            year=year,
+                session,
-            bounds=boundary_metric.bounds,
+                year=year,
-            path=raster_path,
+                partitions=partition_boundaries,
-            timeout=args.request_timeout,
+                output_dir=args.output_dir,
-        )
+                destination=raster_path,
+                timeout=args.request_timeout,
+                force=args.force,
+            )
+        else:
+            raster_profile = download_raster(
+                session,
+                year=year,
+                bounds=boundary_metric.bounds,
+                path=raster_path,
+                timeout=args.request_timeout,
+            )
     else:
         raster_profile = validate_raster(raster_path)
         raster_profile["request_url"] = requests.Request(
@@ -554,6 +690,11 @@ def prepare_snapshot(
         "output_crs": OUTPUT_CRS,
         "source_resolution_metres": SOURCE_RESOLUTION_METRES,
         "boundary_path": str(args.boundary_path),
+        "partition_boundaries_path": (
+            str(getattr(args, "partition_boundaries_path", None))
+            if getattr(args, "partition_boundaries_path", None)
+            else None
+        ),
         **identity,
         "scope_key": args.scope_key,
         "raster_path": str(raster_path),
@@ -681,6 +822,11 @@ def build_provenance_metadata(args: argparse.Namespace, snapshot: PreparedSnapsh
         "raw_raster_path": str(snapshot.raster_path),
         "polygon_artifact_path": str(snapshot.vector_path),
         "manifest_path": str(snapshot.manifest_path),
+        "partition_boundaries_path": (
+            str(getattr(args, "partition_boundaries_path", None))
+            if getattr(args, "partition_boundaries_path", None)
+            else None
+        ),
         "raster_sha256": snapshot.raster_sha256,
         "vector_sha256": snapshot.vector_sha256,
         "source_crs": SOURCE_CRS,
@@ -752,6 +898,7 @@ def main() -> int:
     try:
         boundary = load_boundary(args.boundary_path)
         boundary_metric = metric_boundary(boundary)
+        partition_boundaries = load_partition_boundaries(args.partition_boundaries_path)
         prepared: list[PreparedSnapshot] = []
         with build_session() as source_session:
             verify_coverages(source_session, years, args.request_timeout)
@@ -763,6 +910,7 @@ def main() -> int:
                             args=args,
                             boundary=boundary,
                             boundary_metric=boundary_metric,
+                            partition_boundaries=partition_boundaries,
                             year=year,
                             theme=theme,
                         )

scripts/provision_regional_timeseries.py

@@ -50,7 +50,7 @@ def parse_args() -> argparse.Namespace:
     return parser.parse_args()
 
 
-def resolve_boundary(scope: GeographicScope, scope_output_root: Path) -> tuple[Path, Path]:
+def resolve_boundary(scope: GeographicScope, scope_output_root: Path) -> tuple[Path, Path, Path]:
     scope_dir = scope_output_root / scope.key
     manifest_path = scope_dir / f"{scope.key.replace('-', '_')}_scope_manifest.json"
     if not manifest_path.is_file():
@@ -65,12 +65,20 @@ def resolve_boundary(scope: GeographicScope, scope_output_root: Path) -> tuple[P
     ):
         raise RuntimeError(f"Official scope manifest at {manifest_path} is incomplete or inconsistent")
     boundary_path = scope_dir / str(manifest.get("boundary_filename") or "")
+    members_path = scope_dir / str(manifest.get("municipalities_filename") or "")
     if not boundary_path.is_file():
         raise RuntimeError(f"Official scope boundary referenced by {manifest_path} is missing")
-    return boundary_path, manifest_path
+    if not members_path.is_file():
+        raise RuntimeError(f"Official scope member boundaries referenced by {manifest_path} are missing")
+    return boundary_path, members_path, manifest_path
 
 
-def build_operator_commands(args: argparse.Namespace, scope: GeographicScope, boundary_path: Path) -> list[tuple[str, list[str]]]:
+def build_operator_commands(
+    args: argparse.Namespace,
+    scope: GeographicScope,
+    boundary_path: Path,
+    members_path: Path | None = None,
+) -> list[tuple[str, list[str]]]:
     scripts_dir = Path(__file__).resolve().parent
     scope_output = args.output_root / scope.key
     common_flags = ["--fetch-only"] if args.fetch_only else []
@@ -108,6 +116,7 @@ def build_operator_commands(args: argparse.Namespace, scope: GeographicScope, bo
             )
         )
     if not args.skip_landuse:
+        partition_flags = ["--partition-boundaries-path", str(members_path)] if members_path else []
         commands.append(
             (
                 "forest",
@@ -132,6 +141,7 @@ def build_operator_commands(args: argparse.Namespace, scope: GeographicScope, bo
                     "forest",
                     "--boundary-path",
                     str(boundary_path),
+                    *partition_flags,
                     "--output-dir",
                     str(scope_output / "landuse"),
                     "--request-timeout",
@@ -169,8 +179,8 @@ def main() -> int:
     try:
         if args.max_landuse_features <= 0:
             raise ValueError("max-landuse-features must be greater than zero")
-        boundary_path, manifest_path = resolve_boundary(scope, args.scope_output_root)
+        boundary_path, members_path, manifest_path = resolve_boundary(scope, args.scope_output_root)
-        commands = build_operator_commands(args, scope, boundary_path)
+        commands = build_operator_commands(args, scope, boundary_path, members_path)
         results = {label: run_operator(label, command) for label, command in commands}
     except (OSError, RuntimeError, ValueError, KeyError) as exc:
         print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr)
@@ -185,6 +195,7 @@ def main() -> int:
                 "display_name": scope.display_name,
                 "member_count": len(scope.members),
                 "boundary_path": str(boundary_path),
+                "municipality_boundaries_path": str(members_path),
                 "scope_manifest_path": str(manifest_path),
                 "results": results,
             },