glam/scripts/download_dutch_municipalities_geojson.py

#!/usr/bin/env python3
"""
Download Dutch municipality (gemeente) boundaries from PDOK Kadaster API.

Source: Kadaster Bestuurlijke Gebieden (most up-to-date administrative boundaries)
API: https://api.pdok.nl/kadaster/bestuurlijkegebieden/ogc/v1/

Output: frontend/public/data/netherlands_municipalities.geojson
"""

import json
import requests
from pathlib import Path

# PDOK Kadaster Bestuurlijke Gebieden API
BASE_URL = "https://api.pdok.nl/kadaster/bestuurlijkegebieden/ogc/v1"
COLLECTION = "gemeentegebied"

# Output path
OUTPUT_DIR = Path(__file__).parent.parent / "frontend" / "public" / "data"
OUTPUT_FILE = OUTPUT_DIR / "netherlands_municipalities.geojson"


def fetch_all_municipalities():
    """Fetch all municipality boundaries with cursor-based pagination."""
    all_features = []
    limit = 100  # API limit per request
    
    print("Fetching Dutch municipality boundaries from PDOK Kadaster API...")
    
    # Initial URL
    url = f"{BASE_URL}/collections/{COLLECTION}/items?f=json&limit={limit}"
    
    while url:
        response = requests.get(url)
        response.raise_for_status()
        
        data = response.json()
        features = data.get("features", [])
        
        if not features:
            break
            
        all_features.extend(features)
        print(f"  Fetched {len(all_features)} municipalities...")
        
        # Find the "next" link for cursor-based pagination
        next_url = None
        for link in data.get("links", []):
            if link.get("rel") == "next":
                next_url = link.get("href")
                break
        
        url = next_url
    
    return all_features


def simplify_coordinates(coords, tolerance=0.0001):
    """
    Simplify coordinates using Douglas-Peucker algorithm (simplified version).
    This reduces file size by removing redundant points.
    
    For werkgebied display, we don't need sub-meter precision.
    tolerance=0.0001 degrees ≈ 11 meters
    """
    if len(coords) <= 2:
        return coords
    
    # Simple decimation: keep every Nth point and always keep first/last
    simplified = [coords[0]]
    prev = coords[0]
    
    for i, coord in enumerate(coords[1:-1], 1):
        # Calculate distance from previous kept point
        dx = coord[0] - prev[0]
        dy = coord[1] - prev[1]
        dist_sq = dx * dx + dy * dy
        
        # Keep point if it's far enough from previous
        if dist_sq > tolerance * tolerance:
            simplified.append(coord)
            prev = coord
    
    # Always keep last point
    simplified.append(coords[-1])
    
    return simplified


def simplify_geometry(geometry, tolerance=0.0001):
    """Recursively simplify geometry coordinates."""
    if geometry is None:
        return None
    
    geom_type = geometry.get("type")
    coords = geometry.get("coordinates")
    
    if geom_type == "Polygon":
        # Simplify each ring
        new_coords = []
        for ring in coords:
            simplified_ring = simplify_coordinates(ring, tolerance)
            # Polygon rings need at least 4 points
            if len(simplified_ring) >= 4:
                new_coords.append(simplified_ring)
        return {"type": "Polygon", "coordinates": new_coords} if new_coords else None
    
    elif geom_type == "MultiPolygon":
        # Simplify each polygon
        new_coords = []
        for polygon in coords:
            simplified_polygon = []
            for ring in polygon:
                simplified_ring = simplify_coordinates(ring, tolerance)
                if len(simplified_ring) >= 4:
                    simplified_polygon.append(simplified_ring)
            if simplified_polygon:
                new_coords.append(simplified_polygon)
        return {"type": "MultiPolygon", "coordinates": new_coords} if new_coords else None
    
    return geometry


def simplify_properties(features, simplify_geom=True):
    """
    Simplify feature properties to reduce file size.
    Keep only essential properties for werkgebied display.
    """
    simplified = []
    
    for feature in features:
        props = feature.get("properties", {})
        
        # Extract essential properties (note: API uses snake_case)
        simplified_props = {
            "code": props.get("code"),  # Municipality code (e.g., "0363" for Amsterdam)
            "naam": props.get("naam"),  # Municipality name
            "provincieCode": props.get("ligt_in_provincie_code"),  # Province code
            "provincieNaam": props.get("ligt_in_provincie_naam"),  # Province name
        }
        
        # Simplify geometry if requested
        geometry = feature.get("geometry")
        if simplify_geom and geometry:
            geometry = simplify_geometry(geometry, tolerance=0.0002)  # ~22m precision
        
        if geometry:  # Only add if geometry is valid
            simplified.append({
                "type": "Feature",
                "properties": simplified_props,
                "geometry": geometry
            })
    
    return simplified


def create_geojson(features):
    """Create a GeoJSON FeatureCollection."""
    return {
        "type": "FeatureCollection",
        "name": "netherlands_municipalities",
        "crs": {
            "type": "name",
            "properties": {
                "name": "urn:ogc:def:crs:OGC:1.3:CRS84"
            }
        },
        "features": features
    }


def main():
    # Ensure output directory exists
    OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
    
    # Fetch all municipalities
    features = fetch_all_municipalities()
    print(f"\nTotal municipalities fetched: {len(features)}")
    
    # Simplify properties and geometry
    simplified_features = simplify_properties(features, simplify_geom=True)
    
    # Create GeoJSON
    geojson = create_geojson(simplified_features)
    
    # Write to file
    with open(OUTPUT_FILE, "w", encoding="utf-8") as f:
        json.dump(geojson, f, ensure_ascii=False, indent=2)
    
    file_size_mb = OUTPUT_FILE.stat().st_size / (1024 * 1024)
    print(f"\nSaved to: {OUTPUT_FILE}")
    print(f"File size: {file_size_mb:.2f} MB")
    
    # Print sample municipalities
    print("\nSample municipalities:")
    for feature in simplified_features[:5]:
        props = feature["properties"]
        print(f"  - {props['naam']} (code: {props['code']}, province: {props.get('provincieNaam', 'N/A')})")


if __name__ == "__main__":
    main()