Bewerkte lagen samenstellen met Python

Gebruik deze voorbeelden om zelf bronnen, instellingen en bewerkingen te combineren. Elk script toont de stappen met publieke functies onder de productie-API. Volg eerst Aan de slag voor installatie en uitvoering.

Stel bovenaan het script DATA_DIR, AREA en EXPERIMENT in. Downloads worden hergebruikt. Een eigen bestand gebruiken kan door de downloadaanroep te vervangen door het betreffende pad; laag- en veldnamen moeten aansluiten op de API-referentie. Tussenresultaten blijven beschikbaar om te bekijken of verder te bewerken.

De uitvoer staat onder data/processed_data/voorbeelden/alkmaar_centrum/<experiment>/. Kies bij andere bronnen of instellingen een nieuwe experimentnaam. Losse functies kunnen eerder gemaakte bestanden hergebruiken zonder de gewijzigde instellingen toe te passen. Dit is geen automatisch beheerde productierun.

Landelijke bronbestanden blijven groot, ook voor Alkmaar-centrum. Bij de rastervoorbeelden zijn de instellingen onderdeel van de methode: een andere buffer, resolutie of indeling kan andere resultaten geven dan de standaardproductie.

Afwateringseenheden

Lees de HYDAMO-objecten, selecteer primaire en secundaire watergangen en bereid de segmenten voor met prepare_watersysteem. Daarna berekent calculate_afwateringseenheden_tiles de afwateringseenheden.

Instelbaar zijn onder meer de verbindingstolerantie, segmentlengte, tegelgrootte, buffer, rasterresolutie en inbranddiepte. Het voorbeeld gebruikt een kerntegel van 1 km met 2 km buffer; het samengevoegde resultaat wordt op Alkmaar-centrum geknipt. AHN en watersysteem worden ook voor de omgeving opgehaald. HYDAMO blijft een landelijke download. Gebruik de omgeving met PCRaster.

Uitvoer: watersysteem.gpkg, tegelbestanden en afwateringseenheden.gpkg. Bekijk ook de gemelde randproblemen en resterende gaten in het resultaat.

Betekenis van de laag · API-referentie

Download het Python-script.

"""Watersysteem voorbereiden en afwateringseenheden met eigen instellingen berekenen."""

from math import ceil, floor
from pathlib import Path

import geopandas as gpd
from shapely.geometry import box

from waterlagen.administratieve_gebieden import (
    DEFAULT_BESTUURLIJKE_GEBIEDEN_YEAR,
    download_bestuurlijke_gebieden,
)
from waterlagen.afwateringseenheden import (
    calculate_afwateringseenheden_tiles,
    prepare_watersysteem,
    read_hydroobjecten,
    read_puntobjecten,
    write_watersysteem,
)
from waterlagen.ahn import download_ahn
from waterlagen.areas import NAMED_AREAS, NAMED_AREAS_CRS
from waterlagen.datastore import DataStore
from waterlagen.hydamo import download_hydamo
from waterlagen.logger import production_logging
from waterlagen.settings import settings

DATA_DIR = Path("data")
AREA = "alkmaar_centrum"
EXPERIMENT = "experiment_01"  # Nieuwe naam bij andere bronnen of rekeninstellingen.
TILE_SIZE_M = 1000  # Eén kerntegel voor dit vierkante voorbeeldgebied.
TILE_BUFFER_M = 2000
RESOLUTION_M = 2.0
TOLERANCE_M = 2.0
MAX_SEGMENT_LENGTH_M = 500.0
BURN_DEPTH_M = 100.0
MAX_FILL_DEPTH_M = 50.0
WORKERS = 1

if __name__ == "__main__":
    store = DataStore(
        data_dir=DATA_DIR,
        source_data_dir=DATA_DIR / "source_data",
        processed_data_dir=DATA_DIR / "processed_data",
    )
    output_dir = (
        store.processed_data_dir
        / "voorbeelden"
        / AREA
        / EXPERIMENT
        / "afwateringseenheden"
    )
    gebied = (
        gpd.GeoSeries([NAMED_AREAS[AREA]], crs=NAMED_AREAS_CRS)
        .to_crs(settings.crs)
        .iloc[0]
    )
    # Voor dit voorbeeld ligt de tegel-oorsprong op de linkeronderhoek van het gebied.
    # Rond de dekking naar volledige tegels af voordat de buffer wordt toegevoegd.
    xmin, ymin, xmax, ymax = gebied.bounds
    origin_x, origin_y = floor(xmin), floor(ymin)
    tile_bounds = box(
        origin_x,
        origin_y,
        origin_x + ceil((xmax - origin_x) / TILE_SIZE_M) * TILE_SIZE_M,
        origin_y + ceil((ymax - origin_y) / TILE_SIZE_M) * TILE_SIZE_M,
    )
    context = tile_bounds.buffer(TILE_BUFFER_M).envelope

    with production_logging():
        # 1. Bronnen ophalen. HYDAMO en de landsgrens zijn landelijke bestanden.
        hydamo = download_hydamo(download_dir=store.hydamo_dir, overwrite=False)
        landsgrens = download_bestuurlijke_gebieden(
            year=DEFAULT_BESTUURLIJKE_GEBIEDEN_YEAR,
            download_dir=store.administratieve_gebieden_dir,
            overwrite=False,
        )
        ahn = download_ahn(
            ahn_dir=store.ahn_dir,
            poly_mask=context,
            model="dtm",
            ahn_version=4,
            service="ahn_pdok",
            cell_size="05",
            missing_only=True,
            create_vrt=True,
        )

        # 2. Lees het watersysteem inclusief omgeving; pas de selectie hier aan.
        hydro = read_hydroobjecten(hydamo.target_path, spatial_selection=context)
        primair = hydro.loc[hydro.categorieoppwaterlichaam.eq("primair")].copy()
        secundair = hydro.loc[hydro.categorieoppwaterlichaam.eq("secundair")].copy()
        punten = read_puntobjecten(
            hydamo.target_path, layers=["gemaal", "stuw"], spatial_selection=context
        )
        netwerk = prepare_watersysteem(
            primair,
            punten,
            hydroobject_secundair=secundair,
            tolerance=TOLERANCE_M,
            max_length=MAX_SEGMENT_LENGTH_M,
        )
        watersysteem = write_watersysteem(
            hydroobject_primair=primair,
            hydroobject_secundair=secundair,
            watersysteem=netwerk,
            output_path=output_dir / "watersysteem.gpkg",
            secondary_connection_tolerance=TOLERANCE_M,
            overwrite=False,
        )

        # 3. Bereken de tegels en voeg ze samen. Pas parameters bovenaan aan.
        resultaat = calculate_afwateringseenheden_tiles(
            gebied,
            ahn_vrt_path=ahn,
            landsgrens_path=landsgrens.target_path,
            watersysteem_path=watersysteem,
            output_dir=output_dir / "tiles",
            merged_output_path=output_dir / "afwateringseenheden.gpkg",
            tile_size_m=TILE_SIZE_M,
            tile_buffer_m=TILE_BUFFER_M,
            origin_x=origin_x,
            origin_y=origin_y,
            resolution_m=RESOLUTION_M,
            burn_depth_m=BURN_DEPTH_M,
            max_fill_depth_m=MAX_FILL_DEPTH_M,
            workers=WORKERS,
            random_seed=12345,
            clip_to_area=True,
            overwrite=False,
        )
        print(resultaat.merged_path)
        print("Tegels met randproblemen:", resultaat.boundary_issue_tile_ids)

Auto's

Haal BAG en CBS op en koppel de woon-VBO's aan buurten met bouw_vbo_buurt. Die voorbereiding is landelijk en wordt ook door het inwonersvoorbeeld gebruikt. HERBOUW_VOORBEREIDING=True vernieuwt de koppeling na een bronwijziging; het downloadt geen nieuwe bronversies.

Selecteer daarna de complete buurten die Alkmaar-centrum overlappen en behoud alle bijbehorende woon-VBO's, ook buiten het vierkant. Geef die GeoDataFrames aan bereken_personenautos_per_vbo. Het resultaat kan direct worden bekeken, met andere gegevens worden gecombineerd of worden opgeslagen in autos.gpkg. Verklein de VBO-selectie pas na de berekening; anders ontbreekt buurtcontext.

Betekenis van de laag · API-referentie · Gedeelde voorbereiding

Download het Python-script.

"""autos verdelen en het resultaat zelf verder gebruiken."""

from pathlib import Path

import geopandas as gpd

from waterlagen._geopandas import read_file
from waterlagen.administratieve_gebieden import download_wijk_buurtkaart_2025
from waterlagen.areas import NAMED_AREAS, NAMED_AREAS_CRS
from waterlagen.bag import download_bag_light
from waterlagen.cbs import download_buurtgegevens_2025
from waterlagen.datastore import DataStore
from waterlagen.autos import bereken_personenautos_per_vbo
from waterlagen.logger import production_logging
from waterlagen.vbo_buurt import bouw_vbo_buurt, sorteer_op_hilbert

# Instellingen: kies eigen bronpaden in de downloadstappen of hergebruik deze cache.
DATA_DIR = Path("data")
AREA = "alkmaar_centrum"
EXPERIMENT = "experiment_01"  # Kies een nieuwe naam voor een andere proef.
HERBOUW_VOORBEREIDING = False  # True na het vervangen van BAG/CBS-bronnen.

if __name__ == "__main__":
    store = DataStore(
        data_dir=DATA_DIR,
        source_data_dir=DATA_DIR / "source_data",
        processed_data_dir=DATA_DIR / "processed_data",
    )
    output_dir = store.processed_data_dir / "voorbeelden" / AREA / EXPERIMENT
    target = output_dir / "autos.gpkg"
    if target.exists():
        raise FileExistsError(f"Kies een nieuwe experimentnaam: {target} bestaat al")

    with production_logging():
        # 1. Bronnen ophalen; deze downloads en de koppeling zijn landelijk.
        bag = download_bag_light(download_dir=store.bag_dir, overwrite=False)
        kaart = download_wijk_buurtkaart_2025(
            download_dir=store.administratieve_gebieden_dir, overwrite=False
        )
        cijfers = download_buurtgegevens_2025(
            download_dir=store.cbs_dir, overwrite=False
        )
        koppeling = bouw_vbo_buurt(
            bag_path=bag,
            buurtkaart_path=kaart.target_path,
            cbs_buurtgegevens_path=cijfers.target_path,
            bag_vbo_path=store.bag_vbo_path,
            cbs_buurt_path=store.cbs_buurt_path,
            overwrite=HERBOUW_VOORBEREIDING,
        )

        # 2. Kies complete buurten en behoud ALLE bijbehorende woon-VBO's.
        buurten = read_file(koppeling.cbs_buurt_path, layer="cbs_buurt")
        gebied = (
            gpd.GeoSeries([NAMED_AREAS[AREA]], crs=NAMED_AREAS_CRS)
            .to_crs(buurten.crs)
            .iloc[0]
        )
        buurten = buurten.loc[
            buurten.intersects(gebied) & ~buurten.touches(gebied)
        ].copy()
        if buurten.empty:
            raise ValueError("Geen CBS-buurten in het gekozen gebied")
        vbos = read_file(koppeling.bag_vbo_path, layer="bag_vbo")
        vbos = sorteer_op_hilbert(
            vbos.loc[vbos.buurtcode.isin(buurten.buurtcode)].copy()
        )

        # 3. Bereken met GeoDataFrames; resultaat is beschikbaar voor eigen analyses.
        resultaat, ontbrekende_waarden, ongecontroleerde_buurten = (
            bereken_personenautos_per_vbo(vbos, buurten)
        )
        print(resultaat.head())
        print(f"VBO's met ontbrekende CBS-waarden: {ontbrekende_waarden}")
        print(f"Buurten zonder somcontrole: {ongecontroleerde_buurten}")

        # 4. Voeg hier eigen analyses toe. Schrijf daarna het gekozen resultaat.
        output_dir.mkdir(parents=True, exist_ok=True)
        tijdelijk = target.with_name("autos.tmp.gpkg")
        resultaat.to_file(tijdelijk, layer="autos", driver="GPKG", index=False)
        tijdelijk.replace(target)
        print(target)

Hoogtemodel

Voer eerst het landgebruikvoorbeeld uit. Kies met LANDUSE_EXPERIMENT welke uitsnede als invoer dient. Het bronraster, de gebouw-ID's en de voorbereide pandgeometrieën moeten bij die uitsnede blijven.

Stel met DemConfig de interpolatieafstand, gebouwzoekafstanden en het percentiel in. De maximale gebouwzoekafstand mag niet groter zijn dan de BUILDING_CONTEXT_M waarmee het landgebruik is voorbereid. bouw_dem_tiles berekent het DEM op hetzelfde rasterrooster. Het originele AHN en het landgebruik moeten dezelfde resolutie en pixeluitlijning hebben; er wordt niet automatisch herschaald.

Het voorbeeld haalt AHN op voor de uitsnede én de omgeving van volledige panden. Vervang het VRT-pad in de bouwstap om eigen AHN-invoer te gebruiken. Uitvoer: dem/dem.tif, de floatversie, bronrasters en controlebestanden.

Betekenis van de laag · API-referentie

Download het Python-script.

"""Een DEM maken uit eigen AHN-invoer en een gekozen landgebruikuitsnede."""

from pathlib import Path

import geopandas as gpd
from shapely.geometry import box

from waterlagen._geopandas import read_file
from waterlagen.administratieve_gebieden import (
    DEFAULT_BESTUURLIJKE_GEBIEDEN_YEAR,
    download_bestuurlijke_gebieden,
)
from waterlagen.ahn import download_ahn
from waterlagen.areas import NAMED_AREAS, NAMED_AREAS_CRS
from waterlagen.datastore import DataStore
from waterlagen.dem import DemConfig, bouw_dem_tiles
from waterlagen.functioneel_landgebruik.gebouwen import building_paths
from waterlagen.logger import production_logging
from waterlagen.settings import settings

DATA_DIR = Path("data")
AREA = "alkmaar_centrum"
LANDUSE_EXPERIMENT = "experiment_01"  # Voer eerst het landgebruikvoorbeeld uit.
EXPERIMENT = "experiment_01"  # Nieuwe naam bij andere bronnen of rekeninstellingen.
CONFIG = DemConfig(
    interpolation_max_distance_m=250.0,
    building_initial_buffer_m=1.0,
    building_buffer_step_m=1.0,
    building_max_search_distance_m=5.0,
    building_percentile=75.0,
    building_offset_m=0.05,
)
WORKERS = 1

if __name__ == "__main__":
    store = DataStore(
        data_dir=DATA_DIR,
        source_data_dir=DATA_DIR / "source_data",
        processed_data_dir=DATA_DIR / "processed_data",
    )
    example_dir = store.processed_data_dir / "voorbeelden" / AREA
    landuse_tile = (
        example_dir
        / LANDUSE_EXPERIMENT
        / "functioneel_landgebruik"
        / "functioneel_landgebruik.tif"
    )

    with production_logging():
        # 1. Gebruik een landgebruikuitsnede met bronraster, gebouw-ID's en pandgeometrieën.
        # De resolutie en pixeluitlijning moeten overeenkomen met het originele AHN.
        _, gebouwen_path = building_paths(landuse_tile)
        gebouwen = read_file(gebouwen_path, layer="gebouwen").to_crs(settings.crs)
        gebied = (
            gpd.GeoSeries([NAMED_AREAS[AREA]], crs=NAMED_AREAS_CRS)
            .to_crs(settings.crs)
            .iloc[0]
        )
        context = gebied
        if not gebouwen.empty:
            context = context.union(box(*gebouwen.total_bounds))
        context = context.buffer(
            max(
                CONFIG.interpolation_max_distance_m,
                CONFIG.building_max_search_distance_m,
            )
        )

        # 2. Haal AHN op voor de uitsnede, volledige panden en zoekafstanden.
        ahn = download_ahn(
            ahn_dir=store.ahn_dir,
            poly_mask=context,
            model="dtm",
            ahn_version=4,
            service="ahn_pdok",
            cell_size="05",
            missing_only=True,
            create_vrt=True,
        )
        landsgrens = download_bestuurlijke_gebieden(
            year=DEFAULT_BESTUURLIJKE_GEBIEDEN_YEAR,
            download_dir=store.administratieve_gebieden_dir,
            overwrite=False,
        )

        # 3. Bereken met de instellingen bovenaan. Eigen AHN? Geef hier het VRT-pad op.
        resultaat = bouw_dem_tiles(
            example_dir / EXPERIMENT / "dem",
            ahn_vrt_path=ahn,
            landuse_tiles=[landuse_tile],
            config=CONFIG,
            workers=WORKERS,
            building_workers=WORKERS,
            landgebied_path=landsgrens.target_path,
            overwrite=False,
        )
        print(resultaat)

Functioneel landgebruik

De downloadstappen tonen alle gebruikte bronnen. Daarna worden LIWO-vlakken en een index van BAG-panden voorbereid. Met FunctioneelLandgebruikSources kiest u expliciet de invoer voor bouw_functioneel_landgebruik.

Pas de resolutie, aanvulafstand of classificatietabel aan. Het voorbeeld maakt een rechthoek van 1 × 1 km; bij 0,5 m zijn dat 4 miljoen rastercellen, plus tijdelijke rekenruimte. De downloads en een deel van de bronvoorbereiding blijven landelijk. De BGT-download bevat hier de zes benodigde objecttypen; het kleine BGT-downloadvoorbeeld met alleen panden en water is hiervoor onvoldoende.

De pandindex en BUILDING_CONTEXT_M zorgen dat het resultaat ook voor het DEM kan worden gebruikt. Houd dan de resolutie op 0,5 m. Met landgebied_path wordt de landsgrens voor het aanvullen expliciet meegegeven.

Uitvoer: functioneel_landgebruik/functioneel_landgebruik.tif, met bronraster, gebouwgegevens, QGIS-stijl en controlebestand. Dit is een losse rasteruitsnede.

Betekenis van de laag · API-referentie

Download het Python-script.

"""Een instelbare uitsnede van functioneel landgebruik maken, inclusief gebouwen."""

from pathlib import Path

import geopandas as gpd

from waterlagen.administratieve_gebieden import (
    DEFAULT_BESTUURLIJKE_GEBIEDEN_YEAR,
    download_bestuurlijke_gebieden,
)
from waterlagen.areas import NAMED_AREAS, NAMED_AREAS_CRS
from waterlagen.bag import download_bag_light
from waterlagen.bgt import download_bgt
from waterlagen.brp import download_brp
from waterlagen.datastore import DataStore
from waterlagen.functioneel_landgebruik import (
    FunctioneelLandgebruikSources,
    bouw_functioneel_landgebruik,
)
from waterlagen.functioneel_landgebruik.gebouwen import ensure_building_index
from waterlagen.functioneel_landgebruik.osm_drinkwater import download_osm_drinkwater
from waterlagen.hydamo import download_hydamo
from waterlagen.liwo import FILENAME, download_raster, polygonize_selected
from waterlagen.logger import production_logging
from waterlagen.settings import settings
from waterlagen.top10nl import download_top10nl
from waterlagen.waterketen_damo import download_waterketen_damo

DATA_DIR = Path("data")
AREA = "alkmaar_centrum"
EXPERIMENT = "experiment_01"  # Nieuwe naam bij andere bronnen of instellingen.
RESOLUTION_M = 0.5  # Houd 0.5 aan als deze uitsnede daarna naar het DEM gaat.
GAP_FILL_DISTANCE_M = 1.0
BUILDING_CONTEXT_M = 5.0  # Minstens de maximale gebouwzoekafstand van het DEM.
MAPPING_CSV = None  # Of Path("mijn_landgebruik.csv"), met het bestaande tabelschema.
BGT_FEATURES = (
    "waterdeel",
    "ondersteunendwaterdeel",
    "wegdeel",
    "ondersteunendwegdeel",
    "begroeidterreindeel",
    "onbegroeidterreindeel",
)

if __name__ == "__main__":
    store = DataStore(
        data_dir=DATA_DIR,
        source_data_dir=DATA_DIR / "source_data",
        processed_data_dir=DATA_DIR / "processed_data",
    )
    output_dir = (
        store.processed_data_dir
        / "voorbeelden"
        / AREA
        / EXPERIMENT
        / "functioneel_landgebruik"
    )
    gebied = (
        gpd.GeoSeries([NAMED_AREAS[AREA]], crs=NAMED_AREAS_CRS)
        .to_crs(settings.crs)
        .iloc[0]
    )

    with production_logging():
        # 1. Bronnen ophalen/hergebruiken. Dit zijn landelijke downloads!
        # Vervang een downloadaanroep door een eigen Path om lokale invoer te kiezen.
        bag = download_bag_light(download_dir=store.bag_dir, overwrite=False)
        bgt = download_bgt(
            download_dir=store.bgt_dir,
            target_path=store.bgt_dir / "voorbeeld_landgebruik.gpkg",
            featuretypes=BGT_FEATURES,
            overwrite=False,
        )
        brp = download_brp(
            download_dir=store.brp_dir, filename="brp.gpkg", overwrite=False
        )
        top10nl = download_top10nl(download_dir=store.top10nl_dir, overwrite=False)
        liwo = download_raster(
            store.source_data_dir / "liwo" / FILENAME, overwrite=False
        )
        hydamo = download_hydamo(download_dir=store.hydamo_dir, overwrite=False)
        waterketen = download_waterketen_damo(
            download_dir=store.waterketen_damo_dir, overwrite=False
        )
        drinkwater = download_osm_drinkwater(data_store=store, overwrite=False)
        landsgrens = download_bestuurlijke_gebieden(
            year=DEFAULT_BESTUURLIJKE_GEBIEDEN_YEAR,
            download_dir=store.administratieve_gebieden_dir,
            overwrite=False,
        )

        # 2. Maak de LIWO-vlakken en een herbruikbare koppeling van BAG-ID's.
        buitendijks = output_dir / "buitendijks.gpkg"
        if not buitendijks.exists():
            polygonize_selected(liwo, buitendijks)
        gebouw_index = ensure_building_index(bag, output_dir / "gebouw_index.sqlite")

        # 3. Kies expliciet de invoer. Eigen bestanden moeten dezelfde lagen/velden hebben.
        sources = FunctioneelLandgebruikSources(
            bag_gpkg=bag,
            bgt_gpkg=bgt.target_path,
            brp_gpkg=brp.path,
            top10nl_gpkg=top10nl,
            buitendijks_gpkg=buitendijks,
            gemalen_gpkg=hydamo.target_path,
            rwzi_gpkg=waterketen.target_path,
            drinking_water_gpkg=drinkwater,
        )

        # 4. Bereken precies deze rechthoek; de bronvoorbereiding gebruikt ook de omgeving.
        resultaat = bouw_functioneel_landgebruik(
            output_dir / "functioneel_landgebruik.tif",
            bounds=gebied.bounds,
            crs=settings.crs,
            resolution_m=RESOLUTION_M,
            gap_fill_distance_m=GAP_FILL_DISTANCE_M,
            sources=sources,
            landgebied_path=landsgrens.target_path,
            mapping_csv=MAPPING_CSV,
            building_index_path=gebouw_index,
            building_context_m=BUILDING_CONTEXT_M,
            diagnostics_path=output_dir / "nodata.gpkg",
            download_missing_sources=False,
            overwrite=False,
        )
        print(resultaat)

Inwoners

Dit voorbeeld gebruikt dezelfde landelijke BAG/CBS-koppeling en selectie van complete buurten als Auto's. Met bereken_inwoners_per_vbo worden beide verdeelmethoden berekend. De uitkomst is een GeoDataFrame: bekijk de kolommen, voeg eigen analyses toe en kies zelf wat u opslaat. Het voorbeeld schrijft inwoners.gpkg.

Betekenis van de laag · API-referentie · Gedeelde voorbereiding

Download het Python-script.

"""inwoners verdelen en het resultaat zelf verder gebruiken."""

from pathlib import Path

import geopandas as gpd

from waterlagen._geopandas import read_file
from waterlagen.administratieve_gebieden import download_wijk_buurtkaart_2025
from waterlagen.areas import NAMED_AREAS, NAMED_AREAS_CRS
from waterlagen.bag import download_bag_light
from waterlagen.cbs import download_buurtgegevens_2025
from waterlagen.datastore import DataStore
from waterlagen.inwoners import bereken_inwoners_per_vbo
from waterlagen.logger import production_logging
from waterlagen.vbo_buurt import bouw_vbo_buurt, sorteer_op_hilbert

# Instellingen: kies eigen bronpaden in de downloadstappen of hergebruik deze cache.
DATA_DIR = Path("data")
AREA = "alkmaar_centrum"
EXPERIMENT = "experiment_01"  # Kies een nieuwe naam voor een andere proef.
HERBOUW_VOORBEREIDING = False  # True na het vervangen van BAG/CBS-bronnen.

if __name__ == "__main__":
    store = DataStore(
        data_dir=DATA_DIR,
        source_data_dir=DATA_DIR / "source_data",
        processed_data_dir=DATA_DIR / "processed_data",
    )
    output_dir = store.processed_data_dir / "voorbeelden" / AREA / EXPERIMENT
    target = output_dir / "inwoners.gpkg"
    if target.exists():
        raise FileExistsError(f"Kies een nieuwe experimentnaam: {target} bestaat al")

    with production_logging():
        # 1. Bronnen ophalen; deze downloads en de koppeling zijn landelijk.
        bag = download_bag_light(download_dir=store.bag_dir, overwrite=False)
        kaart = download_wijk_buurtkaart_2025(
            download_dir=store.administratieve_gebieden_dir, overwrite=False
        )
        cijfers = download_buurtgegevens_2025(
            download_dir=store.cbs_dir, overwrite=False
        )
        koppeling = bouw_vbo_buurt(
            bag_path=bag,
            buurtkaart_path=kaart.target_path,
            cbs_buurtgegevens_path=cijfers.target_path,
            bag_vbo_path=store.bag_vbo_path,
            cbs_buurt_path=store.cbs_buurt_path,
            overwrite=HERBOUW_VOORBEREIDING,
        )

        # 2. Kies complete buurten en behoud ALLE bijbehorende woon-VBO's.
        buurten = read_file(koppeling.cbs_buurt_path, layer="cbs_buurt")
        gebied = (
            gpd.GeoSeries([NAMED_AREAS[AREA]], crs=NAMED_AREAS_CRS)
            .to_crs(buurten.crs)
            .iloc[0]
        )
        buurten = buurten.loc[
            buurten.intersects(gebied) & ~buurten.touches(gebied)
        ].copy()
        if buurten.empty:
            raise ValueError("Geen CBS-buurten in het gekozen gebied")
        vbos = read_file(koppeling.bag_vbo_path, layer="bag_vbo")
        vbos = sorteer_op_hilbert(
            vbos.loc[vbos.buurtcode.isin(buurten.buurtcode)].copy()
        )

        # 3. Bereken met GeoDataFrames; resultaat is beschikbaar voor eigen analyses.
        resultaat, ontbrekende_waarden, ongecontroleerde_buurten = (
            bereken_inwoners_per_vbo(vbos, buurten)
        )
        print(resultaat.head())
        print(f"VBO's met ontbrekende CBS-waarden: {ontbrekende_waarden}")
        print(f"Buurten zonder somcontrole: {ongecontroleerde_buurten}")

        # 4. Voeg hier eigen analyses toe. Schrijf daarna het gekozen resultaat.
        output_dir.mkdir(parents=True, exist_ok=True)
        tijdelijk = target.with_name("inwoners.tmp.gpkg")
        resultaat.to_file(tijdelijk, layer="inwoners", driver="GPKG", index=False)
        tijdelijk.replace(target)
        print(target)

Redactie: AI