feat: add official modern Mol land-use series

CHANGELOG.md

@@ -7,6 +7,15 @@
 
 # Changelog
 
+## Sprint 188 Official modern Mol land-use series (2026-07-14)
+
+- Added an explicit, reusable MercatorNet WCS operator for official Departement Omgeving land-use snapshots in 2013, 2016, 2019, 2022 and 2025.
+- Validated categorical integer GeoTIFF input at 10 m in EPSG:31370, retained raw rasters/checksums/manifests and polygonized only documented class 12 (`Bos`).
+- Imported forest polygons through the existing DatasetService/VectorFeatureService path with canonical EPSG:4326 geometry and hectare intersection metrics; no startup fetch or direct PostGIS write was added.
+- Kept the modern 2013-2025 series separate from historical 1778/1873/1969 cartography and added a compact temporal-series selector when both are available.
+- Made the map-first forest theme prefer the authoritative modern source and its latest 2025 snapshot while preserving historical access.
+- Added focused raster, CRS, class, provenance, packaging and frontend contract tests plus readiness compilation coverage.
+
 ## Sprint 187 Temporal Mol explorer (2026-07-14)
 
 - Added first-class temporal dataset metadata and immutable dataset-version provenance for uploaded and derived vector/raster datasets.

backend/README.md

@@ -947,12 +947,15 @@ snapshots explicitly:
 ```bash
 docker exec geointel python /app/scripts/provision_mol_population_history.py
 docker exec geointel python /app/scripts/provision_mol_historical_landuse.py
+docker exec geointel python /app/scripts/provision_official_landuse_timeseries.py
 ```
 
 The first command imports Statbel sector population for 2021-2025. The second
-imports Digitaal Vlaanderen historical land use for 1778, 1873 and 1969. Both
+imports Digitaal Vlaanderen historical land use for 1778, 1873 and 1969. The
-are idempotent, use the normal API/DatasetService flow and retain fetched
+third imports the Departement Omgeving 10 m forest class for 2013, 2016, 2019,
-artifacts in persistent operator storage. They never run on app startup.
+2022 and 2025. All commands are idempotent, use the normal
+API/DatasetService flow and retain fetched artifacts in persistent operator
+storage. They never run on app startup.
 
 Historical land-use work can be bounded explicitly:
 
@@ -964,6 +967,9 @@ docker exec geointel python /app/scripts/provision_mol_historical_landuse.py --y
 `POST /api/v1/projects/{project_id}/temporal/compare` compares two snapshots
 inside one EPSG:4326 bbox. Partial statistical sectors are estimates; old map
 editions without stable identities do not produce invented object changes.
+Modern raster-derived forest polygons have the same identity limitation. Their
+area is measured in EPSG:31370 and is exact within the 10 m source
+representation, not a cadastral forest survey.
 
 ## Helpful repository scripts
 

backend/tests/test_sprint188_official_landuse_timeseries.py

@@ -0,0 +1,155 @@
+from __future__ import annotations
+
+import importlib.util
+from pathlib import Path
+from types import SimpleNamespace
+import sys
+
+import numpy as np
+from pyproj import Transformer
+import rasterio
+from rasterio.transform import from_origin
+from shapely.geometry import Polygon, shape
+from shapely.ops import transform
+
+
+ROOT = Path(__file__).resolve().parents[2]
+
+
+def load_provisioner():
+    script_path = ROOT / "scripts" / "provision_official_landuse_timeseries.py"
+    spec = importlib.util.spec_from_file_location("official_landuse_provisioner", script_path)
+    assert spec is not None
+    assert spec.loader is not None
+    module = importlib.util.module_from_spec(spec)
+    sys.modules[spec.name] = module
+    spec.loader.exec_module(module)
+    return module
+
+
+def test_official_landuse_wcs_contract_is_categorical_and_deterministic() -> None:
+    module = load_provisioner()
+
+    params = module.build_wcs_params(2025, (196594.1, 205064.2, 210910.7, 223974.9))
+
+    assert module.SUPPORTED_YEARS == (2013, 2016, 2019, 2022, 2025)
+    assert module.LAND_USE_CLASSES[12] == "Bos"
+    assert params == {
+        "SERVICE": "WCS",
+        "VERSION": "1.0.0",
+        "REQUEST": "GetCoverage",
+        "COVERAGE": "lu:lu_landgebruik_vlaa_2025_v3",
+        "CRS": "EPSG:31370",
+        "BBOX": "196590.000,205060.000,210920.000,223980.000",
+        "RESX": "10",
+        "RESY": "10",
+        "FORMAT": "image/tiff",
+        "RESPONSE_CRS": "EPSG:31370",
+    }
+    assert module.series_key(module.THEMES[0], "Mol") == "department-omgeving:land-use:forest:mol"
+
+
+def test_official_landuse_polygonization_clips_and_preserves_provenance(tmp_path: Path) -> None:
+    module = load_provisioner()
+    raster_path = tmp_path / "landuse.tif"
+    values = np.array(
+        [
+            [1, 1, 1, 1, 1, 1],
+            [1, 12, 12, 1, 1, 1],
+            [1, 12, 12, 1, 12, 1],
+            [1, 1, 1, 1, 12, 1],
+            [1, 1, 1, 1, 1, 1],
+            [1, 1, 1, 1, 1, 1],
+        ],
+        dtype="int32",
+    )
+    with rasterio.open(
+        raster_path,
+        "w",
+        driver="GTiff",
+        width=6,
+        height=6,
+        count=1,
+        dtype="int32",
+        crs="EPSG:31370",
+        transform=from_origin(200000, 210000, 10, 10),
+        nodata=-9999,
+    ) as destination:
+        destination.write(values, 1)
+
+    to_wgs84 = Transformer.from_crs(31370, 4326, always_xy=True)
+    boundary_metric = Polygon(
+        [(200005, 209945), (200055, 209945), (200055, 209995), (200005, 209995), (200005, 209945)]
+    )
+    boundary = transform(to_wgs84.transform, boundary_metric)
+
+    payload, stats = module.polygonize_snapshot(
+        raster_path=raster_path,
+        boundary=boundary,
+        year=2025,
+        theme=module.THEMES[0],
+        municipality_name="Mol",
+        nis_code="13025",
+        scope_key="mol",
+        max_features=100,
+    )
+
+    assert payload["type"] == "FeatureCollection"
+    assert payload["crs"]["properties"]["name"] == "EPSG:4326"
+    assert payload["source_coverage_id"] == "lu:lu_landgebruik_vlaa_2025_v3"
+    assert stats["source_pixel_count"] == 6
+    assert stats["source_pixel_area_m2"] == 600
+    assert stats["feature_count"] == 2
+    assert 0 < stats["polygon_area_m2"] <= 600
+    assert stats["class_histogram"] == {"1": 30, "12": 6}
+    for feature in payload["features"]:
+        geometry = shape(feature["geometry"])
+        properties = feature["properties"]
+        assert geometry.is_valid
+        assert geometry.within(boundary.buffer(1e-9))
+        assert properties["source_name"] == "department_omgeving_land_use"
+        assert properties["land_use_class_ids"] == [12]
+        assert properties["source_resolution_m"] == 10.0
+        assert properties["source_raster_sha256"] == stats["raster_sha256"]
+
+
+def test_official_landuse_metadata_keeps_modern_series_separate(tmp_path: Path) -> None:
+    module = load_provisioner()
+    theme = module.THEMES[0]
+    snapshot = module.PreparedSnapshot(
+        year=2022,
+        theme=theme,
+        raster_path=tmp_path / "source.tif",
+        vector_path=tmp_path / "forest.geojson",
+        manifest_path=tmp_path / "forest.manifest.json",
+        feature_count=42,
+        raster_sha256="a" * 64,
+        vector_sha256="b" * 64,
+    )
+    args = SimpleNamespace(scope_key="mol", municipality_name="Mol", nis_code="13025")
+
+    source = module.build_source_metadata(args, snapshot)
+    provenance = module.build_provenance_metadata(args, snapshot)
+
+    assert source["temporal_series_label"] == "Moderne landgebruikskaart (10 m)"
+    assert source["selection_aggregation"]["method"] == "intersection_area"
+    assert source["identity_stable"] is False
+    assert source["land_use_class_names"] == ["Bos"]
+    assert "10 m" in source["selection_aggregation"]["warning"]
+    assert provenance["operator_explicit_fetch"] is True
+    assert provenance["coverage_id"] == "lu:lu_landgebruik_vlaa_2022_v3"
+    assert "historical-landuse" not in module.series_key(theme, "mol")
+
+
+def test_official_landuse_operator_is_packaged_and_readiness_checked() -> None:
+    readiness = (ROOT / "scripts/run_readiness_check.sh").read_text(encoding="utf-8")
+    dockerfile = (ROOT / "deploy/unraid/Dockerfile.all-in-one").read_text(encoding="utf-8")
+    workspace = (ROOT / "frontend/src/components/map/MapWorkspace.tsx").read_text(encoding="utf-8")
+
+    assert "py_compile scripts/provision_official_landuse_timeseries.py" in readiness
+    assert "COPY scripts/provision_official_landuse_timeseries.py" in dockerfile
+    assert "department_omgeving_land_use' ? 90_000" in workspace
+    assert "activeTemporalSeriesGroups.length > 1" in workspace
+    assert "Moderne landgebruikskaart (10 m)" in (ROOT / "scripts/provision_official_landuse_timeseries.py").read_text(
+        encoding="utf-8"
+    )

deploy/unraid/Dockerfile.all-in-one

@@ -76,6 +76,7 @@ COPY scripts/provision_mol_municipality_workspace.py /app/scripts/provision_mol_
 COPY scripts/provision_mol_context_layers.py /app/scripts/provision_mol_context_layers.py
 COPY scripts/provision_mol_population_history.py /app/scripts/provision_mol_population_history.py
 COPY scripts/provision_mol_historical_landuse.py /app/scripts/provision_mol_historical_landuse.py
+COPY scripts/provision_official_landuse_timeseries.py /app/scripts/provision_official_landuse_timeseries.py
 COPY scripts/export_operator_yolo_tile_dataset.py /app/scripts/export_operator_yolo_tile_dataset.py
 COPY scripts/audit_operator_yolo_dataset_quality.py /app/scripts/audit_operator_yolo_dataset_quality.py
 COPY scripts/render_operator_yolo_label_qa_contact_sheets.py /app/scripts/render_operator_yolo_label_qa_contact_sheets.py

docs/CODEX_EXECUTION_LOG.md

@@ -7819,3 +7819,27 @@ Remaining source limitations:
 - The latest official population snapshot in this workspace is 2025.
 - The forest/green evolution series currently represents the available historical land-use editions through 1969; it must not be presented as current forest cover.
 - Partial-sector population results remain area-weighted estimates because no finer authoritative population surface has been ingested.
+
+## Sprint 188 Official modern Mol land-use series (2026-07-14)
+
+Implemented:
+- Added `provision_official_landuse_timeseries.py` for explicit MercatorNet WCS subsets of the official Departement Omgeving version 3 land-use maps for 2013, 2016, 2019, 2022 and 2025.
+- Validated one-band integer GeoTIFF input, EPSG:31370, 10 m resolution and the documented 1-19 class domain before processing.
+- Preserved every raw raster, exact request URL, catalogue URL, raster/vector checksum, class histogram and processing manifest in operator storage.
+- Polygonized documented class 12 (`Bos`) in metric CRS, clipped it to the official Mol boundary, normalized it to EPSG:4326 and uploaded it through the existing API/DatasetService/vector-feature path.
+- Added source-governed hectare aggregation metadata, 10 m/non-cadastral limitations and unstable raster-polygon identity declarations.
+- Kept `department-omgeving:land-use:forest:mol` separate from the historical land-use series and added a compact frontend series selector.
+
+Source validation:
+- Live WCS capabilities exposed all five expected coverages through `lu:lu_landgebruik_vlaa_<year>_v3`.
+- A full-Mol fetch-only run produced 3,768 / 3,823 / 3,953 / 3,615 / 3,676 forest polygons for 2013 / 2016 / 2019 / 2022 / 2025 without truncation.
+- Polygonized forest area was 3,723.17 / 3,615.39 / 3,592.27 / 3,648.98 / 3,626.56 ha. These are measurements within the official 10 m representation, not cadastral forest areas.
+- Every raster cell touching the municipality is considered before exact vector clipping. The 916 NoData edge cells in each WCS subset were explicitly excluded and recorded rather than assigned a class.
+
+Validation:
+- New focused backend suite passed 4 tests.
+- Frontend TypeScript typecheck passed after temporal-series selection wiring.
+- Full readiness and live Tower/PostGIS/browser validation remain the final steps of this pass.
+
+Next:
+- Deploy and provision all five modern snapshots on Tower, then verify current forest selection, both temporal series and a drawn rectangle in the internal browser.

docs/DATASET_STRATEGY.md

@@ -61,6 +61,33 @@ OSM is used for fast, broad, fallback context.
 ### OSM caveat
 OSM is community-maintained and may be incomplete. UI and reports must describe it as contextual/fallback data, not official ground truth.
 
+## Official modern land-use strategy
+
+The Departement Omgeving version 3 land-use maps are the authoritative modern
+source series for categorical land use. GeoIntel retrieves only explicit
+boundary subsets from MercatorNet WCS; no whole-Flanders raster is fetched at
+startup or during an interactive map query.
+
+- Source years: 2013, 2016, 2019, 2022 and 2025.
+- Native representation: one integer category per 10 m cell in EPSG:31370.
+- Preserved evidence: raw GeoTIFF, WCS request, catalogue URL, checksum, class
+  histogram and processing manifest.
+- Operational representation: source cells polygonized in EPSG:31370, clipped
+  to the requested boundary and normalized to EPSG:4326 `vector_features`.
+- Boundary behavior: retain every source cell touching the boundary, then
+  perform the exact metric geometry intersection; exclude and report NoData.
+- Selection metric: exact polygon intersection area within the source's 10 m
+  representation, reported in hectares with a non-cadastral warning.
+- Identity: raster polygons are not stable source objects; temporal object
+  lineage is unavailable.
+
+Modern land-use series must remain separate from the 1778, 1873 and 1969
+historical cartographic series. Comparisons are valid within each persisted
+series and must retain the source methodology warning. Extension from Mol to
+Kempen reuses the same operator with an explicit approved boundary, project,
+area and scope key; an ambiguous regional label is never converted into an
+invented boundary.
+
 ## User-uploaded raster strategy
 V1 must support controlled local datasets because public raster access and model compatibility can be difficult.
 

docs/DATA_SOURCES.md

@@ -63,6 +63,41 @@ The operator uses standards-compliant WFS 2.0 XML POST requests. This keeps
 the spatial/class filters server-side without exposing a long XML filter in a
 GET query, which the public gateway rejects.
 
+### Mol modern land use
+
+`scripts/provision_official_landuse_timeseries.py` uses the public Departement
+Omgeving/MercatorNet WCS to retrieve the harmonized version 3 land-use maps for
+2013, 2016, 2019, 2022 and 2025. The source is a categorical 10 m GeoTIFF in
+Belgian Lambert 72 (`EPSG:31370`) with 19 documented classes. GeoIntel currently
+derives only class `12` (`Bos`) as the operational forest theme.
+
+Every source raster is clipped against the explicit official boundary,
+validated for integer classes, CRS and resolution, checksummed and retained in
+persistent operator storage. Forest cells are polygonized and clipped in
+`EPSG:31370`, then transformed to canonical `EPSG:4326` and uploaded through
+DatasetService/VectorFeatureService. The raw raster remains the provenance
+artifact; the persisted polygons provide rectangle selection, PostGIS area
+aggregation and temporal comparison without a parallel query path.
+All cells touching the requested boundary are considered before the exact
+metric geometry clip; NoData cells are excluded and counted in the manifest.
+
+The modern series key is
+`department-omgeving:land-use:forest:mol`. It remains separate from
+`digitaal-vlaanderen:historical-landuse:forest:mol`: the 1778-1969
+cartographic editions and the harmonized 2013-2025 10 m land-use maps are not
+presented as one continuous measurement method. Raster-derived polygon
+identities are unstable, so GeoIntel compares hectares and never invents
+added/removed forest objects.
+
+Official catalogues:
+
+- https://www.vlaanderen.be/datavindplaats/catalogus/landgebruik-vlaanderen-toestand-2013
+- https://www.vlaanderen.be/datavindplaats/catalogus/landgebruik-vlaanderen-toestand-2016
+- https://www.vlaanderen.be/datavindplaats/catalogus/landgebruik-vlaanderen-toestand-2019
+- https://www.vlaanderen.be/datavindplaats/catalogus/landgebruik-vlaanderen-toestand-2022
+- https://www.vlaanderen.be/datavindplaats/catalogus/landgebruik-vlaanderen-toestand-2025
+- https://www.vlaanderen.be/statistiek-vlaanderen/ruimtegebruik/landgebruik/metadata-landgebruik
+
 ## OSM
 
 - Naam: OpenStreetMap

docs/TODO.md

@@ -6,8 +6,9 @@
 - [x] Add drag rectangle selection with automatic persisted theme queries.
 - [x] Separate exact bbox intersection totals from the bounded map preview.
 - [x] Add official Mol GRB roads, water and parcels provisioning support.
-- [ ] Select and validate an official Mol population/statistical-sector source before enabling the population theme.
+- [x] Select, validate and provision official Statbel 2021-2025 population/statistical-sector snapshots for Mol.
-- [ ] Select and validate an authoritative Flemish land-cover source before enabling the forest/green theme.
+- [x] Select, validate and provision official Departement Omgeving 2013-2025 forest snapshots from the 10 m land-use map.
+- [ ] Define the exact administrative Kempen scope before provisioning municipality or regional copies of the proven source series.
 
 This file now starts with the current implementation status. Older preparation/backlog sections are preserved below as historical planning context and should not be treated as the live sprint board without checking `docs/CODEX_EXECUTION_LOG.md`.
 

frontend/README.md

@@ -10,6 +10,10 @@ current, while `Evolution` lets the operator compare an earlier and later
 snapshot from the same series over a drawn rectangle. Results show units,
 absolute/percentage change, estimate status and source limitations.
 Added/removed/modified overlays only appear for stable source identities.
+When a theme has multiple valid methodologies, a compact series selector keeps
+them explicit. Forest therefore defaults to the official modern 2013-2025
+10 m series, while the separate 1778-1969 historical map series remains
+selectable and is never merged into the same trend.
 
 Selection results use dataset-specific PostGIS summaries. Object layers show
 intersecting counts, population shows inhabitants with partial-sector

frontend/src/components/map/MapWorkspace.tsx

@@ -7,6 +7,7 @@ import { useTemporalComparison } from '../../hooks/useTemporalComparison'
 
 const DEFAULT_SELECTED_FEATURE_FILENAME = 'selected-feature.geojson'
 const DEFAULT_AREA_SELECTION_FILENAME = 'area-selection.geojson'
+const EMPTY_TEMPORAL_SERIES: DatasetCreateResponse[] = []
 
 type DataThemeId = 'buildings' | 'population' | 'forest' | 'water' | 'roads' | 'parcels'
 
@@ -18,6 +19,12 @@ interface DataTheme {
   tokens: string[]
 }
 
+interface TemporalSeriesGroup {
+  key: string
+  label: string
+  items: DatasetCreateResponse[]
+}
+
 const DATA_THEMES: DataTheme[] = [
   {
     id: 'buildings',
@@ -90,6 +97,7 @@ function pickThemeDataset(datasets: DatasetCreateResponse[], theme: DataTheme):
     const score = (dataset: DatasetCreateResponse) =>
       (dataset.reference_layer_name && theme.tokens.includes(dataset.reference_layer_name.toLowerCase()) ? 1_000_000 : 0) +
       (dataset.source_name === 'grb' ? 100_000 : 0) +
+      (dataset.source_name === 'department_omgeving_land_use' ? 90_000 : 0) +
       (dataset.dataset_role === 'reference' ? 10_000 : 0) +
       (dataset.observed_at ? new Date(dataset.observed_at).getTime() / 100_000_000 : 0) +
       (dataset.feature_count ?? dataset.vector_summary?.feature_count ?? 0)
@@ -98,7 +106,21 @@ function pickThemeDataset(datasets: DatasetCreateResponse[], theme: DataTheme):
   return candidates[0] ?? null
 }
 
-function pickThemeTemporalSeries(datasets: DatasetCreateResponse[], theme: DataTheme): DatasetCreateResponse[] {
+function temporalSeriesLabel(items: DatasetCreateResponse[]): string {
+  const configuredLabel = items.find((item) => typeof item.source_metadata?.['temporal_series_label'] === 'string')
+    ?.source_metadata?.['temporal_series_label']
+  if (typeof configuredLabel === 'string' && configuredLabel.trim()) {
+    return configuredLabel
+  }
+  const first = items[0]
+  const source = first?.source_name ?? first?.source ?? 'Tijdreeks'
+  const firstYear = first?.observed_at ? new Date(first.observed_at).getUTCFullYear() : null
+  const last = items[items.length - 1]
+  const lastYear = last?.observed_at ? new Date(last.observed_at).getUTCFullYear() : null
+  return firstYear && lastYear ? `${source} (${firstYear}-${lastYear})` : source
+}
+
+function listThemeTemporalSeries(datasets: DatasetCreateResponse[], theme: DataTheme): TemporalSeriesGroup[] {
   const groups = new Map<string, DatasetCreateResponse[]>()
   for (const dataset of datasets) {
     if (!datasetMatchesTheme(dataset, theme) || !dataset.temporal_series_key || !dataset.observed_at) {
@@ -108,16 +130,21 @@ function pickThemeTemporalSeries(datasets: DatasetCreateResponse[], theme: DataT
     items.push(dataset)
     groups.set(dataset.temporal_series_key, items)
   }
-  return Array.from(groups.values())
+  return Array.from(groups.entries())
-    .filter((items) => items.length >= 2)
+    .filter(([, items]) => items.length >= 2)
+    .map(([key, items]) => {
+      const ordered = [...items].sort(
+        (left, right) => new Date(left.observed_at ?? 0).getTime() - new Date(right.observed_at ?? 0).getTime(),
+      )
+      return { key, label: temporalSeriesLabel(ordered), items: ordered }
+    })
     .sort((left, right) => {
-      if (right.length !== left.length) {
+      if (right.items.length !== left.items.length) {
-        return right.length - left.length
+        return right.items.length - left.items.length
       }
-      const latest = (items: DatasetCreateResponse[]) => Math.max(...items.map((item) => new Date(item.observed_at ?? 0).getTime()))
+      const latest = (group: TemporalSeriesGroup) => Math.max(...group.items.map((item) => new Date(item.observed_at ?? 0).getTime()))
       return latest(right) - latest(left)
-    })[0]
+    })
-    ?.sort((left, right) => new Date(left.observed_at ?? 0).getTime() - new Date(right.observed_at ?? 0).getTime()) ?? []
 }
 
 function formatObservationDate(value: string | null | undefined): string {
@@ -496,6 +523,7 @@ export function MapWorkspace({
     clearTemporalComparison,
   } = useTemporalComparison(selectedProjectId)
   const [analysisMode, setAnalysisMode] = useState<'current' | 'evolution'>('current')
+  const [selectedTemporalSeriesKey, setSelectedTemporalSeriesKey] = useState('')
   const [earlierDatasetId, setEarlierDatasetId] = useState('')
   const [laterDatasetId, setLaterDatasetId] = useState('')
   const [bboxSelectionMode, setBboxSelectionMode] = useState(false)
@@ -536,10 +564,13 @@ export function MapWorkspace({
   )
   const activeTheme = DATA_THEMES.find((theme) => theme.id === activeThemeId) ?? DATA_THEMES[0]
   const activeThemeDataset = themeDatasetMap[activeTheme.id]
-  const activeTemporalSeries = useMemo(
+  const activeTemporalSeriesGroups = useMemo(
-    () => pickThemeTemporalSeries(availableMapDatasets, activeTheme),
+    () => listThemeTemporalSeries(availableMapDatasets, activeTheme),
     [activeTheme, availableMapDatasets],
   )
+  const activeTemporalSeriesGroup = activeTemporalSeriesGroups.find((group) => group.key === selectedTemporalSeriesKey)
+    ?? activeTemporalSeriesGroups[0]
+  const activeTemporalSeries = activeTemporalSeriesGroup?.items ?? EMPTY_TEMPORAL_SERIES
   const themeResults = useMemo(
     () =>
       themeInsights.flatMap((insight) => {
@@ -592,6 +623,14 @@ export function MapWorkspace({
     setBboxInput(bboxToInputState(mapSelectionBbox))
   }, [mapSelectionBbox])
 
+  useEffect(() => {
+    setSelectedTemporalSeriesKey((current) =>
+      activeTemporalSeriesGroups.some((group) => group.key === current)
+        ? current
+        : activeTemporalSeriesGroups[0]?.key ?? '',
+    )
+  }, [activeTemporalSeriesGroups])
+
   useEffect(() => {
     const first = activeTemporalSeries[0]
     const last = activeTemporalSeries[activeTemporalSeries.length - 1]
@@ -912,7 +951,7 @@ export function MapWorkspace({
               <span>{analysisMode === 'evolution' ? 'Tijdreeks' : 'Actieve bron'}</span>
               <strong>
                 {analysisMode === 'evolution'
-                  ? activeTemporalSeries[0]?.temporal_series_key ?? 'Geen tijdreeks beschikbaar'
+                  ? activeTemporalSeriesGroup?.label ?? 'Geen tijdreeks beschikbaar'
                   : activeThemeDataset?.name ?? 'Geen databron beschikbaar'}
               </strong>
               <small>
@@ -928,6 +967,22 @@ export function MapWorkspace({
 
             {analysisMode === 'evolution' ? (
               <div className="geo-time-controls" aria-label="Meetmomenten vergelijken">
+                {activeTemporalSeriesGroups.length > 1 ? (
+                  <label className="geo-series-control">
+                    Reeks
+                    <select
+                      value={activeTemporalSeriesGroup?.key ?? ''}
+                      onChange={(event) => {
+                        setSelectedTemporalSeriesKey(event.target.value)
+                        clearTemporalComparison()
+                      }}
+                    >
+                      {activeTemporalSeriesGroups.map((group) => (
+                        <option key={group.key} value={group.key}>{group.label}</option>
+                      ))}
+                    </select>
+                  </label>
+                ) : null}
                 <label>
                   Van
                   <select value={earlierDatasetId} onChange={(event) => { setEarlierDatasetId(event.target.value); clearTemporalComparison() }} disabled={activeTemporalSeries.length < 2}>
@@ -1208,7 +1263,7 @@ export function MapWorkspace({
           <span>
             <strong>Bron:</strong>{' '}
             {analysisMode === 'evolution'
-              ? activeTemporalSeries[0]?.temporal_series_key ?? 'geen vergelijkbare tijdreeks'
+              ? activeTemporalSeriesGroup?.label ?? 'geen vergelijkbare tijdreeks'
               : activeThemeDataset
                 ? `${activeThemeDataset.source_name ?? activeThemeDataset.source} · ${activeThemeDataset.name}`
                 : 'niet beschikbaar'}

frontend/src/styles/app.css

@@ -5675,6 +5675,10 @@ section {
   font-size: 0.68rem;
 }
 
+.geo-time-controls .geo-series-control {
+  grid-column: 1 / -1;
+}
+
 .geo-time-controls button {
   grid-column: 1 / -1;
   min-height: 2.25rem;

scripts/README.md

@@ -1212,9 +1212,47 @@ normal dataset API. Artifacts and manifests are retained below
 datasets; use `--force` only for an explicit source refresh. Use
 `--layers roads,water` or `--fetch-only` for a bounded operator run.
 
-The context provisioner does not add population or forest values. Those themes
+The context provisioner itself does not add population or forest values. Use
-remain visibly unavailable until an authoritative/statistically appropriate
+the dedicated official-source operators below; zero is never substituted for
-source is imported; zero is never substituted for missing source data.
+missing source data.
+
+## Official Mol temporal sources
+
+After the municipality workspace exists, provision official Statbel population
+and both independently modelled land-use series:
+
+```bash
+docker exec -it geointel python3 /app/scripts/provision_mol_population_history.py
+docker exec -it geointel python3 /app/scripts/provision_mol_historical_landuse.py
+docker exec -it geointel python3 /app/scripts/provision_official_landuse_timeseries.py
+```
+
+The modern land-use command checks the MercatorNet WCS capabilities, downloads
+only the Mol bounding subset of each 10 m `EPSG:31370` raster, validates the
+categorical integer grid, clips against the official boundary and polygonizes
+class `12` (`Bos`). Raw rasters, vector artifacts and checksum manifests are
+stored under `/app/storage/operator-data/official-landuse/mol`. The resulting
+2013, 2016, 2019, 2022 and 2025 vectors are uploaded through the canonical API
+as `department-omgeving:land-use:forest:mol`.
+
+Prepare and inspect artifacts without changing the database:
+
+```bash
+docker exec -it geointel python3 \
+  /app/scripts/provision_official_landuse_timeseries.py --fetch-only
+```
+
+Use `--force` only to refetch and rebuild local artifacts. Existing persisted
+snapshots remain immutable and are reused by year/series. To use the operator
+for another approved region, pass all scope inputs explicitly, for example
+`--boundary-path`, `--project-name`, `--area-name`, `--municipality-name`,
+`--nis-code`, `--scope-key` and `--output-dir`. GeoIntel does not infer what
+"Kempen" means administratively.
+
+The 2013-2025 series is methodologically separate from the historical
+1778/1873/1969 series. The map-first Evolution view exposes a series selector
+when both exist; it never calculates one continuous trend across those source
+families.
 
 ## Tower deployment
 

scripts/provision_official_landuse_timeseries.py

@@ -0,0 +1,817 @@
+"""Provision official modern Flemish land-use snapshots for an explicit area.
+
+The operator downloads categorical 10 metre GeoTIFF subsets from the public
+Departement Omgeving MercatorNet WCS, clips them to a supplied boundary and
+polygonizes only explicitly supported classes. Raw rasters and checksum
+manifests remain provenance artifacts. Vector output is imported through the
+normal GeoIntel dataset API and is never written directly to PostGIS.
+
+Mol is the safe default. Other scopes must provide their own boundary, project,
+area name and identity explicitly.
+"""
+
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+import math
+import os
+import sys
+from dataclasses import dataclass
+from datetime import datetime, timezone
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+import rasterio
+import requests
+from pyproj import Transformer
+from rasterio.features import geometry_mask, shapes
+from requests.adapters import HTTPAdapter
+from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, mapping, shape
+from shapely.ops import transform, unary_union
+from shapely.validation import make_valid
+from urllib3.util.retry import Retry
+
+
+DEFAULT_API_URL = "http://127.0.0.1:8000"
+DEFAULT_PROJECT_NAME = "Mol Municipality Workbench"
+DEFAULT_AREA_NAME = "Gemeente Mol"
+DEFAULT_MUNICIPALITY_NAME = "Mol"
+DEFAULT_NIS_CODE = "13025"
+DEFAULT_SCOPE_KEY = "mol"
+DEFAULT_BOUNDARY_PATH = Path("/app/storage/operator-data/mol-municipality/mol_municipality_boundary.geojson")
+DEFAULT_OUTPUT_DIR = Path("/app/storage/operator-data/official-landuse/mol")
+
+WCS_URL = "https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs"
+WCS_VERSION = "1.0.0"
+SOURCE_CRS = "EPSG:31370"
+OUTPUT_CRS = "EPSG:4326"
+SOURCE_RESOLUTION_METRES = 10.0
+SUPPORTED_YEARS = (2013, 2016, 2019, 2022, 2025)
+ATTRIBUTION = "Bron: Landgebruik Vlaanderen, Departement Omgeving"
+SERIES_LABEL = "Moderne landgebruikskaart (10 m)"
+
+CATALOGUE_URLS = {
+    year: f"https://www.vlaanderen.be/datavindplaats/catalogus/landgebruik-vlaanderen-toestand-{year}"
+    for year in SUPPORTED_YEARS
+}
+
+LAND_USE_CLASSES = {
+    1: "Huizen en tuinen",
+    2: "Industrie",
+    3: "Commerciele doeleinden",
+    4: "Diensten",
+    5: "Transportinfrastructuur",
+    6: "Recreatie",
+    7: "Landbouwgebouwen en -infrastructuur",
+    8: "Overige bebouwde terreinen",
+    9: "Overige onbebouwde terreinen",
+    10: "Actieve groeves",
+    11: "Luchthavens",
+    12: "Bos",
+    13: "Akker",
+    14: "Grasland in landbouwgebruik",
+    15: "Struikgewas",
+    16: "Braakliggend en duinen",
+    17: "Water",
+    18: "Moeras",
+    19: "Overige graslanden",
+}
+
+
+@dataclass(frozen=True)
+class ThemeDefinition:
+    key: str
+    label: str
+    class_ids: tuple[int, ...]
+
+
+@dataclass(frozen=True)
+class PreparedSnapshot:
+    year: int
+    theme: ThemeDefinition
+    raster_path: Path
+    vector_path: Path
+    manifest_path: Path
+    feature_count: int
+    raster_sha256: str
+    vector_sha256: str
+
+
+THEMES = (ThemeDefinition("forest", "Bos", (12,)),)
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description="Provision official 2013-2025 Flemish land-use snapshots.")
+    parser.add_argument("--base-url", default=os.environ.get("GEOINTEL_INTERNAL_API_URL", DEFAULT_API_URL))
+    parser.add_argument("--project-name", default=DEFAULT_PROJECT_NAME)
+    parser.add_argument("--area-name", default=DEFAULT_AREA_NAME, help="Case-insensitive fragment identifying the persisted Area.")
+    parser.add_argument("--municipality-name", default=DEFAULT_MUNICIPALITY_NAME)
+    parser.add_argument("--nis-code", default=DEFAULT_NIS_CODE)
+    parser.add_argument("--scope-key", default=DEFAULT_SCOPE_KEY)
+    parser.add_argument("--years", default=",".join(str(year) for year in SUPPORTED_YEARS))
+    parser.add_argument("--themes", default="forest")
+    parser.add_argument(
+        "--boundary-path",
+        type=Path,
+        default=Path(os.environ.get("OFFICIAL_LANDUSE_BOUNDARY_PATH", DEFAULT_BOUNDARY_PATH)),
+    )
+    parser.add_argument(
+        "--output-dir",
+        type=Path,
+        default=Path(os.environ.get("OFFICIAL_LANDUSE_OUTPUT_DIR", DEFAULT_OUTPUT_DIR)),
+    )
+    parser.add_argument("--request-timeout", type=int, default=240)
+    parser.add_argument("--import-timeout", type=int, default=1800)
+    parser.add_argument("--max-features", type=int, default=100000)
+    parser.add_argument("--force", action="store_true", help="Refetch and rebuild local source artifacts; persisted datasets stay immutable.")
+    parser.add_argument("--fetch-only", action="store_true", help="Prepare and verify artifacts without changing GeoIntel persistence.")
+    return parser.parse_args()
+
+
+def utc_now() -> str:
+    return datetime.now(timezone.utc).isoformat()
+
+
+def sha256_file(path: Path) -> str:
+    digest = hashlib.sha256()
+    with path.open("rb") as handle:
+        for chunk in iter(lambda: handle.read(1024 * 1024), b""):
+            digest.update(chunk)
+    return digest.hexdigest()
+
+
+def write_json_atomic(path: Path, payload: dict[str, Any], *, pretty: bool = False) -> None:
+    temporary = path.with_suffix(f"{path.suffix}.partial")
+    temporary.write_text(
+        json.dumps(
+            payload,
+            ensure_ascii=False,
+            indent=2 if pretty else None,
+            separators=None if pretty else (",", ":"),
+            sort_keys=pretty,
+        ),
+        encoding="utf-8",
+    )
+    temporary.replace(path)
+
+
+def build_session() -> requests.Session:
+    retry = Retry(
+        total=5,
+        connect=5,
+        read=5,
+        status=5,
+        backoff_factor=1.0,
+        status_forcelist=(429, 500, 502, 503, 504),
+        allowed_methods=frozenset({"GET"}),
+        raise_on_status=True,
+    )
+    session = requests.Session()
+    session.headers.update({"User-Agent": "GeoIntel-Official-Landuse-Operator/1.0"})
+    adapter = HTTPAdapter(max_retries=retry)
+    session.mount("https://", adapter)
+    session.mount("http://", adapter)
+    return session
+
+
+def coverage_id(year: int) -> str:
+    return f"lu:lu_landgebruik_vlaa_{year}_v3"
+
+
+def series_key(theme: ThemeDefinition, scope_key: str) -> str:
+    return f"department-omgeving:land-use:{theme.key}:{scope_key.strip().lower()}"
+
+
+def load_boundary(path: Path):
+    if not path.exists():
+        raise RuntimeError(f"Boundary artifact is missing at {path}")
+    payload = json.loads(path.read_text(encoding="utf-8"))
+    features = payload.get("features") or []
+    if len(features) != 1:
+        raise RuntimeError("Boundary artifact must contain exactly one feature")
+    boundary = normalize_polygonal(shape(features[0].get("geometry")))
+    if boundary is None:
+        raise RuntimeError("Boundary artifact is empty, invalid or non-polygonal")
+    return boundary
+
+
+def normalize_polygonal(geometry):
+    if geometry is None or geometry.is_empty:
+        return None
+    if not geometry.is_valid:
+        geometry = make_valid(geometry)
+    if isinstance(geometry, (Polygon, MultiPolygon)):
+        return geometry
+    if isinstance(geometry, GeometryCollection):
+        polygons = [part for part in geometry.geoms if isinstance(part, (Polygon, MultiPolygon)) and not part.is_empty]
+        if not polygons:
+            return None
+        merged = unary_union(polygons)
+        return merged if isinstance(merged, (Polygon, MultiPolygon)) and not merged.is_empty else None
+    return None
+
+
+def metric_boundary(boundary):
+    transformer = Transformer.from_crs(OUTPUT_CRS, SOURCE_CRS, always_xy=True)
+    projected = normalize_polygonal(transform(transformer.transform, boundary))
+    if projected is None:
+        raise RuntimeError("Boundary could not be projected to EPSG:31370")
+    return projected
+
+
+def snapped_bounds(bounds: tuple[float, float, float, float]) -> tuple[float, float, float, float]:
+    min_x, min_y, max_x, max_y = bounds
+    resolution = SOURCE_RESOLUTION_METRES
+    return (
+        math.floor(min_x / resolution) * resolution,
+        math.floor(min_y / resolution) * resolution,
+        math.ceil(max_x / resolution) * resolution,
+        math.ceil(max_y / resolution) * resolution,
+    )
+
+
+def build_wcs_params(year: int, bounds: tuple[float, float, float, float]) -> dict[str, str]:
+    min_x, min_y, max_x, max_y = snapped_bounds(bounds)
+    return {
+        "SERVICE": "WCS",
+        "VERSION": WCS_VERSION,
+        "REQUEST": "GetCoverage",
+        "COVERAGE": coverage_id(year),
+        "CRS": SOURCE_CRS,
+        "BBOX": f"{min_x:.3f},{min_y:.3f},{max_x:.3f},{max_y:.3f}",
+        "RESX": str(int(SOURCE_RESOLUTION_METRES)),
+        "RESY": str(int(SOURCE_RESOLUTION_METRES)),
+        "FORMAT": "image/tiff",
+        "RESPONSE_CRS": SOURCE_CRS,
+    }
+
+
+def verify_coverages(session: requests.Session, years: list[int], timeout: int) -> None:
+    response = session.get(
+        WCS_URL,
+        params={"SERVICE": "WCS", "VERSION": WCS_VERSION, "REQUEST": "GetCapabilities"},
+        timeout=timeout,
+    )
+    response.raise_for_status()
+    missing = [coverage_id(year) for year in years if coverage_id(year) not in response.text]
+    if missing:
+        raise RuntimeError(f"Official WCS is missing expected coverages: {', '.join(missing)}")
+
+
+def validate_raster(path: Path) -> dict[str, Any]:
+    try:
+        with rasterio.open(path) as dataset:
+            epsg = dataset.crs.to_epsg() if dataset.crs else None
+            resolution = (abs(float(dataset.res[0])), abs(float(dataset.res[1])))
+            if epsg != 31370:
+                raise RuntimeError(f"Expected EPSG:31370 source raster, received {dataset.crs}")
+            if dataset.count != 1:
+                raise RuntimeError(f"Expected one categorical raster band, received {dataset.count}")
+            if any(abs(value - SOURCE_RESOLUTION_METRES) > 0.01 for value in resolution):
+                raise RuntimeError(f"Expected 10 metre source resolution, received {resolution}")
+            if not np.issubdtype(np.dtype(dataset.dtypes[0]), np.integer):
+                raise RuntimeError(f"Expected integer land-use classes, received {dataset.dtypes[0]}")
+            return {
+                "width": dataset.width,
+                "height": dataset.height,
+                "dtype": dataset.dtypes[0],
+                "nodata": dataset.nodata,
+                "crs": SOURCE_CRS,
+                "resolution_metres": SOURCE_RESOLUTION_METRES,
+                "bounds": list(dataset.bounds),
+            }
+    except rasterio.errors.RasterioError as exc:
+        raise RuntimeError(f"Official WCS response is not a readable GeoTIFF: {exc}") from exc
+
+
+def download_raster(
+    session: requests.Session,
+    *,
+    year: int,
+    bounds: tuple[float, float, float, float],
+    path: Path,
+    timeout: int,
+) -> dict[str, Any]:
+    response = session.get(WCS_URL, params=build_wcs_params(year, bounds), timeout=timeout, stream=True)
+    response.raise_for_status()
+    content_type = str(response.headers.get("content-type") or "").lower()
+    if "tiff" not in content_type:
+        preview = response.content[:500].decode("utf-8", errors="replace")
+        raise RuntimeError(f"Official WCS returned {content_type or 'unknown content'} instead of GeoTIFF: {preview}")
+    temporary = path.with_suffix(f"{path.suffix}.partial")
+    try:
+        with temporary.open("wb") as handle:
+            for chunk in response.iter_content(chunk_size=1024 * 1024):
+                if chunk:
+                    handle.write(chunk)
+        profile = validate_raster(temporary)
+        temporary.replace(path)
+    except Exception:
+        temporary.unlink(missing_ok=True)
+        raise
+    return {**profile, "request_url": response.url, "retrieved_at": utc_now()}
+
+
+def polygonize_snapshot(
+    *,
+    raster_path: Path,
+    boundary,
+    year: int,
+    theme: ThemeDefinition,
+    municipality_name: str,
+    nis_code: str,
+    scope_key: str,
+    max_features: int,
+) -> tuple[dict[str, Any], dict[str, Any]]:
+    boundary_metric = metric_boundary(boundary)
+    to_output = Transformer.from_crs(SOURCE_CRS, OUTPUT_CRS, always_xy=True)
+    raster_sha256 = sha256_file(raster_path)
+
+    with rasterio.open(raster_path) as dataset:
+        validate_raster(raster_path)
+        values = dataset.read(1)
+        inside = geometry_mask(
+            [mapping(boundary_metric)],
+            out_shape=values.shape,
+            transform=dataset.transform,
+            invert=True,
+            # Exact clipping happens after polygonization, so retain every
+            # classified source cell that overlaps the requested boundary.
+            all_touched=True,
+        )
+        valid_inside = inside.copy()
+        nodata_count = 0
+        if dataset.nodata is not None:
+            nodata = np.isclose(values, dataset.nodata)
+            nodata_count = int(np.count_nonzero(inside & nodata))
+            valid_inside &= ~nodata
+        available_values, available_counts = np.unique(values[valid_inside], return_counts=True)
+        class_histogram = {str(int(value)): int(count) for value, count in zip(available_values, available_counts)}
+        unknown_classes = sorted(int(value) for value in available_values if int(value) not in LAND_USE_CLASSES)
+        if unknown_classes:
+            raise RuntimeError(f"Land-use raster {year} contains undocumented classes: {unknown_classes}")
+
+        class_mask = np.isin(values, np.asarray(theme.class_ids)) & valid_inside
+        source_pixel_count = int(np.count_nonzero(class_mask))
+        if source_pixel_count == 0:
+            raise RuntimeError(f"Land-use raster {year} contains no {theme.label} cells inside the boundary")
+
+        features: list[dict[str, Any]] = []
+        polygon_area_m2 = 0.0
+        for raw_geometry, raw_value in shapes(
+            class_mask.astype("uint8"),
+            mask=class_mask,
+            transform=dataset.transform,
+            connectivity=4,
+        ):
+            if int(raw_value) != 1:
+                continue
+            geometry_metric = normalize_polygonal(shape(raw_geometry).intersection(boundary_metric))
+            if geometry_metric is None or geometry_metric.area <= 0:
+                continue
+            geometry_output = normalize_polygonal(transform(to_output.transform, geometry_metric))
+            if geometry_output is None:
+                continue
+            if len(features) >= max_features:
+                raise RuntimeError(
+                    f"Land-use {theme.key} {year} exceeds the {max_features} feature safety limit; refusing truncation"
+                )
+            feature_hash = hashlib.sha256(
+                f"{year}:{theme.key}:".encode("utf-8") + geometry_metric.wkb
+            ).hexdigest()[:24]
+            feature_id = f"land-use-{year}-{theme.key}-{feature_hash}"
+            area_m2 = float(geometry_metric.area)
+            polygon_area_m2 += area_m2
+            features.append(
+                {
+                    "type": "Feature",
+                    "id": feature_id,
+                    "geometry": mapping(geometry_output),
+                    "properties": {
+                        "source_name": "department_omgeving_land_use",
+                        "source_feature_id": feature_id,
+                        "reference_layer_name": theme.key,
+                        "layer_type": theme.key,
+                        "authority_level": "authoritative",
+                        "coverage_scope": scope_key,
+                        "municipality": municipality_name,
+                        "nis_code": nis_code,
+                        "observation_year": year,
+                        "land_use_class_ids": list(theme.class_ids),
+                        "land_use_class_names": [LAND_USE_CLASSES[class_id] for class_id in theme.class_ids],
+                        "source_resolution_m": SOURCE_RESOLUTION_METRES,
+                        "polygon_area_m2": round(area_m2, 3),
+                        "source_coverage_id": coverage_id(year),
+                        "source_raster_sha256": raster_sha256,
+                        "attribution": ATTRIBUTION,
+                    },
+                }
+            )
+
+    payload = {
+        "type": "FeatureCollection",
+        "name": f"{theme.label} - {municipality_name} {year}",
+        "crs": {"type": "name", "properties": {"name": OUTPUT_CRS}},
+        "municipality": municipality_name,
+        "nis_code": nis_code,
+        "scope_key": scope_key,
+        "observation_year": year,
+        "source_coverage_id": coverage_id(year),
+        "source_resolution_m": SOURCE_RESOLUTION_METRES,
+        "attribution": ATTRIBUTION,
+        "features": features,
+    }
+    stats = {
+        "feature_count": len(features),
+        "source_pixel_count": source_pixel_count,
+        "source_pixel_area_m2": source_pixel_count * SOURCE_RESOLUTION_METRES**2,
+        "polygon_area_m2": polygon_area_m2,
+        "nodata_pixels_inside_boundary": nodata_count,
+        "class_histogram": class_histogram,
+        "raster_sha256": raster_sha256,
+    }
+    return payload, stats
+
+
+def existing_snapshot(
+    *,
+    year: int,
+    theme: ThemeDefinition,
+    raster_path: Path,
+    vector_path: Path,
+    manifest_path: Path,
+) -> PreparedSnapshot | None:
+    if not raster_path.exists() or not vector_path.exists() or not manifest_path.exists():
+        return None
+    try:
+        manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
+        raster_sha256 = sha256_file(raster_path)
+        vector_sha256 = sha256_file(vector_path)
+        if (
+            manifest.get("year") != year
+            or manifest.get("theme") != theme.key
+            or manifest.get("coverage_id") != coverage_id(year)
+            or manifest.get("class_ids") != list(theme.class_ids)
+            or manifest.get("raster_sha256") != raster_sha256
+            or manifest.get("vector_sha256") != vector_sha256
+        ):
+            return None
+        validate_raster(raster_path)
+        feature_count = int(manifest["feature_count"])
+        if feature_count <= 0:
+            return None
+        return PreparedSnapshot(
+            year=year,
+            theme=theme,
+            raster_path=raster_path,
+            vector_path=vector_path,
+            manifest_path=manifest_path,
+            feature_count=feature_count,
+            raster_sha256=raster_sha256,
+            vector_sha256=vector_sha256,
+        )
+    except (OSError, ValueError, KeyError, RuntimeError, rasterio.errors.RasterioError):
+        return None
+
+
+def prepare_snapshot(
+    session: requests.Session,
+    *,
+    args: argparse.Namespace,
+    boundary,
+    boundary_metric,
+    year: int,
+    theme: ThemeDefinition,
+) -> PreparedSnapshot:
+    stem = f"{args.scope_key}_land_use_{theme.key}_{year}"
+    raster_path = args.output_dir / f"{stem}.tif"
+    vector_path = args.output_dir / f"{stem}.geojson"
+    manifest_path = args.output_dir / f"{stem}.manifest.json"
+    if not args.force:
+        prepared = existing_snapshot(
+            year=year,
+            theme=theme,
+            raster_path=raster_path,
+            vector_path=vector_path,
+            manifest_path=manifest_path,
+        )
+        if prepared:
+            return prepared
+
+    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,
+        )
+    else:
+        raster_profile = validate_raster(raster_path)
+        raster_profile["request_url"] = requests.Request(
+            "GET", WCS_URL, params=build_wcs_params(year, boundary_metric.bounds)
+        ).prepare().url
+        raster_profile["retrieved_at"] = None
+
+    payload, stats = polygonize_snapshot(
+        raster_path=raster_path,
+        boundary=boundary,
+        year=year,
+        theme=theme,
+        municipality_name=args.municipality_name,
+        nis_code=args.nis_code,
+        scope_key=args.scope_key,
+        max_features=args.max_features,
+    )
+    write_json_atomic(vector_path, payload)
+    vector_sha256 = sha256_file(vector_path)
+    manifest = {
+        "schema_version": 1,
+        "year": year,
+        "theme": theme.key,
+        "class_ids": list(theme.class_ids),
+        "class_names": [LAND_USE_CLASSES[class_id] for class_id in theme.class_ids],
+        "coverage_id": coverage_id(year),
+        "catalogue_url": CATALOGUE_URLS[year],
+        "wcs_url": WCS_URL,
+        "wcs_version": WCS_VERSION,
+        "wcs_request_url": raster_profile.get("request_url"),
+        "source_crs": SOURCE_CRS,
+        "output_crs": OUTPUT_CRS,
+        "source_resolution_metres": SOURCE_RESOLUTION_METRES,
+        "boundary_path": str(args.boundary_path),
+        "municipality": args.municipality_name,
+        "nis_code": args.nis_code,
+        "scope_key": args.scope_key,
+        "raster_path": str(raster_path),
+        "vector_path": str(vector_path),
+        "raster_profile": raster_profile,
+        "raster_sha256": stats["raster_sha256"],
+        "vector_sha256": vector_sha256,
+        "feature_count": stats["feature_count"],
+        "source_pixel_count": stats["source_pixel_count"],
+        "source_pixel_area_m2": stats["source_pixel_area_m2"],
+        "polygon_area_m2": stats["polygon_area_m2"],
+        "nodata_pixels_inside_boundary": stats["nodata_pixels_inside_boundary"],
+        "class_histogram": stats["class_histogram"],
+        "generated_at": utc_now(),
+    }
+    write_json_atomic(manifest_path, manifest, pretty=True)
+    return PreparedSnapshot(
+        year=year,
+        theme=theme,
+        raster_path=raster_path,
+        vector_path=vector_path,
+        manifest_path=manifest_path,
+        feature_count=int(stats["feature_count"]),
+        raster_sha256=str(stats["raster_sha256"]),
+        vector_sha256=vector_sha256,
+    )
+
+
+def response_data(response: requests.Response) -> Any:
+    try:
+        payload = response.json()
+    except ValueError as exc:
+        raise RuntimeError(f"GeoIntel API returned non-JSON ({response.status_code}): {response.text[:300]}") from exc
+    if not response.ok:
+        raise RuntimeError(f"GeoIntel API failed ({response.status_code}): {json.dumps(payload, ensure_ascii=False)[:800]}")
+    if not isinstance(payload, dict) or "data" not in payload:
+        raise RuntimeError("GeoIntel API response does not use the canonical data envelope")
+    return payload["data"]
+
+
+def locate_workspace(session: requests.Session, base_url: str, args: argparse.Namespace):
+    projects = response_data(session.get(f"{base_url}/api/v1/projects", params={"limit": 200}, timeout=args.import_timeout))
+    project = next((item for item in projects.get("items") or [] if item.get("name") == args.project_name), None)
+    if not project:
+        raise RuntimeError(f"Project {args.project_name!r} is missing")
+    project_id = str(project["id"])
+    areas = response_data(
+        session.get(f"{base_url}/api/v1/projects/{project_id}/areas", params={"limit": 200}, timeout=args.import_timeout)
+    )
+    area_fragment = args.area_name.strip().casefold()
+    matches = [item for item in areas.get("items") or [] if area_fragment in str(item.get("name") or "").casefold()]
+    if len(matches) != 1:
+        raise RuntimeError(f"Expected one Area matching {args.area_name!r}, received {len(matches)}")
+    datasets = response_data(
+        session.get(f"{base_url}/api/v1/projects/{project_id}/datasets", params={"limit": 500}, timeout=args.import_timeout)
+    )
+    return project_id, str(matches[0]["id"]), list(datasets.get("items") or [])
+
+
+def build_source_metadata(args: argparse.Namespace, snapshot: PreparedSnapshot) -> dict[str, Any]:
+    return {
+        "provider": "Departement Omgeving",
+        "source_title": f"Landgebruik - Vlaanderen - toestand {snapshot.year}",
+        "catalogue_url": CATALOGUE_URLS[snapshot.year],
+        "coverage_id": coverage_id(snapshot.year),
+        "authority_level": "authoritative",
+        "coverage_scope": args.scope_key,
+        "municipality": args.municipality_name,
+        "nis_code": args.nis_code,
+        "attribution": ATTRIBUTION,
+        "license_note": "Publieke Vlaamse overheidsdata; raadpleeg de toegangs- en gebruiksvoorwaarden in de bronmetadata.",
+        "methodology_version": "3",
+        "source_resolution_metres": SOURCE_RESOLUTION_METRES,
+        "source_crs": SOURCE_CRS,
+        "polygon_crs": OUTPUT_CRS,
+        "land_use_class_ids": list(snapshot.theme.class_ids),
+        "land_use_class_names": [LAND_USE_CLASSES[class_id] for class_id in snapshot.theme.class_ids],
+        "temporal_series_label": SERIES_LABEL,
+        "observation_date_precision": "year",
+        "identity_stable": False,
+        "identity_limitation": "Raster-derived polygons can split or merge between source editions; object lineage is not inferred.",
+        "selection_aggregation": {
+            "method": "intersection_area",
+            "label": "Oppervlakte",
+            "unit": "ha",
+            "is_estimate": False,
+            "warning": "Oppervlakte is exact binnen de officiele 10 m rasterrepresentatie en is niet perceelsnauwkeurig.",
+        },
+        "comparison_limitation": "Compare editions as 10 m land-use states; source inputs and methodology can evolve between publication years.",
+    }
+
+
+def build_provenance_metadata(args: argparse.Namespace, snapshot: PreparedSnapshot) -> dict[str, Any]:
+    return {
+        "operator_tool": "provision_official_landuse_timeseries.py",
+        "operator_explicit_fetch": True,
+        "wcs_url": WCS_URL,
+        "wcs_version": WCS_VERSION,
+        "coverage_id": coverage_id(snapshot.year),
+        "catalogue_url": CATALOGUE_URLS[snapshot.year],
+        "raw_raster_path": str(snapshot.raster_path),
+        "polygon_artifact_path": str(snapshot.vector_path),
+        "manifest_path": str(snapshot.manifest_path),
+        "raster_sha256": snapshot.raster_sha256,
+        "vector_sha256": snapshot.vector_sha256,
+        "source_crs": SOURCE_CRS,
+        "output_crs": OUTPUT_CRS,
+        "source_resolution_metres": SOURCE_RESOLUTION_METRES,
+        "generated_at": utc_now(),
+    }
+
+
+def upload_snapshot(
+    session: requests.Session,
+    *,
+    base_url: str,
+    project_id: str,
+    area_id: str,
+    args: argparse.Namespace,
+    snapshot: PreparedSnapshot,
+) -> dict[str, Any]:
+    observed_at = f"{snapshot.year}-01-01T00:00:00Z"
+    with snapshot.vector_path.open("rb") as handle:
+        response = session.post(
+            f"{base_url}/api/v1/projects/{project_id}/datasets/upload",
+            data={
+                "dataset_type": "vector",
+                "source": "operator_official_import",
+                "dataset_role": "reference",
+                "source_name": "department_omgeving_land_use",
+                "reference_layer_name": snapshot.theme.key,
+                "source_metadata_json": json.dumps(build_source_metadata(args, snapshot), ensure_ascii=False),
+                "provenance_metadata_json": json.dumps(build_provenance_metadata(args, snapshot), ensure_ascii=False),
+                "area_id": area_id,
+                "temporal_series_key": series_key(snapshot.theme, args.scope_key),
+                "observed_at": observed_at,
+                "valid_from": observed_at,
+                "temporal_granularity": "year",
+                "source_version": f"{snapshot.year}-v3",
+            },
+            files={"file": (snapshot.vector_path.name, handle, "application/geo+json")},
+            timeout=args.import_timeout,
+        )
+    return response_data(response)
+
+
+def main() -> int:
+    args = parse_args()
+    try:
+        years = sorted({int(value.strip()) for value in args.years.split(",") if value.strip()})
+    except ValueError:
+        print(json.dumps({"status": "error", "message": "Years must be comma-separated integers"}), file=sys.stderr)
+        return 2
+    requested_themes = {value.strip().lower() for value in args.themes.split(",") if value.strip()}
+    definitions = [definition for definition in THEMES if definition.key in requested_themes]
+    unsupported_years = [year for year in years if year not in SUPPORTED_YEARS]
+    unsupported_themes = requested_themes - {definition.key for definition in THEMES}
+    if unsupported_years or unsupported_themes or not years or not definitions:
+        print(
+            json.dumps(
+                {"status": "error", "message": f"Unsupported years={unsupported_years}, themes={sorted(unsupported_themes)}"}
+            ),
+            file=sys.stderr,
+        )
+        return 2
+    if not args.scope_key.strip() or not args.project_name.strip() or not args.area_name.strip():
+        print(json.dumps({"status": "error", "message": "scope-key, project-name and area-name are required"}), file=sys.stderr)
+        return 2
+
+    args.output_dir.mkdir(parents=True, exist_ok=True)
+    results: list[dict[str, Any]] = []
+    try:
+        boundary = load_boundary(args.boundary_path)
+        boundary_metric = metric_boundary(boundary)
+        prepared: list[PreparedSnapshot] = []
+        with build_session() as source_session:
+            verify_coverages(source_session, years, args.request_timeout)
+            for year in years:
+                for theme in definitions:
+                    prepared.append(
+                        prepare_snapshot(
+                            source_session,
+                            args=args,
+                            boundary=boundary,
+                            boundary_metric=boundary_metric,
+                            year=year,
+                            theme=theme,
+                        )
+                    )
+
+        if args.fetch_only:
+            results = [
+                {
+                    "year": item.year,
+                    "theme": item.theme.key,
+                    "status": "prepared",
+                    "feature_count": item.feature_count,
+                    "raster_path": str(item.raster_path),
+                    "vector_path": str(item.vector_path),
+                    "manifest_path": str(item.manifest_path),
+                }
+                for item in prepared
+            ]
+        else:
+            base_url = args.base_url.rstrip("/")
+            with requests.Session() as api_session:
+                project_id, area_id, existing = locate_workspace(api_session, base_url, args)
+                for item in prepared:
+                    key = series_key(item.theme, args.scope_key)
+                    observed_date = f"{item.year}-01-01"
+                    dataset = next(
+                        (
+                            candidate
+                            for candidate in existing
+                            if candidate.get("temporal_series_key") == key
+                            and str(candidate.get("observed_at") or "").startswith(observed_date)
+                        ),
+                        None,
+                    )
+                    if dataset:
+                        results.append(
+                            {
+                                "year": item.year,
+                                "theme": item.theme.key,
+                                "status": "existing",
+                                "dataset_id": dataset["id"],
+                                "feature_count": dataset.get("feature_count"),
+                            }
+                        )
+                        continue
+                    dataset = upload_snapshot(
+                        api_session,
+                        base_url=base_url,
+                        project_id=project_id,
+                        area_id=area_id,
+                        args=args,
+                        snapshot=item,
+                    )
+                    existing.append(dataset)
+                    results.append(
+                        {
+                            "year": item.year,
+                            "theme": item.theme.key,
+                            "status": "imported",
+                            "dataset_id": dataset["id"],
+                            "feature_count": dataset.get("feature_count"),
+                        }
+                    )
+    except (OSError, RuntimeError, requests.RequestException, ValueError, KeyError, json.JSONDecodeError) as exc:
+        print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr)
+        return 1
+
+    print(
+        json.dumps(
+            {
+                "status": "ok",
+                "scope": args.scope_key,
+                "municipality": args.municipality_name,
+                "series": [series_key(theme, args.scope_key) for theme in definitions],
+                "snapshots": results,
+            },
+            ensure_ascii=False,
+            indent=2,
+        )
+    )
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

scripts/run_readiness_check.sh

@@ -44,6 +44,9 @@ ${PYTHON_BIN} -m py_compile backend/scripts/yolo_preflight.py
 ${PYTHON_BIN} -m py_compile scripts/prepare_operator_real_data_samples.py
 ${PYTHON_BIN} -m py_compile scripts/provision_mol_municipality_workspace.py
 ${PYTHON_BIN} -m py_compile scripts/provision_mol_context_layers.py
+${PYTHON_BIN} -m py_compile scripts/provision_mol_population_history.py
+${PYTHON_BIN} -m py_compile scripts/provision_mol_historical_landuse.py
+${PYTHON_BIN} -m py_compile scripts/provision_official_landuse_timeseries.py
 ${PYTHON_BIN} -m py_compile scripts/export_operator_yolo_dataset.py
 ${PYTHON_BIN} -m py_compile scripts/export_operator_yolo_tile_dataset.py
 ${PYTHON_BIN} -m py_compile scripts/audit_operator_yolo_dataset_quality.py