diff --git a/CHANGELOG.md b/CHANGELOG.md index 690f8127..86e233a5 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/backend/README.md b/backend/README.md index 6518cc6c..35835bd2 100644 --- a/backend/README.md +++ b/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 diff --git a/backend/tests/test_sprint194_regional_timeseries.py b/backend/tests/test_sprint194_regional_timeseries.py index d89cffbd..78c13e35 100644 --- a/backend/tests/test_sprint194_regional_timeseries.py +++ b/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") diff --git a/docs/CODEX_EXECUTION_LOG.md b/docs/CODEX_EXECUTION_LOG.md index 8e8324c5..65bd7070 100644 --- a/docs/CODEX_EXECUTION_LOG.md +++ b/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: diff --git a/docs/DATA_SOURCES.md b/docs/DATA_SOURCES.md index 457cc02b..de65a38b 100644 --- a/docs/DATA_SOURCES.md +++ b/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: diff --git a/scripts/README.md b/scripts/README.md index 53bf89ea..5557dafc 100644 --- a/scripts/README.md +++ b/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 diff --git a/scripts/provision_official_landuse_timeseries.py b/scripts/provision_official_landuse_timeseries.py index baa2bde5..cd152b94 100644 --- a/scripts/provision_official_landuse_timeseries.py +++ b/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( - session, - year=year, - bounds=boundary_metric.bounds, - path=raster_path, - timeout=args.request_timeout, - ) + if partition_boundaries: + raster_profile = download_partitioned_raster( + session, + year=year, + partitions=partition_boundaries, + 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, ) diff --git a/scripts/provision_regional_timeseries.py b/scripts/provision_regional_timeseries.py index 641fd850..3fd6bc5d 100644 --- a/scripts/provision_regional_timeseries.py +++ b/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) - commands = build_operator_commands(args, scope, boundary_path) + boundary_path, members_path, manifest_path = resolve_boundary(scope, args.scope_output_root) + 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, },