feat: operationalize Flemish land and nature themes

backend/app/api/routes/datasets.py

@@ -32,6 +32,7 @@ from app.schemas import (
     ThematicRasterAcquireRequest,
     ThematicRasterSelectionRequest,
     GrbAcquireRequest,
+    OfficialVectorAcquireRequest,
     VectorBBoxResponse,
     VectorBufferRequest,
     VectorClipRequest,
@@ -51,6 +52,7 @@ from app.services.source_freshness_service import SourceFreshnessService
 from app.services.source_catalog_probe_service import SourceCatalogProbeService
 from app.services.grb_refresh_plan_service import GrbRefreshPlanService
 from app.services.grb_acquisition_service import GrbAcquisitionService
+from app.services.official_vector_acquisition_service import OfficialVectorAcquisitionService
 from app.services.orthophoto_acquisition_service import OrthophotoAcquisitionService
 from app.services.dhmv_acquisition_service import DhmvAcquisitionService
 from app.services.terrain_analysis_service import TerrainAnalysisService
@@ -219,6 +221,30 @@ def list_grb_products(project_id: UUID, db: Session = Depends(get_db)):
     return envelope({"items": items, "total": len(items)})
 
 
+@router.post("/datasets/official-vector/acquire", response_model=dict)
+def acquire_bounded_official_vector(
+    project_id: UUID,
+    payload: OfficialVectorAcquireRequest,
+    db: Session = Depends(get_db),
+):
+    job = JobService.run_sync_job(
+        db=db,
+        project_id=project_id,
+        job_type="vector.official.acquire",
+        parameters=payload.model_dump(mode="json"),
+        operation=lambda: OfficialVectorAcquisitionService.acquire(db, project_id, payload),
+    )
+    return envelope(job)
+
+
+@router.get("/datasets/official-vector/products", response_model=dict)
+def list_official_vector_products(project_id: UUID, db: Session = Depends(get_db)):
+    if not db.get(Project, project_id):
+        raise AppError(code="PROJECT_NOT_FOUND", message="Project not found", status_code=404)
+    items = OfficialVectorAcquisitionService.list_products()
+    return envelope({"items": items, "total": len(items)})
+
+
 @router.post("/datasets/flood-hazard/acquire", response_model=dict)
 def acquire_bounded_flood_hazard(
     project_id: UUID,

backend/app/core/config.py

@@ -87,6 +87,66 @@ class Settings(BaseSettings):
         validation_alias="GRB_MAX_TOTAL_RESPONSE_MB",
     )
     grb_cache_ttl_hours: int = Field(default=24, ge=0, le=8760, validation_alias="GRB_CACHE_TTL_HOURS")
+    official_vector_enabled: bool = Field(default=True, validation_alias="OFFICIAL_VECTOR_ENABLED")
+    bwk_wfs_url: str = Field(
+        default="https://geo.api.vlaanderen.be/BWK/wfs",
+        validation_alias="BWK_WFS_URL",
+    )
+    dov_soil_wfs_url: str = Field(
+        default="https://www.dov.vlaanderen.be/geoserver/wfs",
+        validation_alias="DOV_SOIL_WFS_URL",
+    )
+    official_vector_min_side_m: float = Field(
+        default=10.0,
+        gt=0,
+        validation_alias="OFFICIAL_VECTOR_MIN_SIDE_M",
+    )
+    official_vector_max_side_m: float = Field(
+        default=20_000.0,
+        gt=0,
+        validation_alias="OFFICIAL_VECTOR_MAX_SIDE_M",
+    )
+    official_vector_page_size: int = Field(
+        default=1000,
+        ge=1,
+        le=2000,
+        validation_alias="OFFICIAL_VECTOR_PAGE_SIZE",
+    )
+    official_vector_max_pages: int = Field(
+        default=200,
+        ge=1,
+        le=1000,
+        validation_alias="OFFICIAL_VECTOR_MAX_PAGES",
+    )
+    official_vector_max_features: int = Field(
+        default=100_000,
+        ge=1,
+        validation_alias="OFFICIAL_VECTOR_MAX_FEATURES",
+    )
+    official_vector_timeout_seconds: int = Field(
+        default=180,
+        ge=1,
+        le=600,
+        validation_alias="OFFICIAL_VECTOR_TIMEOUT_SECONDS",
+    )
+    official_vector_max_response_mb: int = Field(
+        default=20,
+        ge=1,
+        le=100,
+        validation_alias="OFFICIAL_VECTOR_MAX_RESPONSE_MB",
+    )
+    official_vector_max_total_response_mb: int = Field(
+        default=256,
+        ge=1,
+        le=2048,
+        validation_alias="OFFICIAL_VECTOR_MAX_TOTAL_RESPONSE_MB",
+    )
+    official_vector_cache_ttl_hours: int = Field(
+        default=24,
+        ge=0,
+        le=8760,
+        validation_alias="OFFICIAL_VECTOR_CACHE_TTL_HOURS",
+    )
     dhmv_enabled: bool = Field(default=True, validation_alias="DHMV_ENABLED")
     dhmv_wcs_url: str = Field(
         default="https://geo.api.vlaanderen.be/DHMV/wcs",

backend/app/schemas/__init__.py

@@ -18,6 +18,11 @@ from .source_catalog import (
 )
 from .grb_refresh import GrbRefreshLayerPlan, GrbRefreshPlan, GrbRefreshPlanSummary
 from .grb import GrbAcquireRequest, GrbAcquisitionResult, GrbProductRead
+from .official_vector import (
+    OfficialVectorAcquireRequest,
+    OfficialVectorAcquisitionResult,
+    OfficialVectorProductRead,
+)
 from .detection import (
     DetectionListResponse,
     DetectionModelCapability,
@@ -165,6 +170,9 @@ __all__ = [
     "GrbAcquireRequest",
     "GrbAcquisitionResult",
     "GrbProductRead",
+    "OfficialVectorAcquireRequest",
+    "OfficialVectorAcquisitionResult",
+    "OfficialVectorProductRead",
     "DetectionListResponse",
     "DetectionModelCapability",
     "DetectionModelsResponse",

backend/app/schemas/official_vector.py

@@ -0,0 +1,54 @@
+from __future__ import annotations
+
+from uuid import UUID
+
+from pydantic import BaseModel
+
+from .operations import VectorSelectionBBox
+
+
+class OfficialVectorAcquireRequest(BaseModel):
+    bbox: VectorSelectionBBox
+    area_id: UUID | None = None
+    product_key: str
+    force_refresh: bool = False
+
+
+class OfficialVectorProductRead(BaseModel):
+    key: str
+    display_name: str
+    theme: str
+    provider: str
+    source_name: str
+    reference_layer_name: str
+    service_type: str
+    collection: str
+    geometry_types: list[str]
+    source_crs: str
+    source_version: str
+    observation_label: str
+    authority_level: str
+    catalog_url: str
+    attribution: str
+    license_note: str
+    limitation_message: str
+
+
+class OfficialVectorAcquisitionResult(BaseModel):
+    output_dataset_id: UUID
+    reused: bool
+    product_key: str
+    display_name: str
+    theme: str
+    provider: str
+    source_name: str
+    reference_layer_name: str
+    service_type: str
+    collection: str
+    feature_count: int
+    candidate_feature_count: int
+    page_count: int
+    bbox_epsg4326: list[float]
+    source_version: str
+    attribution: str
+    limitation_message: str

backend/app/schemas/thematic_raster.py

@@ -30,6 +30,7 @@ class ThematicRasterProductRead(BaseModel):
     license_note: str
     legend_min_label: str
     legend_max_label: str
+    included_source_values: list[int]
     limitation_message: str
 
 

backend/app/services/thematic_raster_acquisition_service.py

@@ -45,6 +45,7 @@ class ThematicRasterProduct:
     legend_min_label: str
     legend_max_label: str
     limitation_message: str
+    included_source_values: tuple[int, ...] = ()
 
 
 class ThematicRasterAcquisitionService:
@@ -100,6 +101,44 @@ class ThematicRasterAcquisitionService:
                     "synoniem met natuur, bos, publieke toegankelijkheid of planologische bestemming."
                 ),
             ),
+            ThematicRasterProduct(
+                key="forest_land_use_2025",
+                display_name="Bos volgens Landgebruik Vlaanderen 2025",
+                theme="forest",
+                metric_kind="binary_area",
+                coverage_id="lu:lu_landgebruik_vlaa_2025_v3",
+                native_resolution_m=10.0,
+                source_value_unit="class_0_1",
+                observation_year=2025,
+                source_version="Toestand 2025 versie 3",
+                catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/landgebruik-vlaanderen-toestand-2025",
+                legend_min_label="Geen bosklasse",
+                legend_max_label="Bos",
+                limitation_message=(
+                    "10 m-afleiding van bronklasse 12 (bos) uit Landgebruik Vlaanderen 2025. De oppervlakte is "
+                    "resolutiegebonden en vormt geen juridische bosgrens, boomtelling, kroonbedekking of houtvolume."
+                ),
+                included_source_values=(12,),
+            ),
+            ThematicRasterProduct(
+                key="agricultural_land_use_2025",
+                display_name="Akker en landbouwgrasland 2025",
+                theme="agriculture",
+                metric_kind="binary_area",
+                coverage_id="lu:lu_landgebruik_vlaa_2025_v3",
+                native_resolution_m=10.0,
+                source_value_unit="class_0_1",
+                observation_year=2025,
+                source_version="Toestand 2025 versie 3",
+                catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/landgebruik-vlaanderen-toestand-2025",
+                legend_min_label="Ander landgebruik",
+                legend_max_label="Akker of landbouwgrasland",
+                limitation_message=(
+                    "10 m-afleiding van bronklassen 13 (akker) en 14 (grasland in landbouwgebruik). Dit is werkelijk "
+                    "landgebruik en geen ALZ-perceelaangifte, eigendomsgrens, teeltregister of juridische bestemming."
+                ),
+                included_source_values=(13, 14),
+            ),
             ThematicRasterProduct(
                 key="population_density_2019",
                 display_name="Inwonersdichtheid per hectare 2019",
@@ -176,6 +215,7 @@ class ThematicRasterAcquisitionService:
                 license_note=ThematicRasterAcquisitionService.LICENSE_NOTE,
                 legend_min_label=product.legend_min_label,
                 legend_max_label=product.legend_max_label,
+                included_source_values=list(product.included_source_values),
                 limitation_message=product.limitation_message,
             ).model_dump()
             for product in ThematicRasterAcquisitionService._products().values()
@@ -461,7 +501,29 @@ class ThematicRasterAcquisitionService:
                 invalid = np.ma.getmaskarray(band) | ~np.isfinite(raw)
                 if source.nodata is not None:
                     invalid |= np.isclose(raw, float(source.nodata))
-                normalized = np.ma.array(raw, mask=invalid)
+                source_values = np.ma.array(raw, mask=invalid).compressed().astype("float64")
+                if product.included_source_values:
+                    rounded = np.rint(source_values)
+                    if not np.allclose(source_values, rounded, atol=0.0001):
+                        raise AppError(
+                            code="THEMATIC_RASTER_INVALID_VALUES",
+                            message="Categorical land-use coverage contains non-integer source classes",
+                            status_code=502,
+                        )
+                    if source_values.size and (
+                        float(source_values.min()) < 0
+                        or float(source_values.max()) > 255
+                    ):
+                        raise AppError(
+                            code="THEMATIC_RASTER_INVALID_VALUES",
+                            message="Categorical land-use coverage contains source classes outside the governed range",
+                            status_code=502,
+                        )
+                    source_classes = np.where(invalid, 0, np.rint(raw)).astype("int16")
+                    binary = np.isin(source_classes, product.included_source_values).astype("float32")
+                    normalized = np.ma.array(binary, mask=invalid)
+                else:
+                    normalized = np.ma.array(raw, mask=invalid)
                 values = normalized.compressed().astype("float64")
                 ThematicRasterAcquisitionService._validate_values(values, product)
                 profile = source.profile.copy()
@@ -482,6 +544,9 @@ class ThematicRasterAcquisitionService:
                     "maximum_value": float(values.max()),
                     "p02_value": float(np.percentile(values, 2)),
                     "p98_value": float(np.percentile(values, 98)),
+                    "included_source_values": list(product.included_source_values),
+                    "source_minimum_value": float(source_values.min()),
+                    "source_maximum_value": float(source_values.max()),
                 }
         except AppError:
             raise
@@ -563,6 +628,7 @@ class ThematicRasterAcquisitionService:
                 "analysis_resolution_m": product.native_resolution_m,
                 "source_crs": ThematicRasterAcquisitionService.SOURCE_CRS,
                 "source_value_unit": product.source_value_unit,
+                "included_source_values": list(product.included_source_values),
                 "observation_year": product.observation_year,
                 "observation_date_precision": "year",
                 "valid_pixel_count": validation["valid_pixel_count"],

backend/app/services/thematic_raster_analysis_service.py

@@ -68,6 +68,10 @@ class ThematicRasterAnalysisService:
     @staticmethod
     def _unsupported_metrics(product: ThematicRasterProduct) -> list[str]:
         if product.metric_kind == "binary_area":
+            if product.theme == "forest":
+                return ["tree_count", "canopy_cover", "timber_volume", "legal_forest_boundary"]
+            if product.theme == "agriculture":
+                return ["declared_parcel_area", "crop_declaration", "ownership", "cadastral_area"]
             return ["object_count", "parcel_area", "current_land_use"]
         if product.metric_kind == "population_density":
             return ["current_population", "household_count", "address_level_population"]
@@ -152,7 +156,12 @@ class ThematicRasterAnalysisService:
             positive_count = int(np.count_nonzero(values >= 0.5))
             positive_area_ha = positive_count * cell_area_m2 / 10_000.0
             positive_share = positive_count / max(1, valid_cell_count) * 100.0
-            label = "Ruimtebeslag" if product.theme == "space_occupation" else "Open ruimte"
+            label = {
+                "space_occupation": "Ruimtebeslag",
+                "open_space": "Open ruimte",
+                "forest": "Bos",
+                "agriculture": "Akker en landbouwgrasland",
+            }[product.theme]
             metrics = [
                 metric(f"{product.theme}_area_ha", f"{label} in selectie", positive_area_ha, "ha", "positive_source_cells_times_cell_area"),
                 metric(f"{product.theme}_share_pct", f"Aandeel {label.lower()}", positive_share, "%", "positive_source_cells_divided_by_valid_selected_cells"),
@@ -216,6 +225,8 @@ class ThematicRasterAnalysisService:
         palettes = {
             "space_occupation": np.asarray([[251, 231, 211], [190, 62, 51]], dtype="float64"),
             "open_space": np.asarray([[221, 238, 219], [38, 122, 70]], dtype="float64"),
+            "forest": np.asarray([[223, 237, 226], [43, 117, 72]], dtype="float64"),
+            "agriculture": np.asarray([[245, 237, 204], [166, 122, 35]], dtype="float64"),
             "population": np.asarray([[238, 231, 246], [103, 58, 151]], dtype="float64"),
             "accessibility": np.asarray([[233, 241, 244], [15, 118, 110]], dtype="float64"),
             "services": np.asarray([[255, 244, 191], [182, 109, 22]], dtype="float64"),