dc-l3-vienna-overpass-historical

Rerun

"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass,
compare to find changes, and output vienna_boundary_changes.geojson.
"""

import json
import sys
import time
import urllib.request
import urllib.error
import re

import geopandas as gpd
import pandas as pd
import shapely
import numpy as np

# ---------------------------------------------------------------------------
# 1. Overpass query helper
# ---------------------------------------------------------------------------
OVERPASS_URL = "https://overpass-api.de/api/interpreter"


def overpass_query(query, date=None):
    """Run an Overpass query, optionally at a historical date."""
    full_query = query
    if date:
        # prepend the date setting
        full_query = f"[date:\"{date}\"];\n" + query
    data = f"data={urllib.request.quote(full_query)}"
    req = urllib.request.Request(OVERPASS_URL, data.encode("utf-8"))
    req.add_header("Content-Type", "application/x-www-form-urlencoded")
    # Retry with backoff
    for attempt in range(5):
        try:
            with urllib.request.urlopen(req, timeout=120) as resp:
                return json.loads(resp.read().decode("utf-8"))
        except urllib.error.HTTPError as e:
            if e.code == 429:
                time.sleep(2 ** attempt)
                continue
            raise
        except Exception as e:
            if attempt < 4:
                time.sleep(2 ** attempt)
                continue
            raise
    return None


def osm_json_to_gdf(osm_json):
    """Convert Overpass JSON output (nodes, ways, relations) into a GeoDataFrame
    with MultiPolygon geometries for relations tagged as boundary/administrative.
    
    We look for relations with type=boundary or boundary=administrative.
    """
    # Parse elements
    elements = osm_json.get("elements", [])
    
    # Collect nodes by id
    nodes = {}
    ways = {}
    relations = []
    
    for el in elements:
        etype = el.get("type")
        eid = el.get("id")
        if etype == "node":
            nodes[eid] = (el.get("lon"), el.get("lat"))
        elif etype == "way":
            ways[eid] = el.get("nodes", [])
        elif etype == "relation":
            relations.append(el)
    
    # Build linestrings for each way
    def way_to_coords(way_id):
        node_ids = ways.get(way_id, [])
        coords = []
        for nid in node_ids:
            if nid in nodes:
                coords.append(nodes[nid])
        return coords
    
    # For each relation tagged as boundary=administrative (admin_level=9)
    # build the geometry from its members
    features = []
    
    for rel in relations:
        tags = rel.get("tags", {})
        # Check if this is a district boundary relation
        is_boundary = (tags.get("boundary") == "administrative" or 
                       tags.get("type") == "boundary" or
                       tags.get("admin_level") == "9")
        if not is_boundary:
            continue
        
        name = tags.get("name", "")
        # We want district-level boundaries (admin_level=9)
        if tags.get("admin_level") not in (None, "", "9"):
            continue
            
        members = rel.get("members", [])
        
        # Collect outer/inner rings from way members
        outer_rings = []
        inner_rings = []
        
        for mem in members:
            if mem.get("type") != "way":
                continue
            role = mem.get("role", "")
            ref = mem.get("ref")
            if ref not in ways:
                continue
            coords = way_to_coords(ref)
            if len(coords) < 3:
                continue
            if role == "inner":
                inner_rings.append(coords)
            else:  # outer or no role
                outer_rings.append(coords)
        
        # Build polygons
        polygons = []
        for outer in outer_rings:
            if len(outer) < 3:
                continue
            # Close ring if not closed
            if outer[0] != outer[-1]:
                outer = outer + [outer[0]]
            try:
                outer_ls = shapely.LineString(outer)
                if not outer_ls.is_valid or outer_ls.is_empty:
                    continue
                outer_ring = shapely.polygon.orient(shapely.Polygon(outer_ls))
                
                # Find inners that are contained in this outer
                inners_for_this = []
                remaining_inners = []
                for inner in inner_rings:
                    if inner[0] != inner[-1]:
                        inner = inner + [inner[0]]
                    try:
                        inner_ls = shapely.LineString(inner)
                        if inner_ls.is_valid and not inner_ls.is_empty:
                            ipoly = shapely.Polygon(inner_ls)
                            if outer_ring.contains(ipoly.representative_point()) or outer_ring.intersects(ipoly):
                                inners_for_this.append(ipoly)
                            else:
                                remaining_inners.append(inner)
                        else:
                            remaining_inners.append(inner)
                    except Exception:
                        remaining_inners.append(inner)
                
                inner_rings = remaining_inners
                
                if inners_for_this:
                    poly = shapely.Polygon(outer_ls, [list(i.exterior.coords) for i in inners_for_this])
                else:
                    poly = outer_ls
                
                if poly.is_valid and not poly.is_empty:
                    polygons.append(poly)
            except Exception:
                continue
        
        if polygons:
            geom = shapely.unary_union(polygons) if len(polygons) > 1 else polygons[0]
            if geom.geom_type == "Polygon":
                geom = shapely.geometry.MultiPolygon([geom])
            elif geom.geom_type == "MultiPolygon":
                pass
            else:
                continue
            
            features.append({
                "geometry": geom,
                "name": name,
                "osm_id": rel.get("id"),
                "admin_level": tags.get("admin_level", "9"),
            })
    
    if not features:
        return gpd.GeoDataFrame()
    
    gdf = gpd.GeoDataFrame(features, crs="EPSG:4326")
    return gdf


def fetch_vienna_districts(date=None):
    """Fetch all Vienna district boundary relations (admin_level=9) via Overpass."""
    # Use a bounding box around Vienna to limit results, and filter by admin_level=9
    query = """
    [out:json][timeout:120];
    (
      relation["admin_level"="9"]["boundary"="administrative"](48.1,16.1,48.35,16.6);
    );
    out body;
    >;
    out skel qt;
    """
    return overpass_query(query, date=date)


# ---------------------------------------------------------------------------
# 2. Name normalisation
# ---------------------------------------------------------------------------
# Vienna districts and their expected normalised names
VIENNA_DISTRICTS = {
    1: "Innere Stadt",
    2: "Leopoldstadt",
    3: "Landstraße",
    4: "Wieden",
    5: "Margareten",
    6: "Mariahilf",
    7: "Neubau",
    8: "Josefstadt",
    9: "Alsergrund",
    10: "Favoriten",
    11: "Simmering",
    12: "Meidling",
    13: "Hietzing",
    14: "Penzing",
    15: "Rudolfsheim-Fünfhaus",
    16: "Ottakring",
    17: "Hernals",
    18: "Währing",
    19: "Döbling",
    20: "Brigittenau",
    21: "Floridsdorf",
    22: "Donaustadt",
    23: "Liesing",
}


def normalise_name(name):
    """Normalise a district name to its canonical form."""
    if not name or pd.isna(name):
        return None
    # Strip whitespace, handle known variations
    n = name.strip()
    # Some common variations in casing/spelling
    lookup = {
        "innere stadt": "Innere Stadt",
        "1. bezirk": "Innere Stadt",
        "1. bezirk innere stadt": "Innere Stadt",
        "leopoldstadt": "Leopoldstadt",
        "2. bezirk": "Leopoldstadt",
        "landstraße": "Landstraße",
        "landstrasse": "Landstraße",
        "3. bezirk": "Landstraße",
        "wieden": "Wieden",
        "4. bezirk": "Wieden",
        "margareten": "Margareten",
        "5. bezirk": "Margareten",
        "margarethen": "Margareten",
        "mariahilf": "Mariahilf",
        "6. bezirk": "Mariahilf",
        "neubau": "Neubau",
        "7. bezirk": "Neubau",
        "josefstadt": "Josefstadt",
        "8. bezirk": "Josefstadt",
        "alsergrund": "Alsergrund",
        "9. bezirk": "Alsergrund",
        "favoriten": "Favoriten",
        "10. bezirk": "Favoriten",
        "simmering": "Simmering",
        "11. bezirk": "Simmering",
        "meidling": "Meidling",
        "12. bezirk": "Meidling",
        "hietzing": "Hietzing",
        "13. bezirk": "Hietzing",
        "penzing": "Penzing",
        "14. bezirk": "Penzing",
        "rudolfsheim-fünfhaus": "Rudolfsheim-Fünfhaus",
        "rudolfsheim-fuenfhaus": "Rudolfsheim-Fünfhaus",
        "rudolfsheim funfhaus": "Rudolfsheim-Fünfhaus",
        "rudolfsheim fünfhaus": "Rudolfsheim-Fünfhaus",
        "15. bezirk": "Rudolfsheim-Fünfhaus",
        "ottakring": "Ottakring",
        "16. bezirk": "Ottakring",
        "hernals": "Hernals",
        "17. bezirk": "Hernals",
        "währing": "Währing",
        "waehring": "Währing",
        "18. bezirk": "Währing",
        "döbling": "Döbling",
        "doebling": "Döbling",
        "19. bezirk": "Döbling",
        "brigittenau": "Brigittenau",
        "20. bezirk": "Brigittenau",
        "floridsdorf": "Floridsdorf",
        "21. bezirk": "Floridsdorf",
        "donaustadt": "Donaustadt",
        "22. bezirk": "Donaustadt",
        "liesing": "Liesing",
        "23. bezirk": "Liesing",
    }
    key = n.lower().replace("  ", " ").replace("-", "-")
    if key in lookup:
        return lookup[key]
    # Try fuzzy matching
    for pattern, canon in lookup.items():
        if pattern in key or key in pattern:
            return canon
    return n


# ---------------------------------------------------------------------------
# 3. Main pipeline
# ---------------------------------------------------------------------------
def main():
    print("Fetching current Vienna district boundaries...")
    current_osm = fetch_vienna_districts()
    current_gdf = osm_json_to_gdf(current_osm)
    print(f"  Got {len(current_gdf)} current district features")
    
    print("Fetching 2014-01-01 Vienna district boundaries...")
    historical_osm = fetch_vienna_districts(date="2014-01-01T00:00:00Z")
    historical_gdf = osm_json_to_gdf(historical_osm)
    print(f"  Got {len(historical_gdf)} historical district features")
    
    if current_gdf.empty or historical_gdf.empty:
        print("ERROR: Empty data returned. Saving raw Overpass responses for debugging.")
        with open("/work/current_raw.json", "w") as f:
            json.dump(current_osm, f, indent=2)
        with open("/work/historical_raw.json", "w") as f:
            json.dump(historical_osm, f, indent=2)
        print("Raw data saved to /work/current_raw.json and /work/historical_raw.json")
        return
    
    # Normalise names
    current_gdf["district_name"] = current_gdf["name"].apply(normalise_name)
    historical_gdf["district_name"] = historical_gdf["name"].apply(normalise_name)
    
    print("\nCurrent district names:")
    for n in sorted(current_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    print("\nHistorical district names:")
    for n in sorted(historical_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    
    # Remove non-district features (like "Amtshaus für den 16. Bezirk")
    valid_names = set(VIENNA_DISTRICTS.values())
    current_gdf = current_gdf[current_gdf["district_name"].isin(valid_names)].copy()
    historical_gdf = historical_gdf[historical_gdf["district_name"].isin(valid_names)].copy()
    
    print(f"\nAfter filtering: {len(current_gdf)} current, {len(historical_gdf)} historical districts")
    
    # Check for matching districts
    current_names = set(current_gdf["district_name"].dropna())
    historical_names = set(historical_gdf["district_name"].dropna())
    
    common = current_names & historical_names
    only_current = current_names - historical_names
    only_historical = historical_names - current_names
    
    print(f"Common districts: {len(common)}")
    if only_current:
        print(f"Only in current: {only_current}")
    if only_historical:
        print(f"Only in historical: {only_historical}")
    
    # Build a unified dissolve of all current and all historical boundaries
    current_union = current_gdf.geometry.unary_union
    historical_union = historical_gdf.geometry.unary_union
    
    print(f"\nCurrent total boundary area: {current_union.area:.4f} sq deg")
    print(f"Historical total boundary area: {historical_union.area:.4f} sq deg")
    
    # Symmetric difference gives changed areas
    added = shapely.difference(current_union, historical_union)
    removed = shapely.difference(historical_union, current_union)
    unchanged = shapely.intersection(current_union, historical_union)
    
    print(f"Added area: {added.area:.6f}")
    print(f"Removed area: {removed.area:.6f}")
    print(f"Unchanged area: {unchanged.area:.6f}")
    
    # Prepare fragments
    fragments = []
    
    # Added fragments
    if not added.is_empty:
        if added.geom_type == "MultiPolygon":
            for poly in added.geoms:
                fragments.append((shapely.geometry.MultiPolygon([poly]), "added_since_2014"))
        elif added.geom_type == "Polygon":
            fragments.append((shapely.geometry.MultiPolygon([added]), "added_since_2014"))
    
    # Removed fragments
    if not removed.is_empty:
        if removed.geom_type == "MultiPolygon":
            for poly in removed.geoms:
                fragments.append((shapely.geometry.MultiPolygon([poly]), "removed_since_2014"))
        elif removed.geom_type == "Polygon":
            fragments.append((shapely.geometry.MultiPolygon([removed]), "removed_since_2014"))
    
    # Unchanged fragments - we can split by district for better attribution
    unchanged_frags = []
    for _, row in current_gdf.iterrows():
        dist_name = row["district_name"]
        dist_geom = row["geometry"]
        unchanged_part = shapely.intersection(dist_geom, unchanged)
        if not unchanged_part.is_empty:
            if unchanged_part.geom_type == "MultiPolygon":
                for poly in unchanged_part.geoms:
                    unchanged_frags.append((shapely.geometry.MultiPolygon([poly]), "unchanged", dist_name))
            elif unchanged_part.geom_type == "Polygon":
                unchanged_frags.append((shapely.geometry.MultiPolygon([unchanged_part]), "unchanged", dist_name))
            else:
                unchanged_frags.append((shapely.geometry.MultiPolygon([unchanged_part]), "unchanged", dist_name))
    
    # For added/removed, determine which district they belong to
    # We'll assign the district_name based on which current/historical district contains/intersects them
    added_with_dist = []
    for geom, ctype in fragments:
        if ctype == "added_since_2014":
            # Find which current district(s) this belongs to
            touching_districts = set()
            for _, row in current_gdf.iterrows():
                if geom.intersects(row["geometry"]):
                    touching_districts.add(row["district_name"])
            dist_name = ", ".join(sorted(touching_districts)) if touching_districts else "unknown"
            added_with_dist.append((geom, ctype, dist_name))
        elif ctype == "removed_since_2014":
            touching_districts = set()
            for _, row in historical_gdf.iterrows():
                if geom.intersects(row["geometry"]):
                    touching_districts.add(row["district_name"])
            dist_name = ", ".join(sorted(touching_districts)) if touching_districts else "unknown"
            added_with_dist.append((geom, ctype, dist_name))
    
    # Combine all fragments
    all_fragments = []
    for geom, ctype, dname in added_with_dist:
        all_fragments.append({"geometry": geom, "change_type": ctype, "district_name": dname})
    for geom, ctype, dname in unchanged_frags:
        all_fragments.append({"geometry": geom, "change_type": ctype, "district_name": dname})
    
    # Build GeoDataFrame
    result_gdf = gpd.GeoDataFrame(all_fragments, crs="EPSG:4326")
    
    if result_gdf.empty:
        print("No fragments produced — all boundaries unchanged.")
        # Create a dummy feature showing all current boundaries as unchanged
        for _, row in current_gdf.iterrows():
            result_gdf = gpd.GeoDataFrame([
                {"geometry": row["geometry"], "change_type": "unchanged", "district_name": row["district_name"]}
            ], crs="EPSG:4326")
    else:
        print(f"\nFragments created: {len(result_gdf)}")
        print(result_gdf["change_type"].value_counts())
    
    # -----------------------------------------------------------------------
    # 4. Compute touches_changed for each district
    # -----------------------------------------------------------------------
    # For each district, check if its geometry borders (touches) any added or
    # removed area.
    
    # Get all added/removed geometries combined
    changed_areas = []
    for _, row in result_gdf.iterrows():
        if row["change_type"] in ("added_since_2014", "removed_since_2014"):
            changed_areas.append(row["geometry"])
    
    if changed_areas:
        changed_union = shapely.unary_union(changed_areas)
    else:
        changed_union = shapely.geometry.MultiPolygon()
    
    # For each district in current data, check if it touches changed areas
    district_touches = {}
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        geom = row["geometry"]
        if not changed_union.is_empty and geom.touches(changed_union):
            district_touches[dname] = True
        else:
            # Also check if the district itself IS an added/removed part
            # i.e. if the district geometry overlaps with changed areas
            if not changed_union.is_empty and geom.intersects(changed_union):
                # Check if it's more than just touching - if it contains changed area
                overlap = shapely.intersection(geom, changed_union)
                if not overlap.is_empty and overlap.area > 0:
                    district_touches[dname] = True
                else:
                    district_touches[dname] = False
            else:
                district_touches[dname] = False
    
    # Assign touches_changed to each fragment
    result_gdf["touches_changed"] = result_gdf["district_name"].apply(
        lambda dn: any(district_touches.get(d.strip(), False) for d in dn.split(", "))
        if dn != "unknown" else False
    )
    
    # Simplify geometries a bit to keep file size manageable
    result_gdf["geometry"] = result_gdf["geometry"].simplify(0.00001, preserve_topology=True)
    
    # Output
    output_path = "/work/vienna_boundary_changes.geojson"
    result_gdf.to_file(output_path, driver="GeoJSON")
    print(f"\nOutput written to {output_path}")
    print(f"Features: {len(result_gdf)}")
    print(f"Columns: {result_gdf.columns.tolist()}")
    print(f"\nSample:\n{result_gdf[['change_type', 'district_name', 'touches_changed']].head(10)}")
    print(f"\nDistrict touch summary:")
    for dname, touches in sorted(district_touches.items()):
        print(f"  {dname}: touches_changed={touches}")


if __name__ == "__main__":
    main()

"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass via osmnx,
compare to find changes, and output vienna_boundary_changes.geojson.

Each feature is a MultiPolygon geometry fragment with:
  - change_type: added_since_2014 | removed_since_2014 | unchanged
  - district_name: normalised district name(s) the fragment belongs to
  - touches_changed: bool, whether that district's boundary touches any changed area
"""

import os
import sys
import json
from collections import OrderedDict

import osmnx as ox
import osmnx._overpass as op
import geopandas as gpd
import pandas as pd
import shapely
import shapely.geometry
import numpy as np

ox.settings.log_console = False

# ---------------------------------------------------------------------------
# 1. Overpass query via osmnx internals
# ---------------------------------------------------------------------------

def fetch_vienna_districts_osmnx(date=None):
    """Fetch Vienna district (admin_level=9) boundary relations via Overpass.
    
    Uses osmnx's _overpass_request to handle HTTP correctly.
    Returns a GeoDataFrame with columns: geometry, name, district_name, osm_id.
    """
    # Bounding box covering Vienna
    bbox = (48.1, 16.1, 48.35, 16.6)  # south, west, north, east
    
    if date:
        query = f"""
[out:json][timeout:180][date:"{date}"];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    else:
        query = f"""
[out:json][timeout:180];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    
    data = OrderedDict()
    data['data'] = query
    result = op._overpass_request(data)
    
    elements = result.get("elements", [])
    
    # Parse: collect nodes, ways, relations
    nodes = {}
    ways = {}
    relations = []
    
    for el in elements:
        etype = el.get("type")
        eid = el.get("id")
        if etype == "node":
            nodes[eid] = (el.get("lon"), el.get("lat"))
        elif etype == "way":
            ways[eid] = el.get("nodes", [])
        elif etype == "relation":
            relations.append(el)
    
    def way_to_coords(way_id):
        node_ids = ways.get(way_id, [])
        return [nodes[nid] for nid in node_ids if nid in nodes]
    
    features = []
    
    for rel in relations:
        tags = rel.get("tags", {})
        name = tags.get("name", "")
        admin_level = tags.get("admin_level", "")
        
        # Only Vienna districts (admin_level=9)
        if admin_level != "9":
            continue
        
        members = rel.get("members", [])
        
        outer_rings = []
        inner_rings = []
        
        for mem in members:
            if mem.get("type") != "way":
                continue
            role = mem.get("role", "")
            ref = mem.get("ref")
            if ref not in ways:
                continue
            coords = way_to_coords(ref)
            if len(coords) < 3:
                continue
            if role == "inner":
                inner_rings.append(coords)
            else:
                outer_rings.append(coords)
        
        polygons = []
        for outer in outer_rings:
            if len(outer) < 3:
                continue
            if outer[0] != outer[-1]:
                outer = outer + [outer[0]]
            try:
                outer_ls = shapely.LineString(outer)
                if not outer_ls.is_valid or outer_ls.is_empty:
                    continue
                
                inners_for_this = []
                remaining = []
                for inner in inner_rings:
                    if inner[0] != inner[-1]:
                        inner = inner + [inner[0]]
                    try:
                        inner_ls = shapely.LineString(inner)
                        if inner_ls.is_valid and not inner_ls.is_empty:
                            ipoly = shapely.Polygon(inner_ls)
                            outer_poly = shapely.Polygon(outer_ls)
                            if outer_poly.contains(ipoly.representative_point()) or outer_poly.intersects(ipoly):
                                inners_for_this.append(ipoly)
                            else:
                                remaining.append(inner)
                        else:
                            remaining.append(inner)
                    except Exception:
                        remaining.append(inner)
                
                inner_rings = remaining
                
                if inners_for_this:
                    poly = shapely.Polygon(outer_ls, [list(i.exterior.coords) for i in inners_for_this])
                else:
                    poly = shapely.Polygon(outer_ls)
                
                # Fix orientation and validate
                poly = shapely.polygon.orient(poly)
                if poly.is_valid and not poly.is_empty:
                    polygons.append(poly)
            except Exception:
                continue
        
        if polygons:
            geom = shapely.unary_union(polygons)
            if geom.geom_type == "Polygon":
                geom = shapely.geometry.MultiPolygon([geom])
            
            if not geom.is_empty:
                features.append({
                    "geometry": geom,
                    "name": name,
                    "osm_id": rel.get("id"),
                })
    
    if not features:
        return gpd.GeoDataFrame()
    
    gdf = gpd.GeoDataFrame(features, crs="EPSG:4326")
    return gdf


# ---------------------------------------------------------------------------
# 2. Name normalisation
# ---------------------------------------------------------------------------

# Canonical Vienna district names
VIENNA_DISTRICTS_CANON = {
    "innere stadt": "Innere Stadt",
    "leopoldstadt": "Leopoldstadt",
    "landstraße": "Landstraße",
    "landstrasse": "Landstraße",
    "wieden": "Wieden",
    "margareten": "Margareten",
    "margarethen": "Margareten",
    "mariahilf": "Mariahilf",
    "neubau": "Neubau",
    "josefstadt": "Josefstadt",
    "alsergrund": "Alsergrund",
    "favoriten": "Favoriten",
    "simmering": "Simmering",
    "meidling": "Meidling",
    "hietzing": "Hietzing",
    "penzing": "Penzing",
    "rudolfsheim-fünfhaus": "Rudolfsheim-Fünfhaus",
    "rudolfsheim-fuenfhaus": "Rudolfsheim-Fünfhaus",
    "ottakring": "Ottakring",
    "hernals": "Hernals",
    "währing": "Währing",
    "waehring": "Währing",
    "döbling": "Döbling",
    "doebling": "Döbling",
    "brigittenau": "Brigittenau",
    "floridsdorf": "Floridsdorf",
    "donaustadt": "Donaustadt",
    "liesing": "Liesing",
}


def normalise_name(name):
    """Normalise a district name to its canonical form."""
    if not name or pd.isna(name):
        return None
    n = name.strip()
    key = n.lower().replace("  ", " ").replace("ß", "ss").replace("-", "-")
    # Direct lookup
    if key in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key]
    # try with ß vs ss variations
    key2 = key.replace("ss", "ß")
    if key2 in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key2]
    # Try fuzzy: if any lookup key is contained in our key
    for pattern, canon in VIENNA_DISTRICTS_CANON.items():
        if pattern in key or key in pattern:
            return canon
    return n


# ---------------------------------------------------------------------------
# 3. Main pipeline
# ---------------------------------------------------------------------------

def main():
    print("=" * 60)
    print("Fetching CURRENT Vienna district boundaries...")
    print("=" * 60)
    current_gdf = fetch_vienna_districts_osmnx()
    print(f"  Got {len(current_gdf)} current district features")
    
    print("\n" + "=" * 60)
    print("Fetching 2014-01-01 Vienna district boundaries...")
    print("=" * 60)
    historical_gdf = fetch_vienna_districts_osmnx(date="2014-01-01T00:00:00Z")
    print(f"  Got {len(historical_gdf)} historical district features")
    
    if current_gdf.empty or historical_gdf.empty:
        print("ERROR: Empty data returned!")
        sys.exit(1)
    
    # Normalise names
    current_gdf["district_name"] = current_gdf["name"].apply(normalise_name)
    historical_gdf["district_name"] = historical_gdf["name"].apply(normalise_name)
    
    print("\nCurrent district names:")
    for n in sorted(current_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    print("\nHistorical district names:")
    for n in sorted(historical_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    
    # Remove non-Vienna-district features (e.g. Klein-Engersdorf, Enzersfeld in historical data)
    valid_names = set(VIENNA_DISTRICTS_CANON.values())
    current_gdf = current_gdf[current_gdf["district_name"].isin(valid_names)].copy()
    historical_gdf = historical_gdf[historical_gdf["district_name"].isin(valid_names)].copy()
    
    print(f"\nAfter filtering to Vienna districts only:")
    print(f"  Current: {len(current_gdf)} districts")
    print(f"  Historical: {len(historical_gdf)} districts")
    
    # Check for missing districts
    current_names = set(current_gdf["district_name"].dropna())
    historical_names = set(historical_gdf["district_name"].dropna())
    
    common = current_names & historical_names
    only_current = current_names - historical_names
    only_historical = historical_names - current_names
    
    print(f"\nCommon districts: {len(common)}")
    if only_current:
        print(f"Only in current: {only_current}")
    if only_historical:
        print(f"Only in historical: {only_historical}")
    
    # Build unified dissolve
    current_union = current_gdf.geometry.unary_union
    historical_union = historical_gdf.geometry.unary_union
    
    print(f"\nCurrent total area: {current_union.area:.6f} sq deg")
    print(f"Historical total area: {historical_union.area:.6f} sq deg")
    
    # Symmetric difference
    added = shapely.difference(current_union, historical_union)
    removed = shapely.difference(historical_union, current_union)
    unchanged = shapely.intersection(current_union, historical_union)
    
    print(f"Added area: {added.area:.6f}" if not added.is_empty else "Added: empty")
    print(f"Removed area: {removed.area:.6f}" if not removed.is_empty else "Removed: empty")
    print(f"Unchanged area: {unchanged.area:.6f}" if not unchanged.is_empty else "Unchanged: empty")
    
    # Collect fragments
    fragments = []
    
    # Helper to split MultiPolygon into individual fragments
    def add_fragments(geom, ctype, dist_name_fn):
        """Add geometry fragments, assigning district name via dist_name_fn."""
        if geom.is_empty:
            return
        if geom.geom_type == "MultiPolygon":
            polys = list(geom.geoms)
        elif geom.geom_type == "Polygon":
            polys = [geom]
        elif geom.geom_type == "GeometryCollection":
            polys = [g for g in geom.geoms if g.geom_type in ("Polygon", "MultiPolygon")]
        else:
            polys = [geom]
        
        for poly in polys:
            if poly.is_empty:
                continue
            if poly.geom_type == "Polygon":
                mp = shapely.geometry.MultiPolygon([poly])
            elif poly.geom_type == "MultiPolygon":
                mp = poly
            else:
                continue
            dname = dist_name_fn(mp)
            fragments.append({
                "geometry": mp,
                "change_type": ctype,
                "district_name": dname,
            })
    
    # Added fragments - assign to current districts that contain/intersect them
    def added_dist_name(geom):
        districts = set()
        for _, row in current_gdf.iterrows():
            if geom.intersects(row["geometry"]):
                districts.add(row["district_name"])
        return ", ".join(sorted(districts)) if districts else "unknown"
    
    add_fragments(added, "added_since_2014", added_dist_name)
    
    # Removed fragments - assign to historical districts
    def removed_dist_name(geom):
        districts = set()
        for _, row in historical_gdf.iterrows():
            if geom.intersects(row["geometry"]):
                districts.add(row["district_name"])
        return ", ".join(sorted(districts)) if districts else "unknown"
    
    add_fragments(removed, "removed_since_2014", removed_dist_name)
    
    # Unchanged fragments - assign per current district
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        dist_geom = row["geometry"]
        unchanged_part = shapely.intersection(dist_geom, unchanged)
        if not unchanged_part.is_empty:
            def dn_fn(geom, dn=dname):
                return dn
            add_fragments(unchanged_part, "unchanged", dn_fn)
    
    result_gdf = gpd.GeoDataFrame(fragments, crs="EPSG:4326")
    
    if result_gdf.empty:
        print("No fragments — all boundaries unchanged. Using current boundaries as unchanged.")
        result_gdf = current_gdf.copy()
        result_gdf["change_type"] = "unchanged"
        result_gdf["district_name"] = result_gdf["district_name"]
        result_gdf = result_gdf[["geometry", "change_type", "district_name"]]
    
    print(f"\nFragments created: {len(result_gdf)}")
    print(result_gdf["change_type"].value_counts())
    
    # -----------------------------------------------------------------------
    # 4. Compute touches_changed for each district
    # -----------------------------------------------------------------------
    
    # Collect all changed geometries
    changed_geoms = []
    for _, row in result_gdf.iterrows():
        if row["change_type"] in ("added_since_2014", "removed_since_2014"):
            changed_geoms.append(row["geometry"])
    
    if changed_geoms:
        changed_union = shapely.unary_union(changed_geoms)
    else:
        changed_union = shapely.geometry.MultiPolygon()
    
    # For each district in current data, check if it touches changed areas
    # We use a small buffer (1e-7 degrees ~ 1cm) to capture near-touching
    district_touches = {}
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        geom = row["geometry"]
        if not changed_union.is_empty:
            # Use touches with a tiny buffer for robustness
            buffered = changed_union.buffer(1e-8)
            if geom.touches(changed_union) or geom.intersects(buffered):
                district_touches[dname] = True
            else:
                district_touches[dname] = False
        else:
            district_touches[dname] = False
    
    # For historical-only districts (if any), also check
    for _, row in historical_gdf.iterrows():
        dname = row["district_name"]
        if dname not in district_touches:
            geom = row["geometry"]
            if not changed_union.is_empty and geom.intersects(changed_union):
                district_touches[dname] = True
            else:
                district_touches[dname] = False
    
    def get_touches(dn_str):
        """Check if any district in a comma-separated string touches changed areas."""
        if pd.isna(dn_str) or dn_str == "unknown":
            return False
        for d in dn_str.split(","):
            d = d.strip()
            if district_touches.get(d, False):
                return True
        return False
    
    result_gdf["touches_changed"] = result_gdf["district_name"].apply(get_touches)
    
    # Simplify geometries slightly to keep file manageable
    result_gdf["geometry"] = result_gdf["geometry"].simplify(0.00001, preserve_topology=True)
    
    # Repair any invalid geometries from simplification
    result_gdf["geometry"] = result_gdf["geometry"].buffer(0)
    
    # Ensure MultiPolygon type
    def ensure_multipolygon(geom):
        if geom.geom_type == "Polygon":
            return shapely.geometry.MultiPolygon([geom])
        return geom
    
    result_gdf["geometry"] = result_gdf["geometry"].apply(ensure_multipolygon)
    
    # -----------------------------------------------------------------------
    # 5. Output
    # -----------------------------------------------------------------------
    output_path = "/work/vienna_boundary_changes.geojson"
    result_gdf.to_file(output_path, driver="GeoJSON")
    
    print(f"\n{'=' * 60}")
    print(f"Output written to {output_path}")
    print(f"Features: {len(result_gdf)}")
    print(f"Columns: {result_gdf.columns.tolist()}")
    print(f"\nSample features:")
    print(result_gdf[["change_type", "district_name", "touches_changed"]].head(15))
    
    print(f"\nDistrict touch summary:")
    for dname, touches in sorted(district_touches.items()):
        print(f"  {dname}: touches_changed={touches}")
    
    # Verify output
    verify = gpd.read_file(output_path)
    print(f"\nVerification: read back {len(verify)} features")
    print(verify[["change_type", "district_name", "touches_changed"]].head(10))
    print(f"Geometry types: {verify.geometry.geom_type.unique()}")


if __name__ == "__main__":
    main()

"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass via osmnx,
compare to find changes, and output vienna_boundary_changes.geojson.

Each feature is a MultiPolygon geometry fragment with:
  - change_type: added_since_2014 | removed_since_2014 | unchanged
  - district_name: normalised district name(s) the fragment belongs to
  - touches_changed: bool, whether that district's boundary touches any changed area
"""

import os
import sys
import json
from collections import OrderedDict

import osmnx as ox
import osmnx._overpass as op
import geopandas as gpd
import pandas as pd
import shapely
import shapely.geometry
import numpy as np

ox.settings.log_console = False

# ---------------------------------------------------------------------------
# 1. Overpass query via osmnx internals
# ---------------------------------------------------------------------------

def assemble_ring_from_ways(way_refs, role, ways_dict, nodes_dict):
    """Assemble a closed ring from a list of way references.
    
    Ways in OSM relations are directed segments that need to be chained
    end-to-end by matching their start/end node IDs.
    
    Returns list of (lon, lat) tuples forming a closed ring, or None.
    """
    if not way_refs:
        return None
    
    # Get all segments with their start/end node ids
    segments = []
    for ref in way_refs:
        if ref not in ways_dict:
            continue
        node_ids = ways_dict[ref]
        if len(node_ids) < 2:
            continue
        coords = [nodes_dict.get(nid) for nid in node_ids]
        coords = [c for c in coords if c is not None]
        if len(coords) < 2:
            continue
        segments.append({
            'ref': ref,
            'coords': coords,
            'start_node': node_ids[0],
            'end_node': node_ids[-1],
        })
    
    if not segments:
        return None
    
    # Chain segments together
    used = [False] * len(segments)
    
    # Start with the first segment
    ring_coords = list(segments[0]['coords'])
    current_end_node = segments[0]['end_node']
    used[0] = True
    
    # Keep connecting segments until we close the ring or can't find more
    changed = True
    while changed:
        changed = False
        for i, seg in enumerate(segments):
            if used[i]:
                continue
            if seg['start_node'] == current_end_node:
                # Append coords (skip first to avoid duplication)
                ring_coords.extend(seg['coords'][1:])
                current_end_node = seg['end_node']
                used[i] = True
                changed = True
                break
            elif seg['end_node'] == current_end_node:
                # Reverse segment
                rev_coords = list(reversed(seg['coords']))
                ring_coords.extend(rev_coords[1:])
                current_end_node = seg['start_node']
                used[i] = True
                changed = True
                break
    
    # Check if the ring is closed
    if ring_coords[0] != ring_coords[-1]:
        ring_coords.append(ring_coords[0])
    
    return ring_coords


def fetch_vienna_districts_osmnx(date=None):
    """Fetch Vienna district (admin_level=9) boundary relations via Overpass.
    
    Returns a GeoDataFrame with columns: geometry, name, district_name, osm_id.
    """
    bbox = (48.1, 16.1, 48.35, 16.6)  # south, west, north, east
    
    if date:
        query = f"""
[out:json][timeout:180][date:"{date}"];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    else:
        query = f"""
[out:json][timeout:180];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    
    data = OrderedDict()
    data['data'] = query
    result = op._overpass_request(data)
    
    elements = result.get("elements", [])
    
    nodes = {}
    ways = {}
    relations = []
    
    for el in elements:
        etype = el.get("type")
        eid = el.get("id")
        if etype == "node":
            nodes[eid] = (el.get("lon"), el.get("lat"))
        elif etype == "way":
            ways[eid] = el.get("nodes", [])
        elif etype == "relation":
            relations.append(el)
    
    features = []
    
    for rel in relations:
        tags = rel.get("tags", {})
        name = tags.get("name", "")
        admin_level = tags.get("admin_level", "")
        
        if admin_level != "9":
            continue
        
        members = rel.get("members", [])
        
        # Collect way refs by role
        outer_refs = []
        inner_refs = []
        
        for mem in members:
            if mem.get("type") != "way":
                continue
            role = mem.get("role", "")
            ref = mem.get("ref")
            if ref not in ways:
                continue
            if role == "inner":
                inner_refs.append(ref)
            else:
                outer_refs.append(ref)
        
        # Assemble each outer ring
        outer_polygons = []
        
        # Group outer ways into separate rings
        # Multiple outer rings mean multipolygon
        unused_outer = list(outer_refs)
        
        while unused_outer:
            # Try to build a ring starting from the first unused way
            ring = assemble_ring_from_ways(unused_outer, 'outer', ways, nodes)
            if ring is None or len(ring) < 4:
                break
            
            # Remove used ways
            # (We need to track which refs were used — for simplicity, remove those
            # that start/end match what we consumed)
            # Better approach: find which refs have coords in the ring
            
            # Create linestring and polygon
            try:
                ls = shapely.LineString(ring)
                if not ls.is_valid or ls.is_empty or ls.length < 1e-10:
                    break
                poly = shapely.Polygon(ls)
                poly = shapely.polygon.orient(poly)
                if poly.is_valid and not poly.is_empty and poly.area > 1e-12:
                    outer_polygons.append(poly)
            except Exception:
                break
            
            # Remove used ways from unused_outer
            # We find which refs contributed by checking coordinates
            ring_coords_set = set(tuple(round(c, 10) for c in coord) for coord in ring)
            new_unused = []
            for ref in unused_outer:
                if ref not in ways:
                    continue
                node_ids = ways[ref]
                # Check if this way contributed to the ring
                contributed = False
                for nid in node_ids:
                    coord = nodes.get(nid)
                    if coord and tuple(round(c, 10) for c in coord) in ring_coords_set:
                        contributed = True
                        break
                if not contributed:
                    new_unused.append(ref)
            
            if len(new_unused) == len(unused_outer):
                # No progress - break to avoid infinite loop
                break
            unused_outer = new_unused
        
        if not outer_polygons:
            continue
        
        # Now handle inner rings
        # For each outer polygon, find contained inner rings
        inner_polygons = []
        for ref in inner_refs:
            if ref not in ways:
                continue
            node_ids = ways[ref]
            coords = [nodes.get(nid) for nid in node_ids]
            coords = [c for c in coords if c is not None]
            if len(coords) < 3:
                continue
            if coords[0] != coords[-1]:
                coords.append(coords[0])
            try:
                inner_ls = shapely.LineString(coords)
                if inner_ls.is_valid and not inner_ls.is_empty and inner_ls.length > 1e-10:
                    inner_poly = shapely.Polygon(inner_ls)
                    inner_poly = shapely.polygon.orient(inner_poly)
                    if inner_poly.is_valid and not inner_poly.is_empty:
                        inner_polygons.append(inner_poly)
            except Exception:
                continue
        
        # Build final polygons with holes
        final_polygons = []
        for outer_poly in outer_polygons:
            # Find inners contained in this outer
            inners_for_this = []
            remaining_inners = []
            for ipoly in inner_polygons:
                try:
                    if outer_poly.contains(ipoly.representative_point()) or outer_poly.intersects(ipoly):
                        inners_for_this.append(ipoly)
                    else:
                        remaining_inners.append(ipoly)
                except Exception:
                    remaining_inners.append(ipoly)
            inner_polygons = remaining_inners
            
            if inners_for_this:
                # Create polygon with holes
                try:
                    outer_coords = list(outer_poly.exterior.coords)
                    inner_coords_list = [list(ip.exterior.coords) for ip in inners_for_this]
                    poly = shapely.Polygon(outer_coords, inner_coords_list)
                    poly = shapely.polygon.orient(poly)
                    if poly.is_valid and not poly.is_empty:
                        final_polygons.append(poly)
                except Exception:
                    final_polygons.append(outer_poly)
            else:
                final_polygons.append(outer_poly)
        
        if final_polygons:
            geom = shapely.unary_union(final_polygons)
            if geom.geom_type == "Polygon":
                geom = shapely.geometry.MultiPolygon([geom])
            
            if not geom.is_empty:
                features.append({
                    "geometry": geom,
                    "name": name,
                    "osm_id": rel.get("id"),
                })
    
    if not features:
        return gpd.GeoDataFrame()
    
    gdf = gpd.GeoDataFrame(features, crs="EPSG:4326")
    return gdf


# ---------------------------------------------------------------------------
# 2. Name normalisation
# ---------------------------------------------------------------------------

# Canonical Vienna district names
VIENNA_DISTRICTS_CANON = {
    "innere stadt": "Innere Stadt",
    "leopoldstadt": "Leopoldstadt",
    "landstraße": "Landstraße",
    "landstrasse": "Landstraße",
    "wieden": "Wieden",
    "margareten": "Margareten",
    "margarethen": "Margareten",
    "mariahilf": "Mariahilf",
    "neubau": "Neubau",
    "josefstadt": "Josefstadt",
    "alsergrund": "Alsergrund",
    "favoriten": "Favoriten",
    "simmering": "Simmering",
    "meidling": "Meidling",
    "hietzing": "Hietzing",
    "penzing": "Penzing",
    "rudolfsheim-fünfhaus": "Rudolfsheim-Fünfhaus",
    "rudolfsheim-fuenfhaus": "Rudolfsheim-Fünfhaus",
    "ottakring": "Ottakring",
    "hernals": "Hernals",
    "währing": "Währing",
    "waehring": "Währing",
    "döbling": "Döbling",
    "doebling": "Döbling",
    "brigittenau": "Brigittenau",
    "floridsdorf": "Floridsdorf",
    "donaustadt": "Donaustadt",
    "liesing": "Liesing",
}


def normalise_name(name):
    """Normalise a district name to its canonical form."""
    if not name or pd.isna(name):
        return None
    n = name.strip()
    key = n.lower().replace("  ", " ").replace("ß", "ss").replace("-", "-")
    if key in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key]
    key2 = key.replace("ss", "ß")
    if key2 in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key2]
    for pattern, canon in VIENNA_DISTRICTS_CANON.items():
        if pattern in key or key in pattern:
            return canon
    return n


# ---------------------------------------------------------------------------
# 3. Main pipeline
# ---------------------------------------------------------------------------

def main():
    print("=" * 60)
    print("Fetching CURRENT Vienna district boundaries...")
    print("=" * 60)
    current_gdf = fetch_vienna_districts_osmnx()
    print(f"  Got {len(current_gdf)} current district features")
    
    print("\n" + "=" * 60)
    print("Fetching 2014-01-01 Vienna district boundaries...")
    print("=" * 60)
    historical_gdf = fetch_vienna_districts_osmnx(date="2014-01-01T00:00:00Z")
    print(f"  Got {len(historical_gdf)} historical district features")
    
    if current_gdf.empty or historical_gdf.empty:
        print("ERROR: Empty data returned!")
        sys.exit(1)
    
    # Normalise names
    current_gdf["district_name"] = current_gdf["name"].apply(normalise_name)
    historical_gdf["district_name"] = historical_gdf["name"].apply(normalise_name)
    
    print("\nCurrent district names:")
    for n in sorted(current_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    print("\nHistorical district names:")
    for n in sorted(historical_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    
    # Remove non-Vienna-district features
    valid_names = set(VIENNA_DISTRICTS_CANON.values())
    current_gdf = current_gdf[current_gdf["district_name"].isin(valid_names)].copy()
    historical_gdf = historical_gdf[historical_gdf["district_name"].isin(valid_names)].copy()
    
    print(f"\nAfter filtering to Vienna districts only:")
    print(f"  Current: {len(current_gdf)} districts")
    print(f"  Historical: {len(historical_gdf)} districts")
    
    current_names = set(current_gdf["district_name"].dropna())
    historical_names = set(historical_gdf["district_name"].dropna())
    
    common = current_names & historical_names
    only_current = current_names - historical_names
    only_historical = historical_names - current_names
    
    print(f"\nCommon districts: {len(common)}")
    if only_current:
        print(f"Only in current: {only_current}")
    if only_historical:
        print(f"Only in historical: {only_historical}")
    
    # Build unified dissolve
    current_union = current_gdf.geometry.unary_union
    historical_union = historical_gdf.geometry.unary_union
    
    print(f"\nCurrent total area: {current_union.area:.6f} sq deg")
    print(f"Historical total area: {historical_union.area:.6f} sq deg")
    
    # Symmetric difference
    added = shapely.difference(current_union, historical_union)
    removed = shapely.difference(historical_union, current_union)
    unchanged = shapely.intersection(current_union, historical_union)
    
    print(f"Added area: {added.area:.6f}" if not added.is_empty else "Added: empty")
    print(f"Removed area: {removed.area:.6f}" if not removed.is_empty else "Removed: empty")
    print(f"Unchanged area: {unchanged.area:.6f}" if not unchanged.is_empty else "Unchanged: empty")
    
    # Collect fragments
    fragments = []
    
    def add_fragments(geom, ctype, dist_name_fn):
        if geom.is_empty:
            return
        if geom.geom_type == "MultiPolygon":
            polys = list(geom.geoms)
        elif geom.geom_type == "Polygon":
            polys = [geom]
        elif geom.geom_type == "GeometryCollection":
            polys = [g for g in geom.geoms if g.geom_type in ("Polygon", "MultiPolygon")]
        else:
            polys = [geom]
        
        for poly in polys:
            if poly.is_empty:
                continue
            if poly.geom_type == "Polygon":
                mp = shapely.geometry.MultiPolygon([poly])
            elif poly.geom_type == "MultiPolygon":
                mp = poly
            else:
                continue
            dname = dist_name_fn(mp)
            fragments.append({
                "geometry": mp,
                "change_type": ctype,
                "district_name": dname,
            })
    
    def added_dist_name(geom):
        districts = set()
        for _, row in current_gdf.iterrows():
            if geom.intersects(row["geometry"]):
                districts.add(row["district_name"])
        return ", ".join(sorted(districts)) if districts else "unknown"
    
    add_fragments(added, "added_since_2014", added_dist_name)
    
    def removed_dist_name(geom):
        districts = set()
        for _, row in historical_gdf.iterrows():
            if geom.intersects(row["geometry"]):
                districts.add(row["district_name"])
        return ", ".join(sorted(districts)) if districts else "unknown"
    
    add_fragments(removed, "removed_since_2014", removed_dist_name)
    
    # Unchanged fragments - per current district
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        dist_geom = row["geometry"]
        unchanged_part = shapely.intersection(dist_geom, unchanged)
        if not unchanged_part.is_empty:
            def dn_fn(geom, dn=dname):
                return dn
            add_fragments(unchanged_part, "unchanged", dn_fn)
    
    result_gdf = gpd.GeoDataFrame(fragments, crs="EPSG:4326")
    
    if result_gdf.empty:
        print("No fragments — all boundaries unchanged.")
        result_gdf = current_gdf.copy()
        result_gdf["change_type"] = "unchanged"
        result_gdf["district_name"] = result_gdf["district_name"]
        result_gdf = result_gdf[["geometry", "change_type", "district_name"]]
    
    print(f"\nFragments created: {len(result_gdf)}")
    if not result_gdf.empty:
        print(result_gdf["change_type"].value_counts())
    
    # -----------------------------------------------------------------------
    # 4. Compute touches_changed for each district
    # -----------------------------------------------------------------------
    
    changed_geoms = []
    for _, row in result_gdf.iterrows():
        if row["change_type"] in ("added_since_2014", "removed_since_2014"):
            changed_geoms.append(row["geometry"])
    
    if changed_geoms:
        changed_union = shapely.unary_union(changed_geoms)
    else:
        changed_union = shapely.geometry.MultiPolygon()
    
    district_touches = {}
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        geom = row["geometry"]
        if not changed_union.is_empty:
            # Use touches or a small buffer for robustness
            touches = geom.touches(changed_union)
            if not touches:
                # Also check if the district itself contains changed areas
                overlap = shapely.intersection(geom, changed_union)
                touches = not overlap.is_empty
            district_touches[dname] = touches
        else:
            district_touches[dname] = False
    
    for _, row in historical_gdf.iterrows():
        dname = row["district_name"]
        if dname not in district_touches:
            geom = row["geometry"]
            if not changed_union.is_empty:
                touches = geom.touches(changed_union)
                if not touches:
                    overlap = shapely.intersection(geom, changed_union)
                    touches = not overlap.is_empty
                district_touches[dname] = touches
            else:
                district_touches[dname] = False
    
    def get_touches(dn_str):
        if pd.isna(dn_str) or dn_str == "unknown":
            return False
        for d in dn_str.split(","):
            d = d.strip()
            if district_touches.get(d, False):
                return True
        return False
    
    result_gdf["touches_changed"] = result_gdf["district_name"].apply(get_touches)
    
    # Simplify geometries slightly
    result_gdf["geometry"] = result_gdf["geometry"].simplify(0.00001, preserve_topology=True)
    result_gdf["geometry"] = result_gdf["geometry"].buffer(0)
    
    # Ensure MultiPolygon type
    def ensure_multipolygon(geom):
        if geom.geom_type == "Polygon":
            return shapely.geometry.MultiPolygon([geom])
        return geom
    
    result_gdf["geometry"] = result_gdf["geometry"].apply(ensure_multipolygon)
    
    # -----------------------------------------------------------------------
    # 5. Output
    # -----------------------------------------------------------------------
    output_path = "/work/vienna_boundary_changes.geojson"
    result_gdf.to_file(output_path, driver="GeoJSON")
    
    print(f"\n{'=' * 60}")
    print(f"Output written to {output_path}")
    print(f"Features: {len(result_gdf)}")
    print(f"Columns: {result_gdf.columns.tolist()}")
    print(f"\nSample features:")
    print(result_gdf[["change_type", "district_name", "touches_changed"]].head(15))
    
    print(f"\nDistrict touch summary:")
    for dname, touches in sorted(district_touches.items()):
        print(f"  {dname}: touches_changed={touches}")
    
    # Verify output
    verify = gpd.read_file(output_path)
    print(f"\nVerification: read back {len(verify)} features")
    print(verify[["change_type", "district_name", "touches_changed"]].head(10))
    print(f"Geometry types: {verify.geometry.geom_type.unique()}")


if __name__ == "__main__":
    main()

"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass via osmnx,
compare to find changes, and output vienna_boundary_changes.geojson.
"""

import os
import sys
import json
from collections import OrderedDict

import osmnx as ox
import osmnx._overpass as op
import geopandas as gpd
import pandas as pd
import shapely
import shapely.geometry
import numpy as np

ox.settings.log_console = False

# ---------------------------------------------------------------------------
# 1. Overpass query
# ---------------------------------------------------------------------------

def assemble_ring_from_ways(way_refs, ways_dict, nodes_dict):
    """Assemble a closed ring from a sequence of OSM way references.
    
    Returns (ring_coords, used_mask) or (None, None) if no ring could be built.
    """
    if not way_refs:
        return None, None
    
    segments = []
    for ref in way_refs:
        if ref not in ways_dict:
            continue
        node_ids = ways_dict[ref]
        if len(node_ids) < 2:
            continue
        coords = [nodes_dict.get(nid) for nid in node_ids]
        coords = [c for c in coords if c is not None]
        if len(coords) < 2:
            continue
        segments.append({
            'ref': ref,
            'coords': list(coords),
            'start_node': node_ids[0],
            'end_node': node_ids[-1],
        })
    
    if not segments:
        return None, None
    
    used = [False] * len(segments)
    ring_coords = list(segments[0]['coords'])
    current_end_node = segments[0]['end_node']
    used[0] = True
    
    changed = True
    while changed:
        changed = False
        for i, seg in enumerate(segments):
            if used[i]:
                continue
            if seg['start_node'] == current_end_node:
                ring_coords.extend(seg['coords'][1:])
                current_end_node = seg['end_node']
                used[i] = True
                changed = True
                break
            elif seg['end_node'] == current_end_node:
                rev_coords = list(reversed(seg['coords']))
                ring_coords.extend(rev_coords[1:])
                current_end_node = seg['start_node']
                used[i] = True
                changed = True
                break
    
    if ring_coords[0] != ring_coords[-1]:
        ring_coords.append(ring_coords[0])
    
    return ring_coords, used


def fetch_vienna_districts_osmnx(date=None):
    """Fetch Vienna district (admin_level=9) boundary relations via Overpass."""
    bbox = (48.1, 16.1, 48.35, 16.6)
    
    if date:
        query = f"""
[out:json][timeout:180][date:"{date}"];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    else:
        query = f"""
[out:json][timeout:180];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    
    data = OrderedDict()
    data['data'] = query
    result = op._overpass_request(data)
    elements = result.get("elements", [])
    
    nodes = {}
    ways = {}
    relations = []
    
    for el in elements:
        etype = el.get("type")
        eid = el.get("id")
        if etype == "node":
            nodes[eid] = (el.get("lon"), el.get("lat"))
        elif etype == "way":
            ways[eid] = el.get("nodes", [])
        elif etype == "relation":
            relations.append(el)
    
    features = []
    
    for rel in relations:
        tags = rel.get("tags", {})
        name = tags.get("name", "")
        admin_level = tags.get("admin_level", "")
        
        if admin_level != "9":
            continue
        
        members = rel.get("members", [])
        
        # Collect way refs by role
        all_outer_refs = []
        inner_refs_coords = []
        
        for mem in members:
            if mem.get("type") != "way":
                continue
            role = mem.get("role", "")
            ref = mem.get("ref")
            if ref not in ways:
                continue
            if role == "inner":
                node_ids = ways[ref]
                coords = [nodes.get(nid) for nid in node_ids]
                coords = [c for c in coords if c is not None]
                if len(coords) >= 3:
                    if coords[0] != coords[-1]:
                        coords.append(coords[0])
                    inner_refs_coords.append(coords)
            else:
                all_outer_refs.append(ref)
        
        # Build outer rings iteratively
        outer_polygons = []
        remaining = list(all_outer_refs)
        
        while remaining:
            ring_coords, used = assemble_ring_from_ways(remaining, ways, nodes)
            if ring_coords is None or len(ring_coords) < 4:
                break
            
            try:
                ls = shapely.LineString(ring_coords)
                if not ls.is_valid or ls.is_empty or ls.length < 1e-10:
                    break
                poly = shapely.Polygon(ring_coords)
                poly = shapely.polygon.orient(poly)
                if poly.is_valid and not poly.is_empty and poly.area > 1e-12:
                    outer_polygons.append(poly)
            except Exception:
                break
            
            # Remove used refs from remaining
            remaining = [r for i, r in enumerate(remaining) 
                        if i >= len(used) or not used[i]]
        
        if not outer_polygons:
            continue
        
        # Handle inner rings
        inner_polygons = []
        for coords in inner_refs_coords:
            try:
                inner_ls = shapely.LineString(coords)
                if inner_ls.is_valid and not inner_ls.is_empty and inner_ls.length > 1e-10:
                    ipoly = shapely.Polygon(coords)
                    ipoly = shapely.polygon.orient(ipoly)
                    if ipoly.is_valid and not ipoly.is_empty:
                        inner_polygons.append(ipoly)
            except Exception:
                continue
        
        # Cut inner rings from outer polygons
        final_polygons = []
        for outer_poly in outer_polygons:
            inners_here = []
            rest = []
            for ipoly in inner_polygons:
                try:
                    if outer_poly.contains(ipoly.representative_point()) or outer_poly.intersects(ipoly):
                        inners_here.append(ipoly)
                    else:
                        rest.append(ipoly)
                except Exception:
                    rest.append(ipoly)
            inner_polygons = rest
            
            if inners_here:
                try:
                    outer_coords = list(outer_poly.exterior.coords)
                    inner_coords_list = [list(ip.exterior.coords) for ip in inners_here]
                    poly = shapely.Polygon(outer_coords, inner_coords_list)
                    poly = shapely.polygon.orient(poly)
                    if poly.is_valid and not poly.is_empty:
                        final_polygons.append(poly)
                    else:
                        final_polygons.append(outer_poly)
                except Exception:
                    final_polygons.append(outer_poly)
            else:
                final_polygons.append(outer_poly)
        
        if final_polygons:
            geom = shapely.unary_union(final_polygons)
            if geom.geom_type == "Polygon":
                geom = shapely.geometry.MultiPolygon([geom])
            
            if not geom.is_empty:
                features.append({
                    "geometry": geom,
                    "name": name,
                    "osm_id": rel.get("id"),
                })
    
    if not features:
        return gpd.GeoDataFrame()
    
    return gpd.GeoDataFrame(features, crs="EPSG:4326")


# ---------------------------------------------------------------------------
# 2. Name normalisation
# ---------------------------------------------------------------------------

VIENNA_DISTRICTS_CANON = {
    "innere stadt": "Innere Stadt",
    "leopoldstadt": "Leopoldstadt",
    "landstraße": "Landstraße",
    "landstrasse": "Landstraße",
    "wieden": "Wieden",
    "margareten": "Margareten",
    "margarethen": "Margareten",
    "mariahilf": "Mariahilf",
    "neubau": "Neubau",
    "josefstadt": "Josefstadt",
    "alsergrund": "Alsergrund",
    "favoriten": "Favoriten",
    "simmering": "Simmering",
    "meidling": "Meidling",
    "hietzing": "Hietzing",
    "penzing": "Penzing",
    "rudolfsheim-fünfhaus": "Rudolfsheim-Fünfhaus",
    "rudolfsheim-fuenfhaus": "Rudolfsheim-Fünfhaus",
    "ottakring": "Ottakring",
    "hernals": "Hernals",
    "währing": "Währing",
    "waehring": "Währing",
    "döbling": "Döbling",
    "doebling": "Döbling",
    "brigittenau": "Brigittenau",
    "floridsdorf": "Floridsdorf",
    "donaustadt": "Donaustadt",
    "liesing": "Liesing",
}


def normalise_name(name):
    if not name or pd.isna(name):
        return None
    n = name.strip()
    key = n.lower().replace("  ", " ").replace("ß", "ss").replace("-", "-")
    if key in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key]
    key2 = key.replace("ss", "ß")
    if key2 in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key2]
    for pattern, canon in VIENNA_DISTRICTS_CANON.items():
        if pattern in key or key in pattern:
            return canon
    return n


# ---------------------------------------------------------------------------
# 3. Main pipeline
# ---------------------------------------------------------------------------

def main():
    print("=" * 60)
    print("Fetching CURRENT Vienna district boundaries...")
    print("=" * 60)
    current_gdf = fetch_vienna_districts_osmnx()
    print(f"  Got {len(current_gdf)} current district features")
    
    print("\n" + "=" * 60)
    print("Fetching 2014-01-01 Vienna district boundaries...")
    print("=" * 60)
    historical_gdf = fetch_vienna_districts_osmnx(date="2014-01-01T00:00:00Z")
    print(f"  Got {len(historical_gdf)} historical district features")
    
    if current_gdf.empty or historical_gdf.empty:
        print("ERROR: Empty data returned!")
        sys.exit(1)
    
    # Normalise names
    current_gdf["district_name"] = current_gdf["name"].apply(normalise_name)
    historical_gdf["district_name"] = historical_gdf["name"].apply(normalise_name)
    
    print("\nCurrent district names:")
    for n in sorted(current_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    print("\nHistorical district names:")
    for n in sorted(historical_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    
    # Filter to Vienna districts only
    valid_names = set(VIENNA_DISTRICTS_CANON.values())
    current_gdf = current_gdf[current_gdf["district_name"].isin(valid_names)].copy()
    historical_gdf = historical_gdf[historical_gdf["district_name"].isin(valid_names)].copy()
    
    print(f"\nAfter filtering:")
    print(f"  Current: {len(current_gdf)} districts")
    print(f"  Historical: {len(historical_gdf)} districts")
    
    current_names = set(current_gdf["district_name"].dropna())
    historical_names = set(historical_gdf["district_name"].dropna())
    
    common = current_names & historical_names
    only_current = current_names - historical_names
    only_historical = historical_names - current_names
    
    print(f"\nCommon districts: {len(common)}")
    if only_current:
        print(f"Only in current: {only_current}")
    if only_historical:
        print(f"Only in historical: {only_historical}")
    
    # Build unified dissolve
    current_union = current_gdf.geometry.unary_union
    historical_union = historical_gdf.geometry.unary_union
    
    print(f"\nCurrent total area: {current_union.area:.8f} sq deg")
    print(f"Historical total area: {historical_union.area:.8f} sq deg")
    
    # Symmetric difference
    added = shapely.difference(current_union, historical_union)
    removed = shapely.difference(historical_union, current_union)
    unchanged = shapely.intersection(current_union, historical_union)
    
    for label, geom in [("Added", added), ("Removed", removed), ("Unchanged", unchanged)]:
        if geom.is_empty:
            print(f"{label}: empty")
        else:
            print(f"{label}: area={geom.area:.8f}, type={geom.geom_type}")
    
    # Collect fragments
    fragments = []
    
    def add_fragments(geom, ctype, dist_name_fn):
        if geom.is_empty:
            return
        if geom.geom_type == "MultiPolygon":
            polys = list(geom.geoms)
        elif geom.geom_type == "Polygon":
            polys = [geom]
        elif geom.geom_type == "GeometryCollection":
            polys = [g for g in geom.geoms if g.geom_type in ("Polygon", "MultiPolygon")]
        else:
            polys = [geom]
        
        for poly in polys:
            if poly.is_empty:
                continue
            if poly.geom_type == "Polygon":
                mp = shapely.geometry.MultiPolygon([poly])
            elif poly.geom_type == "MultiPolygon":
                mp = poly
            else:
                continue
            dname = dist_name_fn(mp)
            fragments.append({
                "geometry": mp,
                "change_type": ctype,
                "district_name": dname,
            })
    
    def added_dist_name(geom):
        districts = set()
        for _, row in current_gdf.iterrows():
            if geom.intersects(row["geometry"]):
                districts.add(row["district_name"])
        return ", ".join(sorted(districts)) if districts else "unknown"
    
    add_fragments(added, "added_since_2014", added_dist_name)
    
    def removed_dist_name(geom):
        districts = set()
        for _, row in historical_gdf.iterrows():
            if geom.intersects(row["geometry"]):
                districts.add(row["district_name"])
        return ", ".join(sorted(districts)) if districts else "unknown"
    
    add_fragments(removed, "removed_since_2014", removed_dist_name)
    
    # Unchanged fragments per current district
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        dist_geom = row["geometry"]
        unchanged_part = shapely.intersection(dist_geom, unchanged)
        if not unchanged_part.is_empty:
            def dn_fn(geom, dn=dname):
                return dn
            add_fragments(unchanged_part, "unchanged", dn_fn)
    
    result_gdf = gpd.GeoDataFrame(fragments, crs="EPSG:4326")
    
    if result_gdf.empty:
        print("No fragments — all boundaries unchanged.")
        result_gdf = current_gdf.copy()
        result_gdf["change_type"] = "unchanged"
        result_gdf = result_gdf[["geometry", "change_type", "district_name"]]
    
    print(f"\nFragments created: {len(result_gdf)}")
    if not result_gdf.empty:
        print(result_gdf["change_type"].value_counts())
    
    # -----------------------------------------------------------------------
    # 4. touches_changed
    # -----------------------------------------------------------------------
    
    changed_geoms = []
    for _, row in result_gdf.iterrows():
        if row["change_type"] in ("added_since_2014", "removed_since_2014"):
            changed_geoms.append(row["geometry"])
    
    if changed_geoms:
        changed_union = shapely.unary_union(changed_geoms)
    else:
        changed_union = shapely.geometry.MultiPolygon()
    
    district_touches = {}
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        geom = row["geometry"]
        if not changed_union.is_empty:
            touches = geom.touches(changed_union)
            if not touches:
                overlap = shapely.intersection(geom, changed_union)
                touches = not overlap.is_empty
            district_touches[dname] = touches
        else:
            district_touches[dname] = False
    
    for _, row in historical_gdf.iterrows():
        dname = row["district_name"]
        if dname not in district_touches:
            geom = row["geometry"]
            if not changed_union.is_empty:
                touches = geom.touches(changed_union)
                if not touches:
                    overlap = shapely.intersection(geom, changed_union)
                    touches = not overlap.is_empty
                district_touches[dname] = touches
            else:
                district_touches[dname] = False
    
    def get_touches(dn_str):
        if pd.isna(dn_str) or dn_str == "unknown":
            return False
        for d in dn_str.split(","):
            d = d.strip()
            if district_touches.get(d, False):
                return True
        return False
    
    result_gdf["touches_changed"] = result_gdf["district_name"].apply(get_touches)
    
    # Simplify and ensure MultiPolygon
    result_gdf["geometry"] = result_gdf["geometry"].simplify(0.00001, preserve_topology=True)
    result_gdf["geometry"] = result_gdf["geometry"].buffer(0)
    
    def ensure_multipolygon(geom):
        if geom.geom_type == "Polygon":
            return shapely.geometry.MultiPolygon([geom])
        return geom
    
    result_gdf["geometry"] = result_gdf["geometry"].apply(ensure_multipolygon)
    
    # -----------------------------------------------------------------------
    # 5. Output
    # -----------------------------------------------------------------------
    output_path = "/work/vienna_boundary_changes.geojson"
    result_gdf.to_file(output_path, driver="GeoJSON")
    
    print(f"\n{'=' * 60}")
    print(f"Output written to {output_path}")
    print(f"Features: {len(result_gdf)}")
    print(f"Columns: {result_gdf.columns.tolist()}")
    print(f"\nSample features:")
    print(result_gdf[["change_type", "district_name", "touches_changed"]].head(15))
    
    print(f"\nDistrict touch summary:")
    for dname, touches in sorted(district_touches.items()):
        print(f"  {dname}: touches_changed={touches}")
    
    # Verify
    verify = gpd.read_file(output_path)
    print(f"\nVerification: read back {len(verify)} features")
    print(verify[["change_type", "district_name", "touches_changed"]].head(10))
    print(f"Geometry types: {verify.geometry.geom_type.unique()}")


if __name__ == "__main__":
    main()

"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass via osmnx,
compare to find changes, and output vienna_boundary_changes.geojson.
"""

import os
import sys
import json
from collections import OrderedDict

import osmnx as ox
import osmnx._overpass as op
import geopandas as gpd
import pandas as pd
import shapely
import shapely.geometry
import numpy as np

ox.settings.log_console = False

# ---------------------------------------------------------------------------
# 1. Overpass query
# ---------------------------------------------------------------------------

def build_ring_from_ways(way_refs, ways_dict, nodes_dict):
    """Assemble a closed ring from a sequence of OSM way references.
    
    Returns (ring_coords, used_refs) or (None, []).
    """
    if not way_refs:
        return None, []
    
    segments = []
    for ref in way_refs:
        if ref not in ways_dict:
            continue
        node_ids = ways_dict[ref]
        if len(node_ids) < 2:
            continue
        coords = [nodes_dict.get(nid) for nid in node_ids]
        coords = [c for c in coords if c is not None]
        if len(coords) < 2:
            continue
        segments.append({
            'ref': ref,
            'coords': list(coords),
            'start_node': node_ids[0],
            'end_node': node_ids[-1],
        })
    
    if not segments:
        return None, []
    
    used = [False] * len(segments)
    
    # Start with longest segment
    longest_idx = max(range(len(segments)), key=lambda i: len(segments[i]['coords']))
    ring_coords = list(segments[longest_idx]['coords'])
    current_end_node = segments[longest_idx]['end_node']
    used[longest_idx] = True
    
    changed = True
    while changed:
        changed = False
        for i, seg in enumerate(segments):
            if used[i]:
                continue
            if seg['start_node'] == current_end_node:
                ring_coords.extend(seg['coords'][1:])
                current_end_node = seg['end_node']
                used[i] = True
                changed = True
                break
            elif seg['end_node'] == current_end_node:
                rev_coords = list(reversed(seg['coords']))
                ring_coords.extend(rev_coords[1:])
                current_end_node = seg['start_node']
                used[i] = True
                changed = True
                break
    
    if ring_coords[0] != ring_coords[-1]:
        ring_coords.append(ring_coords[0])
    
    used_refs = [segments[i]['ref'] for i in range(len(segments)) if used[i]]
    return ring_coords, used_refs


def fetch_vienna_districts_osmnx(date=None):
    """Fetch Vienna district (admin_level=9) boundary relations via Overpass."""
    bbox = (48.1, 16.1, 48.35, 16.6)
    
    if date:
        query = f"""
[out:json][timeout:180][date:"{date}"];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    else:
        query = f"""
[out:json][timeout:180];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    
    data = OrderedDict()
    data['data'] = query
    result = op._overpass_request(data)
    elements = result.get("elements", [])
    
    nodes = {}
    ways = {}
    relations = []
    
    for el in elements:
        etype = el.get("type")
        eid = el.get("id")
        if etype == "node":
            nodes[eid] = (el.get("lon"), el.get("lat"))
        elif etype == "way":
            ways[eid] = el.get("nodes", [])
        elif etype == "relation":
            relations.append(el)
    
    features = []
    
    for rel in relations:
        tags = rel.get("tags", {})
        name = tags.get("name", "")
        admin_level = tags.get("admin_level", "")
        
        if admin_level != "9":
            continue
        
        members = rel.get("members", [])
        
        all_outer_refs = []
        inner_refs_coords = []
        
        for mem in members:
            if mem.get("type") != "way":
                continue
            role = mem.get("role", "")
            ref = mem.get("ref")
            if ref not in ways:
                continue
            if role == "inner":
                node_ids = ways[ref]
                coords = [nodes.get(nid) for nid in node_ids]
                coords = [c for c in coords if c is not None]
                if len(coords) >= 3:
                    if coords[0] != coords[-1]:
                        coords.append(coords[0])
                    inner_refs_coords.append(coords)
            else:
                all_outer_refs.append(ref)
        
        # Build outer polygons iteratively
        outer_polygons = []
        remaining = list(all_outer_refs)
        
        while remaining:
            ring_coords, used_refs = build_ring_from_ways(remaining, ways, nodes)
            if ring_coords is None or len(ring_coords) < 4:
                break
            
            try:
                ls = shapely.LineString(ring_coords)
                if not ls.is_valid or ls.is_empty or ls.length < 1e-10:
                    break
                poly = shapely.Polygon(ring_coords)
                poly = shapely.orient_polygons(poly)
                if poly.is_valid and not poly.is_empty and poly.area > 1e-12:
                    outer_polygons.append(poly)
            except Exception:
                break
            
            remaining = [r for r in remaining if r not in used_refs]
        
        if not outer_polygons:
            continue
        
        # Build inner polygons
        inner_polygons = []
        for coords in inner_refs_coords:
            try:
                ipoly = shapely.Polygon(coords)
                ipoly = shapely.orient_polygons(ipoly)
                if ipoly.is_valid and not ipoly.is_empty:
                    inner_polygons.append(ipoly)
            except Exception:
                continue
        
        # Combine outer with holes
        final_polygons = []
        for outer_poly in outer_polygons:
            inners_here = []
            rest = []
            for ipoly in inner_polygons:
                try:
                    if outer_poly.contains(ipoly.representative_point()) or outer_poly.intersects(ipoly):
                        inners_here.append(ipoly)
                    else:
                        rest.append(ipoly)
                except Exception:
                    rest.append(ipoly)
            inner_polygons = rest
            
            if inners_here:
                try:
                    outer_coords = list(outer_poly.exterior.coords)
                    inner_coords_list = [list(ip.exterior.coords) for ip in inners_here]
                    poly = shapely.Polygon(outer_coords, inner_coords_list)
                    poly = shapely.orient_polygons(poly)
                    if poly.is_valid and not poly.is_empty:
                        final_polygons.append(poly)
                    else:
                        final_polygons.append(outer_poly)
                except Exception:
                    final_polygons.append(outer_poly)
            else:
                final_polygons.append(outer_poly)
        
        if final_polygons:
            geom = shapely.union_all(final_polygons)
            if geom.geom_type == "Polygon":
                geom = shapely.geometry.MultiPolygon([geom])
            
            if not geom.is_empty:
                features.append({
                    "geometry": geom,
                    "name": name,
                    "osm_id": rel.get("id"),
                })
    
    if not features:
        return gpd.GeoDataFrame()
    
    return gpd.GeoDataFrame(features, crs="EPSG:4326")


# ---------------------------------------------------------------------------
# 2. Name normalisation
# ---------------------------------------------------------------------------

VIENNA_DISTRICTS_CANON = {
    "innere stadt": "Innere Stadt",
    "leopoldstadt": "Leopoldstadt",
    "landstraße": "Landstraße",
    "landstrasse": "Landstraße",
    "wieden": "Wieden",
    "margareten": "Margareten",
    "margarethen": "Margareten",
    "mariahilf": "Mariahilf",
    "neubau": "Neubau",
    "josefstadt": "Josefstadt",
    "alsergrund": "Alsergrund",
    "favoriten": "Favoriten",
    "simmering": "Simmering",
    "meidling": "Meidling",
    "hietzing": "Hietzing",
    "penzing": "Penzing",
    "rudolfsheim-fünfhaus": "Rudolfsheim-Fünfhaus",
    "rudolfsheim-fuenfhaus": "Rudolfsheim-Fünfhaus",
    "ottakring": "Ottakring",
    "hernals": "Hernals",
    "währing": "Währing",
    "waehring": "Währing",
    "döbling": "Döbling",
    "doebling": "Döbling",
    "brigittenau": "Brigittenau",
    "floridsdorf": "Floridsdorf",
    "donaustadt": "Donaustadt",
    "liesing": "Liesing",
}


def normalise_name(name):
    if not name or pd.isna(name):
        return None
    n = name.strip()
    key = n.lower().replace("  ", " ").replace("ß", "ss").replace("-", "-")
    if key in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key]
    key2 = key.replace("ss", "ß")
    if key2 in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key2]
    for pattern, canon in VIENNA_DISTRICTS_CANON.items():
        if pattern in key or key in pattern:
            return canon
    return n


# ---------------------------------------------------------------------------
# 3. Main pipeline
# ---------------------------------------------------------------------------

def main():
    print("=" * 60)
    print("Fetching CURRENT Vienna district boundaries...")
    print("=" * 60)
    current_gdf = fetch_vienna_districts_osmnx()
    print(f"  Got {len(current_gdf)} current district features")
    
    print("\n" + "=" * 60)
    print("Fetching 2014-01-01 Vienna district boundaries...")
    print("=" * 60)
    historical_gdf = fetch_vienna_districts_osmnx(date="2014-01-01T00:00:00Z")
    print(f"  Got {len(historical_gdf)} historical district features")
    
    if current_gdf.empty or historical_gdf.empty:
        print("ERROR: Empty data returned!")
        sys.exit(1)
    
    # Normalise names
    current_gdf["district_name"] = current_gdf["name"].apply(normalise_name)
    historical_gdf["district_name"] = historical_gdf["name"].apply(normalise_name)
    
    print("\nCurrent district names:")
    for n in sorted(current_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    print("\nHistorical district names:")
    for n in sorted(historical_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    
    # Filter to Vienna districts only
    valid_names = set(VIENNA_DISTRICTS_CANON.values())
    current_gdf = current_gdf[current_gdf["district_name"].isin(valid_names)].copy()
    historical_gdf = historical_gdf[historical_gdf["district_name"].isin(valid_names)].copy()
    
    print(f"\nAfter filtering:")
    print(f"  Current: {len(current_gdf)} districts")
    print(f"  Historical: {len(historical_gdf)} districts")
    
    # Build per-district comparison
    # Index historical by name
    hist_by_name = {}
    for _, row in historical_gdf.iterrows():
        hist_by_name[row["district_name"]] = row["geometry"]
    
    # Tolerance for snapping (roughly 1-2 meters in degrees)
    TOLERANCE = 1e-7
    
    fragments = []
    
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        current_geom = row["geometry"]
        historical_geom = hist_by_name.get(dname)
        
        if historical_geom is None:
            # District exists now but not in 2014 — fully added
            fragments.append({
                "geometry": current_geom,
                "change_type": "added_since_2014",
                "district_name": dname,
            })
            continue
        
        # Snap current to historical to reduce slivers
        # Use normalize for robust comparison
        current_norm = shapely.normalize(current_geom)
        hist_norm = shapely.normalize(historical_geom)
        
        # Snap historical geometry to current geometry's coordinates
        snapped_hist = shapely.snap(hist_norm, current_norm, TOLERANCE)
        
        # Compute differences
        added_part = shapely.difference(current_norm, snapped_hist)
        removed_part = shapely.difference(snapped_hist, current_norm)
        unchanged_part = shapely.intersection(current_norm, snapped_hist)
        
        # Clean up small slivers (noise)
        MIN_AREA = 1e-10  # ~ 1 sq meter in degrees
        
        def clean_geom(geom):
            """Remove tiny polygons and ensure MultiPolygon."""
            if geom.is_empty:
                return shapely.geometry.MultiPolygon()
            
            # Filter out tiny polygons
            if geom.geom_type == "MultiPolygon":
                polys = [g for g in geom.geoms if g.area >= MIN_AREA]
                if not polys:
                    return shapely.geometry.MultiPolygon()
                if len(polys) == 1:
                    geom = polys[0]
                else:
                    geom = shapely.geometry.MultiPolygon(polys)
            
            if geom.geom_type == "Polygon":
                if geom.area < MIN_AREA:
                    return shapely.geometry.MultiPolygon()
                return shapely.geometry.MultiPolygon([geom])
            
            return geom
        
        added_clean = clean_geom(added_part)
        removed_clean = clean_geom(removed_part)
        unchanged_clean = clean_geom(unchanged_part)
        
        if not added_clean.is_empty:
            if added_clean.geom_type == "MultiPolygon":
                for poly in added_clean.geoms:
                    fragments.append({
                        "geometry": shapely.geometry.MultiPolygon([poly]),
                        "change_type": "added_since_2014",
                        "district_name": dname,
                    })
            else:
                fragments.append({
                    "geometry": added_clean,
                    "change_type": "added_since_2014",
                    "district_name": dname,
                })
        
        if not removed_clean.is_empty:
            if removed_clean.geom_type == "MultiPolygon":
                for poly in removed_clean.geoms:
                    fragments.append({
                        "geometry": shapely.geometry.MultiPolygon([poly]),
                        "change_type": "removed_since_2014",
                        "district_name": dname,
                    })
            else:
                fragments.append({
                    "geometry": removed_clean,
                    "change_type": "removed_since_2014",
                    "district_name": dname,
                })
        
        if not unchanged_clean.is_empty:
            if unchanged_clean.geom_type == "MultiPolygon":
                for poly in unchanged_clean.geoms:
                    fragments.append({
                        "geometry": shapely.geometry.MultiPolygon([poly]),
                        "change_type": "unchanged",
                        "district_name": dname,
                    })
            else:
                fragments.append({
                    "geometry": unchanged_clean,
                    "change_type": "unchanged",
                    "district_name": dname,
                })
    
    result_gdf = gpd.GeoDataFrame(fragments, crs="EPSG:4326")
    
    if result_gdf.empty:
        print("No fragments produced. Using current boundaries as unchanged.")
        result_gdf = current_gdf.copy()
        result_gdf["change_type"] = "unchanged"
        result_gdf = result_gdf[["geometry", "change_type", "district_name"]]
    
    print(f"\nFragments created: {len(result_gdf)}")
    if not result_gdf.empty:
        print(result_gdf["change_type"].value_counts())
    
    # -----------------------------------------------------------------------
    # 4. touches_changed
    # -----------------------------------------------------------------------
    
    changed_geoms = []
    for _, row in result_gdf.iterrows():
        if row["change_type"] in ("added_since_2014", "removed_since_2014"):
            changed_geoms.append(row["geometry"])
    
    if changed_geoms:
        changed_union = shapely.union_all(changed_geoms)
    else:
        changed_union = shapely.geometry.MultiPolygon()
    
    # Compute touches per district (using current boundaries as reference)
    district_touches = {}
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        geom = row["geometry"]
        if not changed_union.is_empty:
            touches = geom.touches(changed_union)
            if not touches:
                overlap = shapely.intersection(geom, changed_union)
                touches = not overlap.is_empty
            district_touches[dname] = touches
        else:
            district_touches[dname] = False
    
    # Also check historical-only districts
    for _, row in historical_gdf.iterrows():
        dname = row["district_name"]
        if dname not in district_touches:
            geom = row["geometry"]
            if not changed_union.is_empty:
                touches = geom.touches(changed_union)
                if not touches:
                    overlap = shapely.intersection(geom, changed_union)
                    touches = not overlap.is_empty
                district_touches[dname] = touches
            else:
                district_touches[dname] = False
    
    def get_touches(dn_str):
        if pd.isna(dn_str) or dn_str == "unknown":
            return False
        for d in dn_str.split(","):
            d = d.strip()
            if district_touches.get(d, False):
                return True
        return False
    
    result_gdf["touches_changed"] = result_gdf["district_name"].apply(get_touches)
    
    # Simplify geometries and ensure MultiPolygon
    result_gdf["geometry"] = result_gdf["geometry"].simplify(0.00001, preserve_topology=True)
    result_gdf["geometry"] = result_gdf["geometry"].buffer(0)
    
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == "Polygon":
            return shapely.geometry.MultiPolygon([geom])
        return geom
    
    result_gdf["geometry"] = result_gdf["geometry"].apply(ensure_multipolygon)
    
    # -----------------------------------------------------------------------
    # 5. Output
    # -----------------------------------------------------------------------
    output_path = "/work/vienna_boundary_changes.geojson"
    result_gdf.to_file(output_path, driver="GeoJSON")
    
    print(f"\n{'=' * 60}")
    print(f"Output written to {output_path}")
    print(f"Features: {len(result_gdf)}")
    print(f"Columns: {result_gdf.columns.tolist()}")
    
    print(f"\nSample features:")
    print(result_gdf[["change_type", "district_name", "touches_changed"]].head(20))
    
    print(f"\nDistrict touch summary:")
    for dname, touches in sorted(district_touches.items()):
        print(f"  {dname}: touches_changed={touches}")
    
    # Verify
    verify = gpd.read_file(output_path)
    print(f"\nVerification: read back {len(verify)} features")
    print(verify["change_type"].value_counts())
    print(f"Geometry types: {verify.geometry.geom_type.unique()}")
    
    # Area stats
    for ct in ["added_since_2014", "removed_since_2014", "unchanged"]:
        sub = verify[verify["change_type"] == ct]
        if not sub.empty:
            areas = sub.geometry.area
            print(f"  {ct}: count={len(sub)}, total_area={areas.sum():.10f}")


if __name__ == "__main__":
    main()

"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass via osmnx,
compare to find changes, and output vienna_boundary_changes.geojson.
"""

import os
import sys
import json
from collections import OrderedDict

import osmnx as ox
import osmnx._overpass as op
import geopandas as gpd
import pandas as pd
import shapely
import shapely.geometry
import numpy as np

ox.settings.log_console = False

# ---------------------------------------------------------------------------
# 1. Overpass query
# ---------------------------------------------------------------------------

def build_ring_from_ways(way_refs, ways_dict, nodes_dict):
    if not way_refs:
        return None, []
    
    segments = []
    for ref in way_refs:
        if ref not in ways_dict:
            continue
        node_ids = ways_dict[ref]
        if len(node_ids) < 2:
            continue
        coords = [nodes_dict.get(nid) for nid in node_ids]
        coords = [c for c in coords if c is not None]
        if len(coords) < 2:
            continue
        segments.append({
            'ref': ref,
            'coords': list(coords),
            'start_node': node_ids[0],
            'end_node': node_ids[-1],
        })
    
    if not segments:
        return None, []
    
    used = [False] * len(segments)
    longest_idx = max(range(len(segments)), key=lambda i: len(segments[i]['coords']))
    ring_coords = list(segments[longest_idx]['coords'])
    current_end_node = segments[longest_idx]['end_node']
    used[longest_idx] = True
    
    changed = True
    while changed:
        changed = False
        for i, seg in enumerate(segments):
            if used[i]:
                continue
            if seg['start_node'] == current_end_node:
                ring_coords.extend(seg['coords'][1:])
                current_end_node = seg['end_node']
                used[i] = True
                changed = True
                break
            elif seg['end_node'] == current_end_node:
                rev_coords = list(reversed(seg['coords']))
                ring_coords.extend(rev_coords[1:])
                current_end_node = seg['start_node']
                used[i] = True
                changed = True
                break
    
    if ring_coords[0] != ring_coords[-1]:
        ring_coords.append(ring_coords[0])
    
    used_refs = [segments[i]['ref'] for i in range(len(segments)) if used[i]]
    return ring_coords, used_refs


def fetch_vienna_districts_osmnx(date=None):
    """Fetch Vienna district (admin_level=9) boundary relations via Overpass."""
    bbox = (48.1, 16.1, 48.35, 16.6)
    
    if date:
        query = f"""
[out:json][timeout:180][date:"{date}"];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    else:
        query = f"""
[out:json][timeout:180];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    
    data = OrderedDict()
    data['data'] = query
    result = op._overpass_request(data)
    elements = result.get("elements", [])
    
    nodes = {}
    ways = {}
    relations = []
    
    for el in elements:
        etype = el.get("type")
        eid = el.get("id")
        if etype == "node":
            nodes[eid] = (el.get("lon"), el.get("lat"))
        elif etype == "way":
            ways[eid] = el.get("nodes", [])
        elif etype == "relation":
            relations.append(el)
    
    features = []
    
    for rel in relations:
        tags = rel.get("tags", {})
        name = tags.get("name", "")
        admin_level = tags.get("admin_level", "")
        
        if admin_level != "9":
            continue
        
        members = rel.get("members", [])
        
        all_outer_refs = []
        inner_refs_coords = []
        
        for mem in members:
            if mem.get("type") != "way":
                continue
            role = mem.get("role", "")
            ref = mem.get("ref")
            if ref not in ways:
                continue
            if role == "inner":
                node_ids = ways[ref]
                coords = [nodes.get(nid) for nid in node_ids]
                coords = [c for c in coords if c is not None]
                if len(coords) >= 3:
                    if coords[0] != coords[-1]:
                        coords.append(coords[0])
                    inner_refs_coords.append(coords)
            else:
                all_outer_refs.append(ref)
        
        outer_polygons = []
        remaining = list(all_outer_refs)
        
        while remaining:
            ring_coords, used_refs = build_ring_from_ways(remaining, ways, nodes)
            if ring_coords is None or len(ring_coords) < 4:
                break
            
            try:
                ls = shapely.LineString(ring_coords)
                if not ls.is_valid or ls.is_empty or ls.length < 1e-10:
                    break
                poly = shapely.Polygon(ring_coords)
                poly = shapely.orient_polygons(poly)
                if poly.is_valid and not poly.is_empty and poly.area > 1e-12:
                    outer_polygons.append(poly)
            except Exception:
                break
            
            remaining = [r for r in remaining if r not in used_refs]
        
        if not outer_polygons:
            continue
        
        inner_polygons = []
        for coords in inner_refs_coords:
            try:
                ipoly = shapely.Polygon(coords)
                ipoly = shapely.orient_polygons(ipoly)
                if ipoly.is_valid and not ipoly.is_empty:
                    inner_polygons.append(ipoly)
            except Exception:
                continue
        
        final_polygons = []
        for outer_poly in outer_polygons:
            inners_here = []
            rest = []
            for ipoly in inner_polygons:
                try:
                    if outer_poly.contains(ipoly.representative_point()) or outer_poly.intersects(ipoly):
                        inners_here.append(ipoly)
                    else:
                        rest.append(ipoly)
                except Exception:
                    rest.append(ipoly)
            inner_polygons = rest
            
            if inners_here:
                try:
                    outer_coords = list(outer_poly.exterior.coords)
                    inner_coords_list = [list(ip.exterior.coords) for ip in inners_here]
                    poly = shapely.Polygon(outer_coords, inner_coords_list)
                    poly = shapely.orient_polygons(poly)
                    if poly.is_valid and not poly.is_empty:
                        final_polygons.append(poly)
                    else:
                        final_polygons.append(outer_poly)
                except Exception:
                    final_polygons.append(outer_poly)
            else:
                final_polygons.append(outer_poly)
        
        if final_polygons:
            geom = shapely.union_all(final_polygons)
            if geom.geom_type == "Polygon":
                geom = shapely.geometry.MultiPolygon([geom])
            
            if not geom.is_empty:
                features.append({
                    "geometry": geom,
                    "name": name,
                    "osm_id": rel.get("id"),
                })
    
    if not features:
        return gpd.GeoDataFrame()
    
    return gpd.GeoDataFrame(features, crs="EPSG:4326")


# ---------------------------------------------------------------------------
# 2. Name normalisation
# ---------------------------------------------------------------------------

VIENNA_DISTRICTS_CANON = {
    "innere stadt": "Innere Stadt",
    "leopoldstadt": "Leopoldstadt",
    "landstraße": "Landstraße",
    "landstrasse": "Landstraße",
    "wieden": "Wieden",
    "margareten": "Margareten",
    "margarethen": "Margareten",
    "mariahilf": "Mariahilf",
    "neubau": "Neubau",
    "josefstadt": "Josefstadt",
    "alsergrund": "Alsergrund",
    "favoriten": "Favoriten",
    "simmering": "Simmering",
    "meidling": "Meidling",
    "hietzing": "Hietzing",
    "penzing": "Penzing",
    "rudolfsheim-fünfhaus": "Rudolfsheim-Fünfhaus",
    "rudolfsheim-fuenfhaus": "Rudolfsheim-Fünfhaus",
    "ottakring": "Ottakring",
    "hernals": "Hernals",
    "währing": "Währing",
    "waehring": "Währing",
    "döbling": "Döbling",
    "doebling": "Döbling",
    "brigittenau": "Brigittenau",
    "floridsdorf": "Floridsdorf",
    "donaustadt": "Donaustadt",
    "liesing": "Liesing",
}


def normalise_name(name):
    if not name or pd.isna(name):
        return None
    n = name.strip()
    key = n.lower().replace("  ", " ").replace("ß", "ss").replace("-", "-")
    if key in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key]
    key2 = key.replace("ss", "ß")
    if key2 in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key2]
    for pattern, canon in VIENNA_DISTRICTS_CANON.items():
        if pattern in key or key in pattern:
            return canon
    return n


# ---------------------------------------------------------------------------
# 3. Main pipeline
# ---------------------------------------------------------------------------

def main():
    print("=" * 60)
    print("Fetching CURRENT Vienna district boundaries...")
    print("=" * 60)
    current_gdf = fetch_vienna_districts_osmnx()
    print(f"  Got {len(current_gdf)} current district features")
    
    print("\n" + "=" * 60)
    print("Fetching 2014-01-01 Vienna district boundaries...")
    print("=" * 60)
    historical_gdf = fetch_vienna_districts_osmnx(date="2014-01-01T00:00:00Z")
    print(f"  Got {len(historical_gdf)} historical district features")
    
    if current_gdf.empty or historical_gdf.empty:
        print("ERROR: Empty data returned!")
        sys.exit(1)
    
    # Normalise names
    current_gdf["district_name"] = current_gdf["name"].apply(normalise_name)
    historical_gdf["district_name"] = historical_gdf["name"].apply(normalise_name)
    
    print("\nCurrent district names:")
    for n in sorted(current_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    print("\nHistorical district names:")
    for n in sorted(historical_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    
    # Filter to Vienna districts only
    valid_names = set(VIENNA_DISTRICTS_CANON.values())
    current_gdf = current_gdf[current_gdf["district_name"].isin(valid_names)].copy()
    historical_gdf = historical_gdf[historical_gdf["district_name"].isin(valid_names)].copy()
    
    print(f"\nAfter filtering to Vienna districts:")
    print(f"  Current: {len(current_gdf)} districts")
    print(f"  Historical: {len(historical_gdf)} districts")
    
    # -----------------------------------------------------------------------
    # Strategy: dissolve to city level first, find changes there,
    # then clip back per-district for attribution
    # -----------------------------------------------------------------------
    
    # First, create a clean dissolved union of each snapshot
    current_union = shapely.union_all(current_gdf.geometry)
    historical_union = shapely.union_all(historical_gdf.geometry)
    
    print(f"\nCurrent union area: {current_union.area:.8f}")
    print(f"Historical union area: {historical_union.area:.8f}")
    
    # Snap historical to current to eliminate coordinate precision differences
    TOLERANCE = 1e-7  # ~1cm in degrees
    snapped_historical = shapely.snap(historical_union, current_union, TOLERANCE)
    
    # Compute symmetric difference at city level
    added = shapely.difference(current_union, snapped_historical)
    removed = shapely.difference(snapped_historical, current_union)
    unchanged = shapely.intersection(current_union, snapped_historical)
    
    print(f"\nAdded: area={added.area:.10f}, empty={added.is_empty}")
    print(f"Removed: area={removed.area:.10f}, empty={removed.is_empty}")
    print(f"Unchanged: area={unchanged.area:.10f}, empty={unchanged.is_empty}")
    
    # Remove tiny slivers
    MIN_AREA = 1e-9  # ~ 0.01 sq m in degrees
    
    def filter_slivers(geom):
        if geom.is_empty:
            return shapely.geometry.MultiPolygon()
        if geom.geom_type == "MultiPolygon":
            polys = [g for g in geom.geoms if g.area >= MIN_AREA]
            if not polys:
                return shapely.geometry.MultiPolygon()
            if len(polys) == 1:
                return shapely.geometry.MultiPolygon([polys[0]])
            return shapely.geometry.MultiPolygon(polys)
        elif geom.geom_type == "Polygon":
            if geom.area >= MIN_AREA:
                return shapely.geometry.MultiPolygon([geom])
            return shapely.geometry.MultiPolygon()
        return geom
    
    added = filter_slivers(added)
    removed = filter_slivers(removed)
    unchanged = filter_slivers(unchanged)
    
    print(f"\nAfter sliver removal:")
    print(f"  Added: area={added.area:.10f}" if not added.is_empty else "  Added: empty")
    print(f"  Removed: area={removed.area:.10f}" if not removed.is_empty else "  Removed: empty")
    print(f"  Unchanged: area={unchanged.area:.10f}" if not unchanged.is_empty else "  Unchanged: empty")
    
    # -----------------------------------------------------------------------
    # Build fragments with district attribution
    # -----------------------------------------------------------------------
    fragments = []
    
    def add_fragments(geom, ctype, source_gdf):
        """Add geometry fragments, assigning district_name by intersection with source_gdf."""
        if geom.is_empty:
            return
        if geom.geom_type == "MultiPolygon":
            polys = list(geom.geoms)
        elif geom.geom_type == "Polygon":
            polys = [geom]
        else:
            polys = []
        
        for poly in polys:
            if poly.area < MIN_AREA:
                continue
            poly_mp = shapely.geometry.MultiPolygon([poly]) if poly.geom_type == "Polygon" else poly
            
            # Find which district(s) this belongs to
            districts = set()
            for _, row in source_gdf.iterrows():
                if poly_mp.intersects(row["geometry"]):
                    districts.add(row["district_name"])
            
            dname = ", ".join(sorted(districts)) if districts else "unknown"
            fragments.append({
                "geometry": poly_mp,
                "change_type": ctype,
                "district_name": dname,
            })
    
    # Added areas: from current districts
    add_fragments(added, "added_since_2014", current_gdf)
    # Removed areas: from historical districts
    add_fragments(removed, "removed_since_2014", historical_gdf)
    
    # Unchanged: clip by each current district separately
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        dist_geom = row["geometry"]
        unchanged_part = shapely.intersection(dist_geom, unchanged)
        if not unchanged_part.is_empty:
            unchanged_part = filter_slivers(unchanged_part)
            if not unchanged_part.is_empty:
                if unchanged_part.geom_type == "MultiPolygon":
                    polys = list(unchanged_part.geoms)
                elif unchanged_part.geom_type == "Polygon":
                    polys = [unchanged_part]
                else:
                    polys = []
                
                for poly in polys:
                    if poly.area < MIN_AREA:
                        continue
                    if poly.geom_type == "Polygon":
                        poly_mp = shapely.geometry.MultiPolygon([poly])
                    else:
                        poly_mp = poly
                    fragments.append({
                        "geometry": poly_mp,
                        "change_type": "unchanged",
                        "district_name": dname,
                    })
    
    result_gdf = gpd.GeoDataFrame(fragments, crs="EPSG:4326")
    
    if result_gdf.empty:
        print("No fragments produced. Using current boundaries as unchanged.")
        result_gdf = current_gdf.copy()
        result_gdf["change_type"] = "unchanged"
        result_gdf = result_gdf[["geometry", "change_type", "district_name"]]
    
    print(f"\nFragments created: {len(result_gdf)}")
    if not result_gdf.empty:
        print(result_gdf["change_type"].value_counts())
    
    # -----------------------------------------------------------------------
    # 4. touches_changed
    # -----------------------------------------------------------------------
    changed_geoms = []
    for _, row in result_gdf.iterrows():
        if row["change_type"] in ("added_since_2014", "removed_since_2014"):
            changed_geoms.append(row["geometry"])
    
    if changed_geoms:
        changed_union = shapely.union_all(changed_geoms)
    else:
        changed_union = shapely.geometry.MultiPolygon()
    
    district_touches = {}
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        geom = row["geometry"]
        if not changed_union.is_empty:
            touches = geom.touches(changed_union)
            if not touches:
                overlap = shapely.intersection(geom, changed_union)
                touches = not overlap.is_empty
            district_touches[dname] = touches
        else:
            district_touches[dname] = False
    
    for _, row in historical_gdf.iterrows():
        dname = row["district_name"]
        if dname not in district_touches:
            geom = row["geometry"]
            if not changed_union.is_empty:
                touches = geom.touches(changed_union)
                if not touches:
                    overlap = shapely.intersection(geom, changed_union)
                    touches = not overlap.is_empty
                district_touches[dname] = touches
            else:
                district_touches[dname] = False
    
    def get_touches(dn_str):
        if pd.isna(dn_str) or dn_str == "unknown":
            return False
        for d in dn_str.split(","):
            d = d.strip()
            if district_touches.get(d, False):
                return True
        return False
    
    result_gdf["touches_changed"] = result_gdf["district_name"].apply(get_touches)
    
    # Simplify and ensure MultiPolygon
    result_gdf["geometry"] = result_gdf["geometry"].simplify(0.00001, preserve_topology=True)
    result_gdf["geometry"] = result_gdf["geometry"].buffer(0)
    
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == "Polygon":
            return shapely.geometry.MultiPolygon([geom])
        return geom
    
    result_gdf["geometry"] = result_gdf["geometry"].apply(ensure_multipolygon)
    
    # -----------------------------------------------------------------------
    # 5. Output
    # -----------------------------------------------------------------------
    output_path = "/work/vienna_boundary_changes.geojson"
    result_gdf.to_file(output_path, driver="GeoJSON")
    
    print(f"\n{'=' * 60}")
    print(f"Output written to {output_path}")
    print(f"Features: {len(result_gdf)}")
    print(f"Columns: {result_gdf.columns.tolist()}")
    
    print(f"\nSample features:")
    print(result_gdf[["change_type", "district_name", "touches_changed"]].head(20))
    
    print(f"\nDistrict touch summary:")
    for dname, touches in sorted(district_touches.items()):
        print(f"  {dname}: touches_changed={touches}")
    
    # Verify
    verify = gpd.read_file(output_path)
    print(f"\nVerification: read back {len(verify)} features")
    print(verify["change_type"].value_counts())
    print(f"Geometry types: {verify.geometry.geom_type.unique()}")
    
    # Area stats
    for ct in ["added_since_2014", "removed_since_2014", "unchanged"]:
        sub = verify[verify["change_type"] == ct]
        if not sub.empty:
            areas = sub.geometry.area
            print(f"  {ct}: count={len(sub)}, total_area={areas.sum():.10f}")


if __name__ == "__main__":
    main()

"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass via osmnx,
compare to find changes, and output vienna_boundary_changes.geojson.
"""

import os
import sys
import json
from collections import OrderedDict

import osmnx as ox
import osmnx._overpass as op
import geopandas as gpd
import pandas as pd
import shapely
import shapely.geometry
import numpy as np

ox.settings.log_console = False

# ---------------------------------------------------------------------------
# 1. Overpass query
# ---------------------------------------------------------------------------

def build_ring_from_ways(way_refs, ways_dict, nodes_dict):
    if not way_refs:
        return None, []
    
    segments = []
    for ref in way_refs:
        if ref not in ways_dict:
            continue
        node_ids = ways_dict[ref]
        if len(node_ids) < 2:
            continue
        coords = [nodes_dict.get(nid) for nid in node_ids]
        coords = [c for c in coords if c is not None]
        if len(coords) < 2:
            continue
        segments.append({
            'ref': ref,
            'coords': list(coords),
            'start_node': node_ids[0],
            'end_node': node_ids[-1],
        })
    
    if not segments:
        return None, []
    
    used = [False] * len(segments)
    longest_idx = max(range(len(segments)), key=lambda i: len(segments[i]['coords']))
    ring_coords = list(segments[longest_idx]['coords'])
    current_end_node = segments[longest_idx]['end_node']
    used[longest_idx] = True
    
    changed = True
    while changed:
        changed = False
        for i, seg in enumerate(segments):
            if used[i]:
                continue
            if seg['start_node'] == current_end_node:
                ring_coords.extend(seg['coords'][1:])
                current_end_node = seg['end_node']
                used[i] = True
                changed = True
                break
            elif seg['end_node'] == current_end_node:
                rev_coords = list(reversed(seg['coords']))
                ring_coords.extend(rev_coords[1:])
                current_end_node = seg['start_node']
                used[i] = True
                changed = True
                break
    
    if ring_coords[0] != ring_coords[-1]:
        ring_coords.append(ring_coords[0])
    
    used_refs = [segments[i]['ref'] for i in range(len(segments)) if used[i]]
    return ring_coords, used_refs


def fetch_vienna_districts_osmnx(date=None):
    bbox = (48.1, 16.1, 48.35, 16.6)
    
    if date:
        query = f"""
[out:json][timeout:180][date:"{date}"];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    else:
        query = f"""
[out:json][timeout:180];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    
    data = OrderedDict()
    data['data'] = query
    result = op._overpass_request(data)
    elements = result.get("elements", [])
    
    nodes = {}
    ways = {}
    relations = []
    
    for el in elements:
        etype = el.get("type")
        eid = el.get("id")
        if etype == "node":
            nodes[eid] = (el.get("lon"), el.get("lat"))
        elif etype == "way":
            ways[eid] = el.get("nodes", [])
        elif etype == "relation":
            relations.append(el)
    
    features = []
    
    for rel in relations:
        tags = rel.get("tags", {})
        name = tags.get("name", "")
        admin_level = tags.get("admin_level", "")
        
        if admin_level != "9":
            continue
        
        members = rel.get("members", [])
        
        all_outer_refs = []
        inner_refs_coords = []
        
        for mem in members:
            if mem.get("type") != "way":
                continue
            role = mem.get("role", "")
            ref = mem.get("ref")
            if ref not in ways:
                continue
            if role == "inner":
                node_ids = ways[ref]
                coords = [nodes.get(nid) for nid in node_ids]
                coords = [c for c in coords if c is not None]
                if len(coords) >= 3:
                    if coords[0] != coords[-1]:
                        coords.append(coords[0])
                    inner_refs_coords.append(coords)
            else:
                all_outer_refs.append(ref)
        
        outer_polygons = []
        remaining = list(all_outer_refs)
        
        while remaining:
            ring_coords, used_refs = build_ring_from_ways(remaining, ways, nodes)
            if ring_coords is None or len(ring_coords) < 4:
                break
            
            try:
                ls = shapely.LineString(ring_coords)
                if not ls.is_valid or ls.is_empty or ls.length < 1e-10:
                    break
                poly = shapely.Polygon(ring_coords)
                poly = shapely.orient_polygons(poly)
                if poly.is_valid and not poly.is_empty and poly.area > 1e-12:
                    outer_polygons.append(poly)
            except Exception:
                break
            
            remaining = [r for r in remaining if r not in used_refs]
        
        if not outer_polygons:
            continue
        
        inner_polygons = []
        for coords in inner_refs_coords:
            try:
                ipoly = shapely.Polygon(coords)
                ipoly = shapely.orient_polygons(ipoly)
                if ipoly.is_valid and not ipoly.is_empty:
                    inner_polygons.append(ipoly)
            except Exception:
                continue
        
        final_polygons = []
        for outer_poly in outer_polygons:
            inners_here = []
            rest = []
            for ipoly in inner_polygons:
                try:
                    if outer_poly.contains(ipoly.representative_point()) or outer_poly.intersects(ipoly):
                        inners_here.append(ipoly)
                    else:
                        rest.append(ipoly)
                except Exception:
                    rest.append(ipoly)
            inner_polygons = rest
            
            if inners_here:
                try:
                    outer_coords = list(outer_poly.exterior.coords)
                    inner_coords_list = [list(ip.exterior.coords) for ip in inners_here]
                    poly = shapely.Polygon(outer_coords, inner_coords_list)
                    poly = shapely.orient_polygons(poly)
                    if poly.is_valid and not poly.is_empty:
                        final_polygons.append(poly)
                    else:
                        final_polygons.append(outer_poly)
                except Exception:
                    final_polygons.append(outer_poly)
            else:
                final_polygons.append(outer_poly)
        
        if final_polygons:
            geom = shapely.union_all(final_polygons)
            if geom.geom_type == "Polygon":
                geom = shapely.geometry.MultiPolygon([geom])
            
            if not geom.is_empty:
                features.append({
                    "geometry": geom,
                    "name": name,
                    "osm_id": rel.get("id"),
                })
    
    if not features:
        return gpd.GeoDataFrame()
    
    return gpd.GeoDataFrame(features, crs="EPSG:4326")


# ---------------------------------------------------------------------------
# 2. Name normalisation
# ---------------------------------------------------------------------------

VIENNA_DISTRICTS_CANON = {
    "innere stadt": "Innere Stadt",
    "leopoldstadt": "Leopoldstadt",
    "landstraße": "Landstraße",
    "landstrasse": "Landstraße",
    "wieden": "Wieden",
    "margareten": "Margareten",
    "margarethen": "Margareten",
    "mariahilf": "Mariahilf",
    "neubau": "Neubau",
    "josefstadt": "Josefstadt",
    "alsergrund": "Alsergrund",
    "favoriten": "Favoriten",
    "simmering": "Simmering",
    "meidling": "Meidling",
    "hietzing": "Hietzing",
    "penzing": "Penzing",
    "rudolfsheim-fünfhaus": "Rudolfsheim-Fünfhaus",
    "rudolfsheim-fuenfhaus": "Rudolfsheim-Fünfhaus",
    "ottakring": "Ottakring",
    "hernals": "Hernals",
    "währing": "Währing",
    "waehring": "Währing",
    "döbling": "Döbling",
    "doebling": "Döbling",
    "brigittenau": "Brigittenau",
    "floridsdorf": "Floridsdorf",
    "donaustadt": "Donaustadt",
    "liesing": "Liesing",
}


def normalise_name(name):
    if not name or pd.isna(name):
        return None
    n = name.strip()
    key = n.lower().replace("  ", " ").replace("ß", "ss").replace("-", "-")
    if key in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key]
    key2 = key.replace("ss", "ß")
    if key2 in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key2]
    for pattern, canon in VIENNA_DISTRICTS_CANON.items():
        if pattern in key or key in pattern:
            return canon
    return n


# ---------------------------------------------------------------------------
# 3. Main pipeline
# ---------------------------------------------------------------------------

def main():
    print("=" * 60)
    print("Fetching CURRENT Vienna district boundaries...")
    print("=" * 60)
    current_gdf = fetch_vienna_districts_osmnx()
    print(f"  Got {len(current_gdf)} current district features")
    
    print("\n" + "=" * 60)
    print("Fetching 2014-01-01 Vienna district boundaries...")
    print("=" * 60)
    historical_gdf = fetch_vienna_districts_osmnx(date="2014-01-01T00:00:00Z")
    print(f"  Got {len(historical_gdf)} historical district features")
    
    if current_gdf.empty or historical_gdf.empty:
        print("ERROR: Empty data returned!")
        sys.exit(1)
    
    # Normalise names
    current_gdf["district_name"] = current_gdf["name"].apply(normalise_name)
    historical_gdf["district_name"] = historical_gdf["name"].apply(normalise_name)
    
    print("\nCurrent district names:")
    for n in sorted(current_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    print("\nHistorical district names:")
    for n in sorted(historical_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    
    # Filter to Vienna districts only
    valid_names = set(VIENNA_DISTRICTS_CANON.values())
    current_gdf = current_gdf[current_gdf["district_name"].isin(valid_names)].copy()
    historical_gdf = historical_gdf[historical_gdf["district_name"].isin(valid_names)].copy()
    
    print(f"\nAfter filtering to Vienna districts:")
    print(f"  Current: {len(current_gdf)} districts")
    print(f"  Historical: {len(historical_gdf)} districts")
    
    # Create per-district snapshots indexed by name
    hist_by_name = {row["district_name"]: row["geometry"] for _, row in historical_gdf.iterrows()}
    
    # -----------------------------------------------------------------------
    # Strategy: compare district-per-district, with coordinate snapping
    # -----------------------------------------------------------------------
    
    TOLERANCE = 1e-6  # ~10cm in degrees
    MIN_AREA = 1e-9   # ~0.01 sq m in degrees
    
    fragments = []
    
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        current_geom = row["geometry"]
        historical_geom = hist_by_name.get(dname)
        
        if historical_geom is None:
            # District exists now but not in 2014
            fragments.append({
                "geometry": current_geom,
                "change_type": "added_since_2014",
                "district_name": dname,
            })
            continue
        
        # Snap geometries to each other to eliminate coordinate noise
        current_norm = shapely.normalize(current_geom)
        hist_norm = shapely.normalize(historical_geom)
        
        # Two-way snapping for best alignment
        snapped_current = shapely.snap(current_norm, hist_norm, TOLERANCE)
        snapped_hist = shapely.snap(hist_norm, snapped_current, TOLERANCE)
        
        # Symmetric difference
        added_part = shapely.difference(snapped_current, snapped_hist)
        removed_part = shapely.difference(snapped_hist, snapped_current)
        unchanged_part = shapely.intersection(snapped_current, snapped_hist)
        
        # Clean slivers
        def clean_and_split(geom, label):
            """Clean geometry, filter slivers, return list of (geom, label) fragments."""
            if geom.is_empty:
                return []
            
            if geom.geom_type == "MultiPolygon":
                polys = list(geom.geoms)
            elif geom.geom_type == "Polygon":
                polys = [geom]
            else:
                return []
            
            result = []
            for poly in polys:
                if poly.is_empty or poly.area < MIN_AREA:
                    continue
                if poly.geom_type == "Polygon":
                    mp = shapely.geometry.MultiPolygon([poly])
                else:
                    mp = poly
                result.append((mp, label))
            return result
        
        for mp, ctype in clean_and_split(added_part, "added_since_2014"):
            fragments.append({
                "geometry": mp,
                "change_type": ctype,
                "district_name": dname,
            })
        
        for mp, ctype in clean_and_split(removed_part, "removed_since_2014"):
            fragments.append({
                "geometry": mp,
                "change_type": ctype,
                "district_name": dname,
            })
        
        for mp, ctype in clean_and_split(unchanged_part, "unchanged"):
            fragments.append({
                "geometry": mp,
                "change_type": ctype,
                "district_name": dname,
            })
    
    # Also check for historical districts that no longer exist
    current_names = set(current_gdf["district_name"])
    for _, row in historical_gdf.iterrows():
        dname = row["district_name"]
        if dname not in current_names:
            # District existed in 2014 but not now — fully removed
            fragments.append({
                "geometry": row["geometry"],
                "change_type": "removed_since_2014",
                "district_name": dname,
            })
    
    result_gdf = gpd.GeoDataFrame(fragments, crs="EPSG:4326")
    
    if result_gdf.empty:
        print("No fragments produced. Using current boundaries as unchanged.")
        result_gdf = current_gdf.copy()
        result_gdf["change_type"] = "unchanged"
        result_gdf = result_gdf[["geometry", "change_type", "district_name"]]
    
    print(f"\nFragments created: {len(result_gdf)}")
    if not result_gdf.empty:
        print(result_gdf["change_type"].value_counts())
    
    # -----------------------------------------------------------------------
    # 4. touches_changed
    # -----------------------------------------------------------------------
    changed_geoms = []
    for _, row in result_gdf.iterrows():
        if row["change_type"] in ("added_since_2014", "removed_since_2014"):
            changed_geoms.append(row["geometry"])
    
    if changed_geoms:
        changed_union = shapely.union_all(changed_geoms)
    else:
        changed_union = shapely.geometry.MultiPolygon()
    
    district_touches = {}
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        geom = row["geometry"]
        if not changed_union.is_empty:
            touches = geom.touches(changed_union)
            if not touches:
                overlap = shapely.intersection(geom, changed_union)
                touches = not overlap.is_empty
            district_touches[dname] = touches
        else:
            district_touches[dname] = False
    
    for _, row in historical_gdf.iterrows():
        dname = row["district_name"]
        if dname not in district_touches:
            geom = row["geometry"]
            if not changed_union.is_empty:
                touches = geom.touches(changed_union)
                if not touches:
                    overlap = shapely.intersection(geom, changed_union)
                    touches = not overlap.is_empty
                district_touches[dname] = touches
            else:
                district_touches[dname] = False
    
    def get_touches(dn_str):
        if pd.isna(dn_str) or dn_str == "unknown":
            return False
        for d in dn_str.split(","):
            d = d.strip()
            if district_touches.get(d, False):
                return True
        return False
    
    result_gdf["touches_changed"] = result_gdf["district_name"].apply(get_touches)
    
    # Simplify and ensure MultiPolygon
    result_gdf["geometry"] = result_gdf["geometry"].simplify(0.00001, preserve_topology=True)
    result_gdf["geometry"] = result_gdf["geometry"].buffer(0)
    
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == "Polygon":
            return shapely.geometry.MultiPolygon([geom])
        return geom
    
    result_gdf["geometry"] = result_gdf["geometry"].apply(ensure_multipolygon)
    
    # -----------------------------------------------------------------------
    # 5. Output
    # -----------------------------------------------------------------------
    output_path = "/work/vienna_boundary_changes.geojson"
    result_gdf.to_file(output_path, driver="GeoJSON")
    
    print(f"\n{'=' * 60}")
    print(f"Output written to {output_path}")
    print(f"Features: {len(result_gdf)}")
    print(f"Columns: {result_gdf.columns.tolist()}")
    
    print(f"\nSample features:")
    print(result_gdf[["change_type", "district_name", "touches_changed"]].head(20))
    
    print(f"\nDistrict touch summary:")
    for dname, touches in sorted(district_touches.items()):
        print(f"  {dname}: touches_changed={touches}")
    
    # Verify
    verify = gpd.read_file(output_path)
    print(f"\nVerification: read back {len(verify)} features")
    print(verify["change_type"].value_counts())
    print(f"Geometry types: {verify.geometry.geom_type.unique()}")
    
    # Area stats
    for ct in ["added_since_2014", "removed_since_2014", "unchanged"]:
        sub = verify[verify["change_type"] == ct]
        if not sub.empty:
            areas = sub.geometry.area
            print(f"  {ct}: count={len(sub)}, total_area={areas.sum():.10f}")


if __name__ == "__main__":
    main()

"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass via osmnx,
compare to find changes, and output vienna_boundary_changes.geojson.

Strategy: dissolve to city-level boundary first, compute symmetric difference,
then attribute fragments back to districts. Use generous tolerance to avoid
slivers from coordinate misalignment.
"""

import os
import sys
import json
from collections import OrderedDict

import osmnx as ox
import osmnx._overpass as op
import geopandas as gpd
import pandas as pd
import shapely
import shapely.geometry
import numpy as np

ox.settings.log_console = False

# ---------------------------------------------------------------------------
# 1. Overpass query
# ---------------------------------------------------------------------------

def build_ring_from_ways(way_refs, ways_dict, nodes_dict):
    if not way_refs:
        return None, []
    
    segments = []
    for ref in way_refs:
        if ref not in ways_dict:
            continue
        node_ids = ways_dict[ref]
        if len(node_ids) < 2:
            continue
        coords = [nodes_dict.get(nid) for nid in node_ids]
        coords = [c for c in coords if c is not None]
        if len(coords) < 2:
            continue
        segments.append({
            'ref': ref,
            'coords': list(coords),
            'start_node': node_ids[0],
            'end_node': node_ids[-1],
        })
    
    if not segments:
        return None, []
    
    used = [False] * len(segments)
    longest_idx = max(range(len(segments)), key=lambda i: len(segments[i]['coords']))
    ring_coords = list(segments[longest_idx]['coords'])
    current_end_node = segments[longest_idx]['end_node']
    used[longest_idx] = True
    
    changed = True
    while changed:
        changed = False
        for i, seg in enumerate(segments):
            if used[i]:
                continue
            if seg['start_node'] == current_end_node:
                ring_coords.extend(seg['coords'][1:])
                current_end_node = seg['end_node']
                used[i] = True
                changed = True
                break
            elif seg['end_node'] == current_end_node:
                rev_coords = list(reversed(seg['coords']))
                ring_coords.extend(rev_coords[1:])
                current_end_node = seg['start_node']
                used[i] = True
                changed = True
                break
    
    if ring_coords[0] != ring_coords[-1]:
        ring_coords.append(ring_coords[0])
    
    used_refs = [segments[i]['ref'] for i in range(len(segments)) if used[i]]
    return ring_coords, used_refs


def fetch_vienna_districts_osmnx(date=None):
    bbox = (48.1, 16.1, 48.35, 16.6)
    
    if date:
        query = f"""
[out:json][timeout:180][date:"{date}"];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    else:
        query = f"""
[out:json][timeout:180];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    
    data = OrderedDict()
    data['data'] = query
    result = op._overpass_request(data)
    elements = result.get("elements", [])
    
    nodes = {}
    ways = {}
    relations = []
    
    for el in elements:
        etype = el.get("type")
        eid = el.get("id")
        if etype == "node":
            nodes[eid] = (el.get("lon"), el.get("lat"))
        elif etype == "way":
            ways[eid] = el.get("nodes", [])
        elif etype == "relation":
            relations.append(el)
    
    features = []
    
    for rel in relations:
        tags = rel.get("tags", {})
        name = tags.get("name", "")
        admin_level = tags.get("admin_level", "")
        
        if admin_level != "9":
            continue
        
        members = rel.get("members", [])
        
        all_outer_refs = []
        inner_refs_coords = []
        
        for mem in members:
            if mem.get("type") != "way":
                continue
            role = mem.get("role", "")
            ref = mem.get("ref")
            if ref not in ways:
                continue
            if role == "inner":
                node_ids = ways[ref]
                coords = [nodes.get(nid) for nid in node_ids]
                coords = [c for c in coords if c is not None]
                if len(coords) >= 3:
                    if coords[0] != coords[-1]:
                        coords.append(coords[0])
                    inner_refs_coords.append(coords)
            else:
                all_outer_refs.append(ref)
        
        outer_polygons = []
        remaining = list(all_outer_refs)
        
        while remaining:
            ring_coords, used_refs = build_ring_from_ways(remaining, ways, nodes)
            if ring_coords is None or len(ring_coords) < 4:
                break
            
            try:
                ls = shapely.LineString(ring_coords)
                if not ls.is_valid or ls.is_empty or ls.length < 1e-10:
                    break
                poly = shapely.Polygon(ring_coords)
                poly = shapely.orient_polygons(poly)
                if poly.is_valid and not poly.is_empty and poly.area > 1e-12:
                    outer_polygons.append(poly)
            except Exception:
                break
            
            remaining = [r for r in remaining if r not in used_refs]
        
        if not outer_polygons:
            continue
        
        inner_polygons = []
        for coords in inner_refs_coords:
            try:
                ipoly = shapely.Polygon(coords)
                ipoly = shapely.orient_polygons(ipoly)
                if ipoly.is_valid and not ipoly.is_empty:
                    inner_polygons.append(ipoly)
            except Exception:
                continue
        
        final_polygons = []
        for outer_poly in outer_polygons:
            inners_here = []
            rest = []
            for ipoly in inner_polygons:
                try:
                    if outer_poly.contains(ipoly.representative_point()) or outer_poly.intersects(ipoly):
                        inners_here.append(ipoly)
                    else:
                        rest.append(ipoly)
                except Exception:
                    rest.append(ipoly)
            inner_polygons = rest
            
            if inners_here:
                try:
                    outer_coords = list(outer_poly.exterior.coords)
                    inner_coords_list = [list(ip.exterior.coords) for ip in inners_here]
                    poly = shapely.Polygon(outer_coords, inner_coords_list)
                    poly = shapely.orient_polygons(poly)
                    if poly.is_valid and not poly.is_empty:
                        final_polygons.append(poly)
                    else:
                        final_polygons.append(outer_poly)
                except Exception:
                    final_polygons.append(outer_poly)
            else:
                final_polygons.append(outer_poly)
        
        if final_polygons:
            geom = shapely.union_all(final_polygons)
            if geom.geom_type == "Polygon":
                geom = shapely.geometry.MultiPolygon([geom])
            
            if not geom.is_empty:
                features.append({
                    "geometry": geom,
                    "name": name,
                    "osm_id": rel.get("id"),
                })
    
    if not features:
        return gpd.GeoDataFrame()
    
    return gpd.GeoDataFrame(features, crs="EPSG:4326")


# ---------------------------------------------------------------------------
# 2. Name normalisation
# ---------------------------------------------------------------------------

VIENNA_DISTRICTS_CANON = {
    "innere stadt": "Innere Stadt",
    "leopoldstadt": "Leopoldstadt",
    "landstraße": "Landstraße",
    "landstrasse": "Landstraße",
    "wieden": "Wieden",
    "margareten": "Margareten",
    "margarethen": "Margareten",
    "mariahilf": "Mariahilf",
    "neubau": "Neubau",
    "josefstadt": "Josefstadt",
    "alsergrund": "Alsergrund",
    "favoriten": "Favoriten",
    "simmering": "Simmering",
    "meidling": "Meidling",
    "hietzing": "Hietzing",
    "penzing": "Penzing",
    "rudolfsheim-fünfhaus": "Rudolfsheim-Fünfhaus",
    "rudolfsheim-fuenfhaus": "Rudolfsheim-Fünfhaus",
    "ottakring": "Ottakring",
    "hernals": "Hernals",
    "währing": "Währing",
    "waehring": "Währing",
    "döbling": "Döbling",
    "doebling": "Döbling",
    "brigittenau": "Brigittenau",
    "floridsdorf": "Floridsdorf",
    "donaustadt": "Donaustadt",
    "liesing": "Liesing",
}


def normalise_name(name):
    if not name or pd.isna(name):
        return None
    n = name.strip()
    key = n.lower().replace("  ", " ").replace("ß", "ss").replace("-", "-")
    if key in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key]
    key2 = key.replace("ss", "ß")
    if key2 in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key2]
    for pattern, canon in VIENNA_DISTRICTS_CANON.items():
        if pattern in key or key in pattern:
            return canon
    return n


# ---------------------------------------------------------------------------
# 3. Main pipeline
# ---------------------------------------------------------------------------

def main():
    print("=" * 60)
    print("Fetching CURRENT Vienna district boundaries...")
    print("=" * 60)
    current_gdf = fetch_vienna_districts_osmnx()
    print(f"  Got {len(current_gdf)} current district features")
    
    print("\n" + "=" * 60)
    print("Fetching 2014-01-01 Vienna district boundaries...")
    print("=" * 60)
    historical_gdf = fetch_vienna_districts_osmnx(date="2014-01-01T00:00:00Z")
    print(f"  Got {len(historical_gdf)} historical district features")
    
    if current_gdf.empty or historical_gdf.empty:
        print("ERROR: Empty data returned!")
        sys.exit(1)
    
    # Normalise names
    current_gdf["district_name"] = current_gdf["name"].apply(normalise_name)
    historical_gdf["district_name"] = historical_gdf["name"].apply(normalise_name)
    
    print("\nCurrent district names:")
    for n in sorted(current_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    print("\nHistorical district names:")
    for n in sorted(historical_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    
    # Filter to Vienna districts only
    valid_names = set(VIENNA_DISTRICTS_CANON.values())
    current_gdf = current_gdf[current_gdf["district_name"].isin(valid_names)].copy()
    historical_gdf = historical_gdf[historical_gdf["district_name"].isin(valid_names)].copy()
    
    print(f"\nAfter filtering to Vienna districts:")
    print(f"  Current: {len(current_gdf)} districts")
    print(f"  Historical: {len(historical_gdf)} districts")
    
    # -----------------------------------------------------------------------
    # Strategy: dissolve, snap city boundary, diff, then attribute back
    # -----------------------------------------------------------------------
    
    # 1. Dissolve each snapshot to city-level boundary
    current_union = shapely.union_all(current_gdf.geometry)
    historical_union = shapely.union_all(historical_gdf.geometry)
    
    print(f"\nCurrent union area: {current_union.area:.10f}")
    print(f"Historical union area: {historical_union.area:.10f}")
    
    # 2. Snap historical to current to eliminate coordinate noise
    TOLERANCE = 1e-6  # ~10 cm in degrees
    snapped_hist = shapely.snap(historical_union, current_union, TOLERANCE)
    
    # 3. Compute symmetric difference
    added = shapely.difference(current_union, snapped_hist)
    removed = shapely.difference(snapped_hist, current_union)
    unchanged = shapely.intersection(current_union, snapped_hist)
    
    # 4. Remove slivers: keep only polygons above area threshold
    MIN_AREA = 1e-9  # ~0.01 sq m in degrees at Vienna latitude
    
    def remove_slivers(geom):
        if geom.is_empty:
            return shapely.geometry.MultiPolygon()
        if geom.geom_type == "MultiPolygon":
            polys = [g for g in geom.geoms if g.area >= MIN_AREA]
            if not polys:
                return shapely.geometry.MultiPolygon()
            if len(polys) == 1:
                return shapely.geometry.MultiPolygon([polys[0]])
            return shapely.geometry.MultiPolygon(polys)
        elif geom.geom_type == "Polygon":
            if geom.area >= MIN_AREA:
                return shapely.geometry.MultiPolygon([geom])
            return shapely.geometry.MultiPolygon()
        return geom
    
    added = remove_slivers(added)
    removed = remove_slivers(removed)
    unchanged = remove_slivers(unchanged)
    
    print(f"\nAfter cleaning:")
    print(f"  Added: area={added.area:.10f}, empty={added.is_empty}")
    print(f"  Removed: area={removed.area:.10f}, empty={removed.is_empty}")
    print(f"  Unchanged: area={unchanged.area:.10f}, empty={unchanged.is_empty}")
    
    # 5. Attribute fragments to districts
    fragments = []
    
    def add_fragments(geom, ctype, source_gdf):
        """Add geometry fragments, attributing district_name by intersection with source_gdf."""
        if geom.is_empty:
            return
        if geom.geom_type == "MultiPolygon":
            polys = list(geom.geoms)
        elif geom.geom_type == "Polygon":
            polys = [geom]
        else:
            return
        
        for poly in polys:
            if poly.area < MIN_AREA:
                continue
            if poly.geom_type == "Polygon":
                mp = shapely.geometry.MultiPolygon([poly])
            elif poly.geom_type == "MultiPolygon":
                mp = poly
            else:
                continue
            
            # Find which district(s) this fragment intersects
            districts = set()
            for _, row in source_gdf.iterrows():
                if mp.intersects(row["geometry"]):
                    districts.add(row["district_name"])
            
            dname = ", ".join(sorted(districts)) if districts else "unknown"
            fragments.append({
                "geometry": mp,
                "change_type": ctype,
                "district_name": dname,
            })
    
    # Added areas attributed from current districts
    add_fragments(added, "added_since_2014", current_gdf)
    # Removed areas attributed from historical districts
    add_fragments(removed, "removed_since_2014", historical_gdf)
    
    # Unchanged: clip unchanged area by each current district
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        dist_geom = row["geometry"]
        unchanged_part = shapely.intersection(dist_geom, unchanged)
        if not unchanged_part.is_empty:
            unchanged_part = remove_slivers(unchanged_part)
            if not unchanged_part.is_empty:
                if unchanged_part.geom_type == "MultiPolygon":
                    polys = list(unchanged_part.geoms)
                elif unchanged_part.geom_type == "Polygon":
                    polys = [unchanged_part]
                else:
                    polys = []
                
                for poly in polys:
                    if poly.area < MIN_AREA:
                        continue
                    if poly.geom_type == "Polygon":
                        mp = shapely.geometry.MultiPolygon([poly])
                    elif poly.geom_type == "MultiPolygon":
                        mp = poly
                    else:
                        continue
                    fragments.append({
                        "geometry": mp,
                        "change_type": "unchanged",
                        "district_name": dname,
                    })
    
    result_gdf = gpd.GeoDataFrame(fragments, crs="EPSG:4326")
    
    if result_gdf.empty:
        print("No fragments produced. Using current boundaries as unchanged.")
        result_gdf = current_gdf.copy()
        result_gdf["change_type"] = "unchanged"
        result_gdf = result_gdf[["geometry", "change_type", "district_name"]]
    
    print(f"\nFragments created: {len(result_gdf)}")
    if not result_gdf.empty:
        print(result_gdf["change_type"].value_counts())
    
    # -----------------------------------------------------------------------
    # 4. touches_changed
    # -----------------------------------------------------------------------
    changed_geoms = []
    for _, row in result_gdf.iterrows():
        if row["change_type"] in ("added_since_2014", "removed_since_2014"):
            changed_geoms.append(row["geometry"])
    
    if changed_geoms:
        changed_union = shapely.union_all(changed_geoms)
    else:
        changed_union = shapely.geometry.MultiPolygon()
    
    district_touches = {}
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        geom = row["geometry"]
        if not changed_union.is_empty:
            # Check if the district touches OR intersects changed areas
            touches = geom.touches(changed_union)
            if not touches:
                overlap = shapely.intersection(geom, changed_union)
                touches = not overlap.is_empty
            district_touches[dname] = touches
        else:
            district_touches[dname] = False
    
    for _, row in historical_gdf.iterrows():
        dname = row["district_name"]
        if dname not in district_touches:
            geom = row["geometry"]
            if not changed_union.is_empty:
                touches = geom.touches(changed_union)
                if not touches:
                    overlap = shapely.intersection(geom, changed_union)
                    touches = not overlap.is_empty
                district_touches[dname] = touches
            else:
                district_touches[dname] = False
    
    def get_touches(dn_str):
        if pd.isna(dn_str) or dn_str == "unknown":
            return False
        for d in dn_str.split(","):
            d = d.strip()
            if district_touches.get(d, False):
                return True
        return False
    
    result_gdf["touches_changed"] = result_gdf["district_name"].apply(get_touches)
    
    # Simplify and ensure MultiPolygon
    result_gdf["geometry"] = result_gdf["geometry"].simplify(0.00001, preserve_topology=True)
    result_gdf["geometry"] = result_gdf["geometry"].buffer(0)
    
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == "Polygon":
            return shapely.geometry.MultiPolygon([geom])
        return geom
    
    result_gdf["geometry"] = result_gdf["geometry"].apply(ensure_multipolygon)
    
    # -----------------------------------------------------------------------
    # 5. Output
    # -----------------------------------------------------------------------
    output_path = "/work/vienna_boundary_changes.geojson"
    result_gdf.to_file(output_path, driver="GeoJSON")
    
    print(f"\n{'=' * 60}")
    print(f"Output written to {output_path}")
    print(f"Features: {len(result_gdf)}")
    print(f"Columns: {result_gdf.columns.tolist()}")
    
    print(f"\nSample features:")
    print(result_gdf[["change_type", "district_name", "touches_changed"]].head(20))
    
    print(f"\nDistrict touch summary:")
    for dname, touches in sorted(district_touches.items()):
        print(f"  {dname}: touches_changed={touches}")
    
    # Verify
    verify = gpd.read_file(output_path)
    print(f"\nVerification: read back {len(verify)} features")
    print(verify["change_type"].value_counts())
    print(f"Geometry types: {verify.geometry.geom_type.unique()}")
    
    # Area stats
    for ct in ["added_since_2014", "removed_since_2014", "unchanged"]:
        sub = verify[verify["change_type"] == ct]
        if not sub.empty:
            areas = sub.geometry.area
            print(f"  {ct}: count={len(sub)}, total_area={areas.sum():.10f}")


if __name__ == "__main__":
    main()

"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass via osmnx,
compare to find changes, and output vienna_boundary_changes.geojson.
"""

import os
import sys
import json
from collections import OrderedDict

import osmnx as ox
import osmnx._overpass as op
import geopandas as gpd
import pandas as pd
import shapely
import shapely.geometry
import numpy as np

ox.settings.log_console = False

# ---------------------------------------------------------------------------
# 1. Overpass query
# ---------------------------------------------------------------------------

def build_ring_from_ways(way_refs, ways_dict, nodes_dict):
    if not way_refs:
        return None, []
    
    segments = []
    for ref in way_refs:
        if ref not in ways_dict:
            continue
        node_ids = ways_dict[ref]
        if len(node_ids) < 2:
            continue
        coords = [nodes_dict.get(nid) for nid in node_ids]
        coords = [c for c in coords if c is not None]
        if len(coords) < 2:
            continue
        segments.append({
            'ref': ref,
            'coords': list(coords),
            'start_node': node_ids[0],
            'end_node': node_ids[-1],
        })
    
    if not segments:
        return None, []
    
    used = [False] * len(segments)
    longest_idx = max(range(len(segments)), key=lambda i: len(segments[i]['coords']))
    ring_coords = list(segments[longest_idx]['coords'])
    current_end_node = segments[longest_idx]['end_node']
    used[longest_idx] = True
    
    changed = True
    while changed:
        changed = False
        for i, seg in enumerate(segments):
            if used[i]:
                continue
            if seg['start_node'] == current_end_node:
                ring_coords.extend(seg['coords'][1:])
                current_end_node = seg['end_node']
                used[i] = True
                changed = True
                break
            elif seg['end_node'] == current_end_node:
                rev_coords = list(reversed(seg['coords']))
                ring_coords.extend(rev_coords[1:])
                current_end_node = seg['start_node']
                used[i] = True
                changed = True
                break
    
    if ring_coords[0] != ring_coords[-1]:
        ring_coords.append(ring_coords[0])
    
    used_refs = [segments[i]['ref'] for i in range(len(segments)) if used[i]]
    return ring_coords, used_refs


def fetch_vienna_districts_osmnx(date=None):
    bbox = (48.1, 16.1, 48.35, 16.6)
    
    if date:
        query = f"""
[out:json][timeout:180][date:"{date}"];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    else:
        query = f"""
[out:json][timeout:180];
relation["admin_level"="9"]["boundary"="administrative"]({bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]});
out body;
>;
out skel qt;
"""
    
    data = OrderedDict()
    data['data'] = query
    result = op._overpass_request(data)
    elements = result.get("elements", [])
    
    nodes = {}
    ways = {}
    relations = []
    
    for el in elements:
        etype = el.get("type")
        eid = el.get("id")
        if etype == "node":
            nodes[eid] = (el.get("lon"), el.get("lat"))
        elif etype == "way":
            ways[eid] = el.get("nodes", [])
        elif etype == "relation":
            relations.append(el)
    
    features = []
    
    for rel in relations:
        tags = rel.get("tags", {})
        name = tags.get("name", "")
        admin_level = tags.get("admin_level", "")
        
        if admin_level != "9":
            continue
        
        members = rel.get("members", [])
        
        all_outer_refs = []
        inner_refs_coords = []
        
        for mem in members:
            if mem.get("type") != "way":
                continue
            role = mem.get("role", "")
            ref = mem.get("ref")
            if ref not in ways:
                continue
            if role == "inner":
                node_ids = ways[ref]
                coords = [nodes.get(nid) for nid in node_ids]
                coords = [c for c in coords if c is not None]
                if len(coords) >= 3:
                    if coords[0] != coords[-1]:
                        coords.append(coords[0])
                    inner_refs_coords.append(coords)
            else:
                all_outer_refs.append(ref)
        
        outer_polygons = []
        remaining = list(all_outer_refs)
        
        while remaining:
            ring_coords, used_refs = build_ring_from_ways(remaining, ways, nodes)
            if ring_coords is None or len(ring_coords) < 4:
                break
            
            try:
                ls = shapely.LineString(ring_coords)
                if not ls.is_valid or ls.is_empty or ls.length < 1e-10:
                    break
                poly = shapely.Polygon(ring_coords)
                poly = shapely.orient_polygons(poly)
                if poly.is_valid and not poly.is_empty and poly.area > 1e-12:
                    outer_polygons.append(poly)
            except Exception:
                break
            
            remaining = [r for r in remaining if r not in used_refs]
        
        if not outer_polygons:
            continue
        
        inner_polygons = []
        for coords in inner_refs_coords:
            try:
                ipoly = shapely.Polygon(coords)
                ipoly = shapely.orient_polygons(ipoly)
                if ipoly.is_valid and not ipoly.is_empty:
                    inner_polygons.append(ipoly)
            except Exception:
                continue
        
        final_polygons = []
        for outer_poly in outer_polygons:
            inners_here = []
            rest = []
            for ipoly in inner_polygons:
                try:
                    if outer_poly.contains(ipoly.representative_point()) or outer_poly.intersects(ipoly):
                        inners_here.append(ipoly)
                    else:
                        rest.append(ipoly)
                except Exception:
                    rest.append(ipoly)
            inner_polygons = rest
            
            if inners_here:
                try:
                    outer_coords = list(outer_poly.exterior.coords)
                    inner_coords_list = [list(ip.exterior.coords) for ip in inners_here]
                    poly = shapely.Polygon(outer_coords, inner_coords_list)
                    poly = shapely.orient_polygons(poly)
                    if poly.is_valid and not poly.is_empty:
                        final_polygons.append(poly)
                    else:
                        final_polygons.append(outer_poly)
                except Exception:
                    final_polygons.append(outer_poly)
            else:
                final_polygons.append(outer_poly)
        
        if final_polygons:
            geom = shapely.union_all(final_polygons)
            if geom.geom_type == "Polygon":
                geom = shapely.geometry.MultiPolygon([geom])
            
            if not geom.is_empty:
                features.append({
                    "geometry": geom,
                    "name": name,
                    "osm_id": rel.get("id"),
                })
    
    if not features:
        return gpd.GeoDataFrame()
    
    return gpd.GeoDataFrame(features, crs="EPSG:4326")


# ---------------------------------------------------------------------------
# 2. Name normalisation
# ---------------------------------------------------------------------------

VIENNA_DISTRICTS_CANON = {
    "innere stadt": "Innere Stadt",
    "leopoldstadt": "Leopoldstadt",
    "landstraße": "Landstraße",
    "landstrasse": "Landstraße",
    "wieden": "Wieden",
    "margareten": "Margareten",
    "margarethen": "Margareten",
    "mariahilf": "Mariahilf",
    "neubau": "Neubau",
    "josefstadt": "Josefstadt",
    "alsergrund": "Alsergrund",
    "favoriten": "Favoriten",
    "simmering": "Simmering",
    "meidling": "Meidling",
    "hietzing": "Hietzing",
    "penzing": "Penzing",
    "rudolfsheim-fünfhaus": "Rudolfsheim-Fünfhaus",
    "rudolfsheim-fuenfhaus": "Rudolfsheim-Fünfhaus",
    "ottakring": "Ottakring",
    "hernals": "Hernals",
    "währing": "Währing",
    "waehring": "Währing",
    "döbling": "Döbling",
    "doebling": "Döbling",
    "brigittenau": "Brigittenau",
    "floridsdorf": "Floridsdorf",
    "donaustadt": "Donaustadt",
    "liesing": "Liesing",
}


def normalise_name(name):
    if not name or pd.isna(name):
        return None
    n = name.strip()
    key = n.lower().replace("  ", " ").replace("ß", "ss").replace("-", "-")
    if key in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key]
    key2 = key.replace("ss", "ß")
    if key2 in VIENNA_DISTRICTS_CANON:
        return VIENNA_DISTRICTS_CANON[key2]
    for pattern, canon in VIENNA_DISTRICTS_CANON.items():
        if pattern in key or key in pattern:
            return canon
    return n


# ---------------------------------------------------------------------------
# 3. Main pipeline
# ---------------------------------------------------------------------------

def main():
    print("=" * 60)
    print("Fetching CURRENT Vienna district boundaries...")
    print("=" * 60)
    current_gdf = fetch_vienna_districts_osmnx()
    print(f"  Got {len(current_gdf)} current district features")
    
    print("\n" + "=" * 60)
    print("Fetching 2014-01-01 Vienna district boundaries...")
    print("=" * 60)
    historical_gdf = fetch_vienna_districts_osmnx(date="2014-01-01T00:00:00Z")
    print(f"  Got {len(historical_gdf)} historical district features")
    
    if current_gdf.empty or historical_gdf.empty:
        print("ERROR: Empty data returned!")
        sys.exit(1)
    
    # Normalise names
    current_gdf["district_name"] = current_gdf["name"].apply(normalise_name)
    historical_gdf["district_name"] = historical_gdf["name"].apply(normalise_name)
    
    print("\nCurrent district names:")
    for n in sorted(current_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    print("\nHistorical district names:")
    for n in sorted(historical_gdf["district_name"].dropna().unique()):
        print(f"  - {n}")
    
    # Filter to Vienna districts only
    valid_names = set(VIENNA_DISTRICTS_CANON.values())
    current_gdf = current_gdf[current_gdf["district_name"].isin(valid_names)].copy()
    historical_gdf = historical_gdf[historical_gdf["district_name"].isin(valid_names)].copy()
    
    print(f"\nAfter filtering to Vienna districts:")
    print(f"  Current: {len(current_gdf)} districts")
    print(f"  Historical: {len(historical_gdf)} districts")
    
    # -----------------------------------------------------------------------
    # Strategy: dissolve city boundary, diff, attribute, merge per district
    # -----------------------------------------------------------------------
    
    # 1. Dissolve to city-level
    current_union = shapely.union_all(current_gdf.geometry)
    historical_union = shapely.union_all(historical_gdf.geometry)
    
    print(f"\nCurrent union area: {current_union.area:.10f}")
    print(f"Historical union area: {historical_union.area:.10f}")
    
    # 2. Snap & compute symmetric difference  
    TOLERANCE = 1e-6
    snapped_hist = shapely.snap(historical_union, current_union, TOLERANCE)
    
    added = shapely.difference(current_union, snapped_hist)
    removed = shapely.difference(snapped_hist, current_union)
    unchanged = shapely.intersection(current_union, snapped_hist)
    
    # 3. Remove slivers
    MIN_AREA = 5e-9  # ~40 sq m in reality at Vienna's latitude
    
    def remove_slivers(geom):
        if geom.is_empty:
            return shapely.geometry.MultiPolygon()
        if geom.geom_type == "MultiPolygon":
            polys = [g for g in geom.geoms if g.area >= MIN_AREA]
            if not polys:
                return shapely.geometry.MultiPolygon()
            return shapely.geometry.MultiPolygon(polys) if len(polys) > 1 else shapely.geometry.MultiPolygon([polys[0]])
        elif geom.geom_type == "Polygon":
            if geom.area >= MIN_AREA:
                return shapely.geometry.MultiPolygon([geom])
            return shapely.geometry.MultiPolygon()
        return geom
    
    added = remove_slivers(added)
    removed = remove_slivers(removed)
    
    print(f"\nAfter cleaning:")
    print(f"  Added: area={added.area:.10f}" if not added.is_empty else "  Added: empty")
    print(f"  Removed: area={removed.area:.10f}" if not removed.is_empty else "  Removed: empty")
    print(f"  Unchanged: area={unchanged.area:.10f}" if not unchanged.is_empty else "  Unchanged: empty")
    
    # 4. Attribute fragments to districts and merge per district
    def get_district_fragments(geom, ctype, source_gdf):
        """Attribute geometry fragments to districts, merging by district name."""
        if geom.is_empty:
            return {}
        
        if geom.geom_type == "MultiPolygon":
            polys = list(geom.geoms)
        elif geom.geom_type == "Polygon":
            polys = [geom]
        else:
            return {}
        
        # Group fragments by district
        district_geoms = {}
        for poly in polys:
            if poly.area < MIN_AREA:
                continue
            if poly.geom_type == "Polygon":
                mp = shapely.geometry.MultiPolygon([poly])
            else:
                mp = poly
            
            # Find which districts this intersects
            districts = set()
            for _, row in source_gdf.iterrows():
                if mp.intersects(row["geometry"]):
                    districts.add(row["district_name"])
            
            dname = ", ".join(sorted(districts)) if districts else "unknown"
            
            if dname not in district_geoms:
                district_geoms[dname] = []
            district_geoms[dname].append(poly)
        
        return district_geoms
    
    added_by_district = get_district_fragments(added, "added_since_2014", current_gdf)
    removed_by_district = get_district_fragments(removed, "removed_since_2014", historical_gdf)
    
    fragments = []
    
    # Merge fragments per district
    for dname, polys in added_by_district.items():
        if len(polys) == 1:
            merged = shapely.geometry.MultiPolygon([polys[0]])
        else:
            merged = shapely.union_all(polys)
            if merged.geom_type == "Polygon":
                merged = shapely.geometry.MultiPolygon([merged])
        fragments.append({
            "geometry": merged,
            "change_type": "added_since_2014",
            "district_name": dname,
        })
    
    for dname, polys in removed_by_district.items():
        if len(polys) == 1:
            merged = shapely.geometry.MultiPolygon([polys[0]])
        else:
            merged = shapely.union_all(polys)
            if merged.geom_type == "Polygon":
                merged = shapely.geometry.MultiPolygon([merged])
        fragments.append({
            "geometry": merged,
            "change_type": "removed_since_2014",
            "district_name": dname,
        })
    
    # Unchanged: per district
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        dist_geom = row["geometry"]
        unchanged_part = shapely.intersection(dist_geom, unchanged)
        if not unchanged_part.is_empty:
            unchanged_part = remove_slivers(unchanged_part)
            if not unchanged_part.is_empty:
                if unchanged_part.geom_type == "MultiPolygon":
                    merged = shapely.union_all(list(unchanged_part.geoms))
                elif unchanged_part.geom_type == "Polygon":
                    merged = unchanged_part
                else:
                    continue
                if merged.geom_type == "Polygon":
                    merged = shapely.geometry.MultiPolygon([merged])
                fragments.append({
                    "geometry": merged,
                    "change_type": "unchanged",
                    "district_name": dname,
                })
    
    result_gdf = gpd.GeoDataFrame(fragments, crs="EPSG:4326")
    
    if result_gdf.empty:
        print("No fragments produced. Using current boundaries as unchanged.")
        result_gdf = current_gdf.copy()
        result_gdf["change_type"] = "unchanged"
        result_gdf = result_gdf[["geometry", "change_type", "district_name"]]
    
    print(f"\nFragments after merging: {len(result_gdf)}")
    if not result_gdf.empty:
        print(result_gdf["change_type"].value_counts())
    
    # -----------------------------------------------------------------------
    # 4. touches_changed
    # -----------------------------------------------------------------------
    changed_geoms = []
    for _, row in result_gdf.iterrows():
        if row["change_type"] in ("added_since_2014", "removed_since_2014"):
            changed_geoms.append(row["geometry"])
    
    if changed_geoms:
        changed_union = shapely.union_all(changed_geoms)
    else:
        changed_union = shapely.geometry.MultiPolygon()
    
    district_touches = {}
    for _, row in current_gdf.iterrows():
        dname = row["district_name"]
        geom = row["geometry"]
        if not changed_union.is_empty:
            touches = geom.touches(changed_union)
            if not touches:
                overlap = shapely.intersection(geom, changed_union)
                touches = not overlap.is_empty
            district_touches[dname] = touches
        else:
            district_touches[dname] = False
    
    for _, row in historical_gdf.iterrows():
        dname = row["district_name"]
        if dname not in district_touches:
            geom = row["geometry"]
            if not changed_union.is_empty:
                touches = geom.touches(changed_union)
                if not touches:
                    overlap = shapely.intersection(geom, changed_union)
                    touches = not overlap.is_empty
                district_touches[dname] = touches
            else:
                district_touches[dname] = False
    
    def get_touches(dn_str):
        if pd.isna(dn_str) or dn_str == "unknown":
            return False
        for d in dn_str.split(","):
            d = d.strip()
            if district_touches.get(d, False):
                return True
        return False
    
    result_gdf["touches_changed"] = result_gdf["district_name"].apply(get_touches)
    
    # Simplify and ensure MultiPolygon
    result_gdf["geometry"] = result_gdf["geometry"].simplify(0.00001, preserve_topology=True)
    result_gdf["geometry"] = result_gdf["geometry"].buffer(0)
    
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == "Polygon":
            return shapely.geometry.MultiPolygon([geom])
        return geom
    
    result_gdf["geometry"] = result_gdf["geometry"].apply(ensure_multipolygon)
    
    # -----------------------------------------------------------------------
    # 5. Output
    # -----------------------------------------------------------------------
    output_path = "/work/vienna_boundary_changes.geojson"
    result_gdf.to_file(output_path, driver="GeoJSON")
    
    print(f"\n{'=' * 60}")
    print(f"Output written to {output_path}")
    print(f"Features: {len(result_gdf)}")
    print(f"Columns: {result_gdf.columns.tolist()}")
    
    print(f"\nAll features:")
    print(result_gdf[["change_type", "district_name", "touches_changed"]])
    
    print(f"\nDistrict touch summary:")
    for dname, touches in sorted(district_touches.items()):
        print(f"  {dname}: touches_changed={touches}")
    
    # Verify
    verify = gpd.read_file(output_path)
    print(f"\nVerification: read back {len(verify)} features")
    print(verify["change_type"].value_counts())
    print(f"Geometry types: {verify.geometry.geom_type.unique()}")
    
    for ct in ["added_since_2014", "removed_since_2014", "unchanged"]:
        sub = verify[verify["change_type"] == ct]
        if not sub.empty:
            areas = sub.geometry.area
            print(f"  {ct}: count={len(sub)}, total_area={areas.sum():.10f}, max={areas.max():.8f}")


if __name__ == "__main__":
    main()