dc-l3-vienna-overpass-historical

Rerun

#!/usr/bin/env python3
"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass,
classify boundary changes, and output vienna_boundary_changes.geojson.
"""

import geopandas as gpd
import pandas as pd
import osmnx as ox
import numpy as np
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union, polygonize
from shapely.validation import make_valid
import re
import time
import json

# ── 1. Fetch current Vienna districts via OSMnx ──────────────────────────────

print("Step 1: Fetching current Vienna districts...")
# Vienna has admin_level=9 districts under the Vienna relation (admin_level=8)
# Use a bounding box around Vienna to be safe
try:
    current_districts = ox.features_from_place(
        "Wien, Österreich",
        tags={"admin_level": "9", "boundary": "administrative"},
    )
    print(f"  → Got {len(current_districts)} current district features")
except Exception as e:
    print(f"  OSMnx failed: {e}, trying Overpass directly...")
    # Fallback: use overpass query directly
    import overpy
    api = overpy.Overpass()
    q = """
    [out:json][timeout:60];
    area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
    rel(area.a)["admin_level"="9"]["boundary"="administrative"];
    out geom;
    """
    result = api.query(q)
    # Build GeoDataFrame manually...
    # For now, raise
    raise

# Get only unique districts (sometimes OSMnx returns overlapping results)
# Keep only the first occurrence per relation id
if "element_type" in current_districts.columns:
    current_districts = current_districts[
        ~current_districts.index.duplicated(keep="first")
    ]

print(f"  → {len(current_districts)} districts after dedup")

# ── 2. Fetch 2014-01-01 Vienna districts ─────────────────────────────────────

print("Step 2: Fetching 2014-01-01 Vienna districts...")

# We need to query Overpass with a date filter. OSMnx doesn't support this
# directly, so we'll use overpy with [date:"2014-01-01T00:00:00Z"]

try:
    import overpy
except ImportError:
    import subprocess
    subprocess.run(["pip", "install", "--user", "overpy"], check=True)
    import overpy

# First, get the relation IDs of Vienna districts
api = overpy.Overpass()
# Cache the current relations to avoid double fetching
current_rels = list(current_districts.index) if hasattr(current_districts.index, '__iter__') else []

# Query for current relations first to get IDs
q_current = """
[out:json][timeout:120];
area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
rel(area.a)["admin_level"="9"]["boundary"="administrative"];
out ids;
"""
result_current = api.query(q_current)
current_rel_ids = [r.id for r in result_current.relations]
print(f"  Found {len(current_rel_ids)} current district relation IDs")

# Now query for the 2014-01-01 versions of these relations
old_features = []
for rid in current_rel_ids:
    q_old = f"""
    [out:json][timeout:60][date:"2014-01-01T00:00:00Z"];
    rel({rid});
    out geom;
    """
    try:
        result = api.query(q_old)
        if result.relations:
            r = result.relations[0]
            tags = dict(r.tags)
            members = r.members
            
            # Build geometry from ways
            # We need to collect all ways and their nodes
            parts = []
            for m in members:
                if m.role == "outer" and hasattr(m, 'geometry') and m.geometry:
                    coords = [(float(g.lon), float(g.lat)) for g in m.geometry]
                    if len(coords) >= 3:
                        parts.append(Polygon(coords))
            
            if parts:
                # Merge all outer parts into a single geometry
                merged = unary_union(parts)
                if isinstance(merged, Polygon):
                    merged = MultiPolygon([merged])
                elif isinstance(merged, MultiPolygon) and len(merged.geoms) == 0:
                    continue
                elif isinstance(merged, MultiPolygon):
                    pass
                else:
                    continue
                
                old_features.append({
                    "relation_id": rid,
                    "geometry": merged,
                    **tags
                })
                print(f"    rel {rid}: {tags.get('name', 'unnamed')}")
        
        time.sleep(0.3)  # Be nice to Overpass
        
    except Exception as e:
        print(f"    rel {rid}: SKIPPED ({e})")
        continue

print(f"  → Got {len(old_features)} historical district features")

if len(old_features) == 0:
    print("WARNING: No 2014 data retrieved. Using simpler approach...")
    # Fallback: try one big query for historical data
    q_bulk = """
    [out:json][timeout:120][date:"2014-01-01T00:00:00Z"];
    area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
    rel(area.a)["admin_level"="9"]["boundary"="administrative"];
    out geom;
    """
    try:
        result = api.query(q_bulk)
        print(f"  Bulk query returned {len(result.relations)} relations")
        for r in result.relations:
            tags = dict(r.tags)
            parts = []
            for m in r.members:
                if hasattr(m, 'role') and m.role == "outer" and hasattr(m, 'geometry') and m.geometry:
                    coords = [(float(g.lon), float(g.lat)) for g in m.geometry]
                    if len(coords) >= 3:
                        parts.append(Polygon(coords))
            if parts:
                merged = unary_union(parts)
                if isinstance(merged, Polygon):
                    merged = MultiPolygon([merged])
                elif isinstance(merged, MultiPolygon):
                    pass
                else:
                    continue
                old_features.append({
                    "relation_id": r.id,
                    "geometry": merged,
                    **tags
                })
    except Exception as e2:
        print(f"  Bulk query also failed: {e2}")

print(f"  → Total old districts: {len(old_features)}")

# ── 3. Build GeoDataFrames ───────────────────────────────────────────────────

print("Step 3: Building GeoDataFrames...")

# Current districts
current_gdf = current_districts.copy()
current_gdf = current_gdf.to_crs("EPSG:4326")
# Ensure we have a clean dataframe
current_gdf = current_gdf.reset_index()
# Keep only useful columns
cols_to_keep = [c for c in ['name', 'geometry', 'admin_level'] if c in current_gdf.columns]
if 'name' not in current_gdf.columns:
    # Try alternative name columns
    for c in ['name:de', 'name:en', 'official_name']:
        if c in current_gdf.columns:
            current_gdf['name'] = current_gdf[c]
            break
    if 'name' not in current_gdf.columns:
        current_gdf['name'] = 'unknown'
current_gdf = current_gdf[cols_to_keep + ['element_type', 'osmid']] if 'element_type' in current_gdf.columns and 'osmid' in current_gdf.columns else current_gdf[cols_to_keep]
current_gdf = current_gdf.dropna(subset=['geometry'])
# Fix invalid geometries
current_gdf['geometry'] = current_gdf['geometry'].apply(lambda g: make_valid(g) if g and not g.is_valid else g)
print(f"  Current: {len(current_gdf)} features")

# Historical districts
old_gdf = gpd.GeoDataFrame(old_features, crs="EPSG:4326")
if 'name' not in old_gdf.columns:
    if 'name:de' in old_gdf.columns:
        old_gdf['name'] = old_gdf['name:de']
    elif 'official_name' in old_gdf.columns:
        old_gdf['name'] = old_gdf['official_name']
    else:
        old_gdf['name'] = 'unknown'
old_gdf = old_gdf.dropna(subset=['geometry'])
old_gdf['geometry'] = old_gdf['geometry'].apply(lambda g: make_valid(g) if g and not g.is_valid else g)
print(f"  Historical (2014): {len(old_gdf)} features")

# ── 4. Normalise district names ──────────────────────────────────────────────

print("Step 4: Normalising district names...")

def normalise_name(name):
    """Normalise district names to handle casing and spelling differences."""
    if not isinstance(name, str):
        return ""
    # Lowercase
    s = name.lower().strip()
    # Remove leading zeros from numbers (e.g. "01." -> "1.")
    s = re.sub(r'\b0+(\d)', r'\1', s)
    # Normalise "bezirk" variations
    s = re.sub(r'bezirk$', 'bezirk', s)
    # Remove extra whitespace
    s = re.sub(r'\s+', ' ', s).strip()
    # Remove trailing dot on ordinal numbers (1. -> 1)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    # Handle "Wien" prefix/suffix
    s = re.sub(r'^wien[-\s]', '', s)
    s = re.sub(r'[-\s]wien$', '', s)
    # Handle common abbreviations
    s = s.replace('st.', 'sankt').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor').replace('dr ', 'doktor ')
    s = s.replace('u.', 'und')
    # Collapse whitespace again
    s = re.sub(r'\s+', ' ', s).strip()
    return s

# Apply normalisation
current_gdf['name_norm'] = current_gdf['name'].apply(normalise_name)
old_gdf['name_norm'] = old_gdf['name'].apply(normalise_name)

print("  Current names:", sorted(current_gdf['name_norm'].unique()))
print("  Historical names:", sorted(old_gdf['name_norm'].unique()))

# ── 5. Match districts and compute changes ───────────────────────────────────

print("Step 5: Matching districts and computing changes...")

# Build a name mapping
# First, try exact match on normalised name
old_by_name = {}
for _, row in old_gdf.iterrows():
    old_by_name[row['name_norm']] = row['geometry']

# Try to match each current district to a historical one
matches = []
unmatched_current = []
unmatched_old_names = set(old_by_name.keys())

for _, row in current_gdf.iterrows():
    nm = row['name_norm']
    if nm in old_by_name:
        matches.append((row, old_by_name[nm], nm))
        unmatched_old_names.discard(nm)
    else:
        unmatched_current.append(row)

# For unmatched, try fuzzy matching
print(f"  Matched: {len(matches)}, Unmatched current: {len(unmatched_current)}, Unmatched old: {len(unmatched_old_names)}")

# Simple fallback: try matching by substrings
if unmatched_current or unmatched_old_names:
    old_names_avail = list(old_by_name.keys())
    for cur_row in list(unmatched_current):
        nm = cur_row['name_norm']
        # Try to find a name that starts with the same district number
        cur_num = re.match(r'^(\d+)', nm)
        if cur_num:
            cur_num = cur_num.group(1)
            for on in old_names_avail:
                old_num = re.match(r'^(\d+)', on)
                if old_num and old_num.group(1) == cur_num:
                    matches.append((cur_row, old_by_name[on], nm))
                    unmatched_current = [r for r in unmatched_current if r['name_norm'] != nm]
                    if on in unmatched_old_names:
                        unmatched_old_names.discard(on)
                    break

print(f"  After fuzzy: Matched: {len(matches)}, Unmatched current: {len(unmatched_current)}, Unmatched old: {len(unmatched_old_names)}")

# ── 6. Compute difference fragments ──────────────────────────────────────────

print("Step 6: Computing geometry difference fragments...")

change_fragments = []

# For matched districts, compute symmetric difference
for cur_row, old_geom, norm_name in matches:
    cur_geom = cur_row['geometry']
    
    # Ensure both are valid
    if not cur_geom.is_valid:
        cur_geom = make_valid(cur_geom)
    if not old_geom.is_valid:
        old_geom = make_valid(old_geom)
    
    # Ensure both are MultiPolygon
    if isinstance(cur_geom, Polygon):
        cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon):
        old_geom = MultiPolygon([old_geom])
    
    # Added: in current but not in old
    try:
        added = cur_geom.difference(old_geom)
    except Exception as e:
        print(f"  Error computing difference for {norm_name}: {e}")
        added = MultiPolygon()
    
    # Removed: in old but not in current
    try:
        removed = old_geom.difference(cur_geom)
    except Exception as e:
        print(f"  Error computing difference for {norm_name}: {e}")
        removed = MultiPolygon()
    
    # Unchanged: intersection
    try:
        unchanged = cur_geom.intersection(old_geom)
    except Exception as e:
        print(f"  Error computing intersection for {norm_name}: {e}")
        unchanged = MultiPolygon()
    
    # Normalise to MultiPolygon
    for geom, ctype in [(added, 'added_since_2014'), (removed, 'removed_since_2014'), (unchanged, 'unchanged')]:
        if geom.is_empty:
            continue
        if isinstance(geom, Polygon):
            geom = MultiPolygon([geom])
        if not isinstance(geom, MultiPolygon):
            continue
        if geom.is_empty:
            continue
        
        change_fragments.append({
            'geometry': geom,
            'change_type': ctype,
            'district_name': norm_name,
            'touches_changed': ctype != 'unchanged'  # Default - we'll refine below
        })

# For unmatched current districts (only exist now), everything is added
for cur_row in unmatched_current:
    geom = cur_row['geometry']
    if geom.is_empty:
        continue
    if isinstance(geom, Polygon):
        geom = MultiPolygon([geom])
    if not isinstance(geom, MultiPolygon):
        continue
    change_fragments.append({
        'geometry': geom,
        'change_type': 'added_since_2014',
        'district_name': cur_row['name_norm'],
        'touches_changed': True
    })

# For unmatched old districts (only existed in 2014), everything is removed
for on in unmatched_old_names:
    geom = old_by_name[on]
    if geom.is_empty:
        continue
    if isinstance(geom, Polygon):
        geom = MultiPolygon([geom])
    if not isinstance(geom, MultiPolygon):
        continue
    change_fragments.append({
        'geometry': geom,
        'change_type': 'removed_since_2014',
        'district_name': on,
        'touches_changed': True
    })

# ── 7. Compute touches_changed ───────────────────────────────────────────────

print("Step 7: Computing touches_changed flags...")

# For each district, determine if its geometry (current or historical) touches
# any added or removed area

# Collect all changed areas (added + removed) as a single geometry
changed_geoms = []
for frag in change_fragments:
    if frag['change_type'] in ('added_since_2014', 'removed_since_2014'):
        changed_geoms.append(frag['geometry'])

changed_union = unary_union(changed_geoms) if changed_geoms else MultiPolygon()
print(f"  Changed area count: {len(changed_geoms)}")

# Build a mapping of district_name -> all geometries (current + historical)
district_geoms = {}
for cur_row, old_geom, norm_name in matches:
    if norm_name not in district_geoms:
        district_geoms[norm_name] = []
    district_geoms[norm_name].append(cur_row['geometry'])
    district_geoms[norm_name].append(old_geom)

for cur_row in unmatched_current:
    nm = cur_row['name_norm']
    if nm not in district_geoms:
        district_geoms[nm] = []
    district_geoms[nm].append(cur_row['geometry'])

for on in unmatched_old_names:
    if on not in district_geoms:
        district_geoms[on] = []
    district_geoms[on].append(old_by_name[on])

# For each fragment, determine if its district touches any changed area
for frag in change_fragments:
    dn = frag['district_name']
    if dn in district_geoms:
        district_union = unary_union([g for g in district_geoms[dn] if not g.is_empty])
        if not district_union.is_empty and not changed_union.is_empty:
            try:
                touches = district_union.touches(changed_union)
                if touches:
                    frag['touches_changed'] = True
                else:
                    # Also check if the fragment itself is a changed area
                    frag['touches_changed'] = (frag['change_type'] in ('added_since_2014', 'removed_since_2014'))
            except Exception:
                frag['touches_changed'] = (frag['change_type'] in ('added_since_2014', 'removed_since_2014'))
        else:
            frag['touches_changed'] = False
    else:
        frag['touches_changed'] = (frag['change_type'] in ('added_since_2014', 'removed_since_2014'))

# ── 8. Explode MultiPolygons into individual features ─────────────────────────

print("Step 8: Exploding MultiPolygons into individual features...")

exploded_fragments = []
for frag in change_fragments:
    geom = frag['geometry']
    if geom.is_empty:
        continue
    if isinstance(geom, Polygon):
        geom = MultiPolygon([geom])
    if not isinstance(geom, MultiPolygon):
        continue
    
    # If MultiPolygon has multiple parts, keep them together as one feature
    # (the spec says MultiPolygon, so we keep the MultiPolygon intact)
    exploded_fragments.append({
        'geometry': geom,
        'change_type': frag['change_type'],
        'district_name': frag['district_name'],
        'touches_changed': frag['touches_changed']
    })

# ── 9. Build output GeoDataFrame and save ────────────────────────────────────

print("Step 9: Saving to GeoJSON...")

if not exploded_fragments:
    print("ERROR: No fragments to output!")
    # Create a dummy empty output
    empty_gdf = gpd.GeoDataFrame({
        'change_type': pd.Series(dtype='str'),
        'district_name': pd.Series(dtype='str'),
        'touches_changed': pd.Series(dtype='bool'),
        'geometry': pd.Series(dtype='object'),
    }, crs="EPSG:4326")
else:
    out_gdf = gpd.GeoDataFrame(exploded_fragments, crs="EPSG:4326")
    out_gdf = out_gdf.set_geometry('geometry')
    
    # Ensure MultiPolygon type
    out_gdf['geometry'] = out_gdf['geometry'].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )
    
    out_gdf.to_file(
        "/work/vienna_boundary_changes.geojson",
        driver="GeoJSON",
        engine="pyogrio"
    )
    print(f"  Saved {len(out_gdf)} features to vienna_boundary_changes.geojson")
    print(f"  Change types: {out_gdf['change_type'].value_counts().to_dict()}")
    print(f"  Districts: {sorted(out_gdf['district_name'].unique())}")

print("Done!")

#!/usr/bin/env python3
"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass,
classify boundary changes, and output vienna_boundary_changes.geojson.
"""

import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid
import re
import time
import sys
import json

# ── Retry helper ──────────────────────────────────────────────────────────────

def overpass_query(query, max_retries=3, retry_delay=3):
    """Run an Overpass query, trying multiple endpoints."""
    import overpy
    
    endpoints = [
        "https://overpass-api.de/api/interpreter",
        "https://lz4.overpass-api.de/api/interpreter",
        "https://overpass.kumi.systems/api/interpreter",
    ]
    
    last_err = None
    for endpoint in endpoints:
        for attempt in range(max_retries):
            try:
                api = overpy.Overpass(url=endpoint)
                result = api.query(query)
                return result
            except Exception as e:
                last_err = e
                msg = str(e).lower()
                # Don't retry if it's a query error (not connection issue)
                if "parse" in msg or "syntax" in msg or "timeout" in msg and "runtime error" in msg:
                    break
                print(f"  Retry {attempt+1}/{max_retries} on {endpoint}: {e}")
                time.sleep(retry_delay * (attempt + 1))
        else:
            continue
        break
    
    raise RuntimeError(f"Overpass query failed after all retries: {last_err}")

# ── 1. Fetch current Vienna districts ─────────────────────────────────────────

print("Step 1: Fetching current Vienna districts...")

# First get the relation IDs of Vienna districts from current data
q_current_ids = """
[out:json][timeout:120];
area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
rel(area.a)["admin_level"="9"]["boundary"="administrative"];
out ids;
"""
result_current = overpass_query(q_current_ids)
current_rel_ids = [r.id for r in result_current.relations]
print(f"  Found {len(current_rel_ids)} current district relation IDs")

# Now fetch full geometry for current districts
current_features = []
for rid in current_rel_ids:
    q_geom = f"""
    [out:json][timeout:60];
    rel({rid});
    (._;>;);
    out geom;
    """
    try:
        result = overpass_query(q_geom)
        if result.relations:
            r = result.relations[0]
            tags = dict(r.tags)
            members = r.members
            
            parts = []
            for m in members:
                if isinstance(m, overpy.RelationMember) and m.role == "outer":
                    # If the member is a way and has geometry
                    if hasattr(m, 'resolve') and m.resolve is not None:
                        way = m.resolve
                        if hasattr(way, 'nodes'):
                            coords = [(float(n.lon), float(n.lat)) for n in way.nodes]
                            if len(coords) >= 3:
                                parts.append(Polygon(coords))
            
            if not parts:
                # Alternative: try getting geometry from the relation's members directly
                for m in members:
                    if isinstance(m, overpy.RelationNode) or (hasattr(m, 'geometry') and m.geometry):
                        # For older overpy, members may have geometry directly
                        pass
                
                # Last resort: extract from the way geometries in the result
                # Collect all ways that are part of the relation with role outer
                for m in r.members:
                    if hasattr(m, 'role') and m.role == 'outer' and hasattr(m, 'geometry') and m.geometry:
                        coords = [(float(g.lon), float(g.lat)) for g in m.geometry]
                        if len(coords) >= 3:
                            parts.append(Polygon(coords))
            
            if parts:
                merged = unary_union(parts)
                if isinstance(merged, Polygon):
                    merged = MultiPolygon([merged])
                if not merged.is_empty and isinstance(merged, MultiPolygon):
                    current_features.append({
                        "relation_id": rid,
                        "geometry": merged,
                        **tags
                    })
                    name = tags.get('name', 'unnamed')
                    print(f"    Current rel {rid}: {name}")
        
        time.sleep(0.4)
        
    except Exception as e:
        print(f"    Current rel {rid}: SKIPPED ({e})")
        continue

print(f"  → Got {len(current_features)} current district features")

if len(current_features) < 10:
    print("  Not enough current features, trying bulk query...")
    q_bulk_current = """
    [out:json][timeout:180];
    area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
    rel(area.a)["admin_level"="9"]["boundary"="administrative"];
    out geom;
    """
    result_bulk = overpass_query(q_bulk_current)
    print(f"  Bulk returned {len(result_bulk.relations)} relations")
    
    # Clear and redo
    current_features = []
    for r in result_bulk.relations:
        tags = dict(r.tags)
        parts = []
        for m in r.members:
            if hasattr(m, 'role') and m.role == 'outer' and hasattr(m, 'geometry') and m.geometry:
                coords = [(float(g.lon), float(g.lat)) for g in m.geometry]
                if len(coords) >= 3:
                    parts.append(Polygon(coords))
        if parts:
            merged = unary_union(parts)
            if isinstance(merged, Polygon):
                merged = MultiPolygon([merged])
            if not merged.is_empty and isinstance(merged, MultiPolygon):
                current_features.append({
                    "relation_id": r.id,
                    "geometry": merged,
                    **tags
                })
    
    print(f"  → Got {len(current_features)} from bulk")

# ── 2. Fetch 2014-01-01 Vienna districts ─────────────────────────────────────

print("\nStep 2: Fetching 2014-01-01 Vienna districts...")

# Try bulk query first for historical data
old_features = []
q_old_bulk = """
[out:json][timeout:180][date:"2014-01-01T00:00:00Z"];
area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
rel(area.a)["admin_level"="9"]["boundary"="administrative"];
out geom;
"""
try:
    result_old = overpass_query(q_old_bulk)
    print(f"  Bulk historical returned {len(result_old.relations)} relations")
    
    for r in result_old.relations:
        tags = dict(r.tags)
        parts = []
        for m in r.members:
            if hasattr(m, 'role') and m.role == 'outer' and hasattr(m, 'geometry') and m.geometry:
                coords = [(float(g.lon), float(g.lat)) for g in m.geometry]
                if len(coords) >= 3:
                    parts.append(Polygon(coords))
        if parts:
            merged = unary_union(parts)
            if isinstance(merged, Polygon):
                merged = MultiPolygon([merged])
            if not merged.is_empty and isinstance(merged, MultiPolygon):
                old_features.append({
                    "relation_id": r.id,
                    "geometry": merged,
                    **tags
                })
                print(f"    Historical rel {r.id}: {tags.get('name', 'unnamed')}")
    
except Exception as e:
    print(f"  Bulk historical failed: {e}")
    # Fallback: try individual queries
    for rid in current_rel_ids:
        q_old = f"""
        [out:json][timeout:60][date:"2014-01-01T00:00:00Z"];
        rel({rid});
        out geom;
        """
        try:
            result = overpass_query(q_old)
            if result.relations:
                r = result.relations[0]
                tags = dict(r.tags)
                parts = []
                for m in r.members:
                    if hasattr(m, 'role') and m.role == 'outer' and hasattr(m, 'geometry') and m.geometry:
                        coords = [(float(g.lon), float(g.lat)) for g in m.geometry]
                        if len(coords) >= 3:
                            parts.append(Polygon(coords))
                if parts:
                    merged = unary_union(parts)
                    if isinstance(merged, Polygon):
                        merged = MultiPolygon([merged])
                    if not merged.is_empty and isinstance(merged, MultiPolygon):
                        old_features.append({
                            "relation_id": rid,
                            "geometry": merged,
                            **tags
                        })
                        print(f"    Historical rel {rid}: {tags.get('name', 'unnamed')}")
            
            time.sleep(0.5)
            
        except Exception as e2:
            print(f"    Historical rel {rid}: SKIPPED ({e2})")
            continue

print(f"  → Got {len(old_features)} historical district features")

# ── 3. Build GeoDataFrames ───────────────────────────────────────────────────

print("\nStep 3: Building GeoDataFrames...")

# Current districts
current_gdf = gpd.GeoDataFrame(current_features, crs="EPSG:4326")
if 'name' not in current_gdf.columns:
    if 'name:de' in current_gdf.columns:
        current_gdf['name'] = current_gdf['name:de']
    elif 'official_name' in current_gdf.columns:
        current_gdf['name'] = current_gdf['official_name']
    else:
        current_gdf['name'] = 'unknown'
current_gdf = current_gdf.dropna(subset=['geometry'])
current_gdf['geometry'] = current_gdf['geometry'].apply(
    lambda g: make_valid(g) if g and not g.is_valid else g
)
print(f"  Current: {len(current_gdf)} features")

# Historical districts
old_gdf = gpd.GeoDataFrame(old_features, crs="EPSG:4326")
if 'name' not in old_gdf.columns:
    if 'name:de' in old_gdf.columns:
        old_gdf['name'] = old_gdf['name:de']
    elif 'official_name' in old_gdf.columns:
        old_gdf['name'] = old_gdf['official_name']
    else:
        old_gdf['name'] = 'unknown'
old_gdf = old_gdf.dropna(subset=['geometry'])
old_gdf['geometry'] = old_gdf['geometry'].apply(
    lambda g: make_valid(g) if g and not g.is_valid else g
)
print(f"  Historical (2014): {len(old_gdf)} features")

if len(current_gdf) < 2 or len(old_gdf) < 2:
    print(f"WARNING: Not enough data. Current={len(current_gdf)}, Old={len(old_gdf)}")
    print("Trying alternative query approach...")
    sys.exit(1)

# ── 4. Normalise district names ──────────────────────────────────────────────

print("\nStep 4: Normalising district names...")

def normalise_name(name):
    """Normalise district names to handle casing and spelling differences."""
    if not isinstance(name, str):
        return ""
    s = name.lower().strip()
    # Remove leading zeros from numbers (e.g. "01." -> "1.")
    s = re.sub(r'\b0+(\d)', r'\1', s)
    # Normalise "bezirk" variations
    s = re.sub(r'bezirk$', 'bezirk', s)
    # Remove extra whitespace
    s = re.sub(r'\s+', ' ', s).strip()
    # Remove trailing dot on ordinal numbers (1. -> 1)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    # Handle "Wien" prefix/suffix
    s = re.sub(r'^wien[-\s]', '', s)
    s = re.sub(r'[-\s]wien$', '', s)
    # Handle common abbreviations
    s = s.replace('st.', 'sankt').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor').replace('dr ', 'doktor ')
    s = s.replace('u.', 'und')
    # Collapse whitespace again
    s = re.sub(r'\s+', ' ', s).strip()
    return s

current_gdf['name_norm'] = current_gdf['name'].apply(normalise_name)
old_gdf['name_norm'] = old_gdf['name'].apply(normalise_name)

print("  Current names:", sorted(current_gdf['name_norm'].tolist()))
print("  Historical names:", sorted(old_gdf['name_norm'].tolist()))

# ── 5. Match districts ───────────────────────────────────────────────────────

print("\nStep 5: Matching districts...")

old_by_name = {}
for _, row in old_gdf.iterrows():
    old_by_name[row['name_norm']] = row['geometry']

matches = []
unmatched_current = []
unmatched_old_names = set(old_by_name.keys())

for _, row in current_gdf.iterrows():
    nm = row['name_norm']
    if nm in old_by_name:
        matches.append((row['geometry'], old_by_name[nm], nm, row.get('name', '')))
        unmatched_old_names.discard(nm)
    else:
        unmatched_current.append(row)

# Try fuzzy matching by district number
if unmatched_current or unmatched_old_names:
    old_names_avail = list(unmatched_old_names)
    for cur_row in list(unmatched_current):
        nm = cur_row['name_norm']
        cur_num = re.match(r'^(\d+)', nm)
        if cur_num:
            cur_num = cur_num.group(1)
            for on in list(unmatched_old_names):
                old_num = re.match(r'^(\d+)', on)
                if old_num and old_num.group(1) == cur_num:
                    matches.append((cur_row['geometry'], old_by_name[on], nm, cur_row.get('name', '')))
                    unmatched_current = [r for r in unmatched_current if r['name_norm'] != nm]
                    unmatched_old_names.discard(on)
                    break

print(f"  Matched: {len(matches)}")
print(f"  Unmatched current: {len(unmatched_current)}")
print(f"  Unmatched old: {len(unmatched_old_names)}")
if unmatched_current:
    print(f"    Current names: {[r['name_norm'] for r in unmatched_current]}")
if unmatched_old_names:
    print(f"    Old names: {list(unmatched_old_names)}")

# ── 6. Compute difference fragments ──────────────────────────────────────────

print("\nStep 6: Computing geometry differences...")

change_fragments = []

for cur_geom, old_geom, norm_name, orig_name in matches:
    if not cur_geom.is_valid:
        cur_geom = make_valid(cur_geom)
    if not old_geom.is_valid:
        old_geom = make_valid(old_geom)
    
    if isinstance(cur_geom, Polygon):
        cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon):
        old_geom = MultiPolygon([old_geom])
    
    try:
        added = cur_geom.difference(old_geom)
    except Exception as e:
        print(f"  Error difference({norm_name}): {e}")
        added = MultiPolygon()
    
    try:
        removed = old_geom.difference(cur_geom)
    except Exception as e:
        print(f"  Error difference({norm_name}): {e}")
        removed = MultiPolygon()
    
    try:
        unchanged = cur_geom.intersection(old_geom)
    except Exception as e:
        print(f"  Error intersection({norm_name}): {e}")
        unchanged = MultiPolygon()
    
    for geom, ctype in [
        (added, 'added_since_2014'),
        (removed, 'removed_since_2014'),
        (unchanged, 'unchanged')
    ]:
        if geom.is_empty:
            continue
        if isinstance(geom, Polygon):
            geom = MultiPolygon([geom])
        if not isinstance(geom, MultiPolygon) or geom.is_empty:
            continue
        
        change_fragments.append({
            'geometry': geom,
            'change_type': ctype,
            'district_name': norm_name,
            'touches_changed': False  # computed below
        })

for cur_row in unmatched_current:
    geom = cur_row['geometry']
    if geom.is_empty:
        continue
    if isinstance(geom, Polygon):
        geom = MultiPolygon([geom])
    if not isinstance(geom, MultiPolygon):
        continue
    change_fragments.append({
        'geometry': geom,
        'change_type': 'added_since_2014',
        'district_name': cur_row['name_norm'],
        'touches_changed': True
    })

for on in unmatched_old_names:
    geom = old_by_name[on]
    if geom.is_empty:
        continue
    if isinstance(geom, Polygon):
        geom = MultiPolygon([geom])
    if not isinstance(geom, MultiPolygon):
        continue
    change_fragments.append({
        'geometry': geom,
        'change_type': 'removed_since_2014',
        'district_name': on,
        'touches_changed': True
    })

print(f"  Total fragments: {len(change_fragments)}")

# ── 7. Compute touches_changed ───────────────────────────────────────────────

print("\nStep 7: Computing touches_changed...")

# Collect all changed geometries
changed_geoms = []
for frag in change_fragments:
    if frag['change_type'] in ('added_since_2014', 'removed_since_2014'):
        changed_geoms.append(frag['geometry'])

changed_union = unary_union(changed_geoms) if changed_geoms else MultiPolygon()
print(f"  Changed area count: {len(changed_geoms)}")

# Build district geometries mapping
district_geoms = {}
for cur_geom, old_geom, norm_name, orig_name in matches:
    if norm_name not in district_geoms:
        district_geoms[norm_name] = []
    district_geoms[norm_name].append(cur_geom)
    district_geoms[norm_name].append(old_geom)

for cur_row in unmatched_current:
    nm = cur_row['name_norm']
    if nm not in district_geoms:
        district_geoms[nm] = []
    district_geoms[nm].append(cur_row['geometry'])

for on in unmatched_old_names:
    if on not in district_geoms:
        district_geoms[on] = []
    district_geoms[on].append(old_by_name[on])

for frag in change_fragments:
    dn = frag['district_name']
    if dn in district_geoms and not changed_union.is_empty:
        district_union = unary_union([g for g in district_geoms[dn] if not g.is_empty])
        if not district_union.is_empty:
            try:
                touches = district_union.touches(changed_union)
                frag['touches_changed'] = touches or (frag['change_type'] in ('added_since_2014', 'removed_since_2014'))
            except Exception:
                frag['touches_changed'] = (frag['change_type'] in ('added_since_2014', 'removed_since_2014'))
        else:
            frag['touches_changed'] = False
    else:
        frag['touches_changed'] = (frag['change_type'] in ('added_since_2014', 'removed_since_2014'))

# ── 8. Final assembly and save ───────────────────────────────────────────────

print("\nStep 8: Saving to GeoJSON...")

if not change_fragments:
    print("ERROR: No fragments generated!")
    empty_gdf = gpd.GeoDataFrame({
        'change_type': pd.Series(dtype='str'),
        'district_name': pd.Series(dtype='str'),
        'touches_changed': pd.Series(dtype='bool'),
        'geometry': pd.Series(dtype='object'),
    }, crs="EPSG:4326")
    empty_gdf.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON")
    sys.exit(1)

out_gdf = gpd.GeoDataFrame(change_fragments, crs="EPSG:4326")
out_gdf = out_gdf.set_geometry('geometry')

# Ensure MultiPolygon
out_gdf['geometry'] = out_gdf['geometry'].apply(
    lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
)

out_gdf.to_file(
    "/work/vienna_boundary_changes.geojson",
    driver="GeoJSON",
    engine="pyogrio"
)

print(f"  Saved {len(out_gdf)} features to vienna_boundary_changes.geojson")
print(f"  Change type counts: {out_gdf['change_type'].value_counts().to_dict()}")
print(f"  Districts: {sorted(out_gdf['district_name'].unique())}")

# Quick sanity check
print(f"  touches_changed true: {out_gdf['touches_changed'].sum()}")
print(f"  touches_changed false: {(~out_gdf['touches_changed']).sum()}")

print("\nDone!")

#!/usr/bin/env python3
"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass,
classify boundary changes, and output vienna_boundary_changes.geojson.
"""
import os, re, sys, json, time, pickle, traceback

import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import MultiPolygon, Polygon, shape
from shapely.ops import unary_union
from shapely.validation import make_valid

# ── Retry helper ──────────────────────────────────────────────────────────────

def robust_overpass(query, max_retries=3, retry_delay=2):
    """Run Overpass query trying multiple endpoints."""
    import overpy
    endpoints = [
        "https://overpass-api.de/api/interpreter",
        "https://lz4.overpass-api.de/api/interpreter",
        "https://overpass.kumi.systems/api/interpreter",
    ]
    for ep in endpoints:
        for attempt in range(max_retries):
            try:
                api = overpy.Overpass(url=ep, read_chunk_size=64*1024)
                return api.query(query)
            except Exception as e:
                msg = str(e).lower()
                print(f"  [{ep.split('/')[2]}] attempt {attempt+1}: {e!r}")
                if "parse" in msg or "syntax" in msg:
                    raise
                time.sleep(retry_delay * (attempt + 1))
    raise RuntimeError("All Overpass endpoints failed")

# ── 1. Fetch current districts (using relation list + full geometry) ─────────

print("=" * 60)
print("STEP 1: Fetching current Vienna districts")
print("=" * 60)

CACHE_FILE = "/work/_cache_data.pkl"

if os.path.exists(CACHE_FILE):
    print("Loading from cache...")
    with open(CACHE_FILE, "rb") as f:
        cache = pickle.load(f)
    current_features = cache["current"]
    old_features = cache["old"]
else:
    # Get current relation IDs
    q_ids = """
    [out:json][timeout:90];
    area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
    rel(area.a)["admin_level"="9"]["boundary"="administrative"];
    out ids;
    """
    res = robust_overpass(q_ids)
    rel_ids = [r.id for r in res.relations]
    print(f"  Relation IDs: {rel_ids}")

    # Bulk fetch current - full geometry
    q_cur = """
    [out:json][timeout:180];
    area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
    rel(area.a)["admin_level"="9"]["boundary"="administrative"];
    out geom;
    """
    res_cur = robust_overpass(q_cur)
    print(f"  Current bulk: {len(res_cur.relations)} relations")

    current_features = []
    for r in res_cur.relations:
        tags = dict(r.tags)
        parts = []
        for m in r.members:
            if getattr(m, 'role', None) == 'outer' and hasattr(m, 'geometry') and m.geometry:
                coords = [(float(g.lon), float(g.lat)) for g in m.geometry]
                if len(coords) >= 3:
                    parts.append(Polygon(coords))
        if parts:
            merged = unary_union(parts)
            if isinstance(merged, Polygon):
                merged = MultiPolygon([merged])
            if not merged.is_empty and isinstance(merged, MultiPolygon):
                current_features.append({
                    "relation_id": r.id,
                    "geometry": merged,
                    **tags
                })
                print(f"    ✓ current rel {r.id}: {tags.get('name','?')}")

    print(f"  → {len(current_features)} current features")

    # ── 2. Fetch 2014-01-01 districts ─────────────────────────────────────────
    print()
    print("=" * 60)
    print("STEP 2: Fetching 2014-01-01 Vienna districts")
    print("=" * 60)

    q_old = """
    [out:json][timeout:180][date:"2014-01-01T00:00:00Z"];
    area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
    rel(area.a)["admin_level"="9"]["boundary"="administrative"];
    out geom;
    """
    try:
        res_old = robust_overpass(q_old)
    except Exception as e:
        print(f"  Bulk historical failed: {e}")
        print("  Trying per-relation...")
        res_old = None
        old_relations = []
        for rid in rel_ids:
            q_one = f"""
            [out:json][timeout:60][date:"2014-01-01T00:00:00Z"];
            rel({rid});
            out geom;
            """
            try:
                r2 = robust_overpass(q_one)
                if r2.relations:
                    old_relations.append(r2.relations[0])
                time.sleep(0.3)
            except Exception as e2:
                print(f"    rel {rid}: {e2}")
        # build a mock result
        class MockResult:
            pass
        res_old = MockResult()
        res_old.relations = old_relations
        print(f"  → {len(old_relations)} via per-rel")

    old_features = []
    for r in (res_old.relations if res_old else []):
        tags = dict(r.tags)
        parts = []
        for m in r.members:
            if getattr(m, 'role', None) == 'outer' and hasattr(m, 'geometry') and m.geometry:
                coords = [(float(g.lon), float(g.lat)) for g in m.geometry]
                if len(coords) >= 3:
                    parts.append(Polygon(coords))
        if parts:
            merged = unary_union(parts)
            if isinstance(merged, Polygon):
                merged = MultiPolygon([merged])
            if not merged.is_empty and isinstance(merged, MultiPolygon):
                old_features.append({
                    "relation_id": r.id,
                    "geometry": merged,
                    **tags
                })
                print(f"    ✓ old rel {r.id}: {tags.get('name','?')}")

    print(f"  → {len(old_features)} old features")

    # Save cache
    with open(CACHE_FILE, "wb") as f:
        pickle.dump({"current": current_features, "old": old_features}, f)
    print("  Cache saved.")

# ── 3. Build GeoDataFrames ───────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 3: Building GeoDataFrames")
print("=" * 60)

def name_from_tags(tags):
    for k in ('name', 'name:de', 'official_name'):
        v = tags.get(k)
        if v:
            return v
    return 'unknown'

for feat in current_features:
    feat['name'] = name_from_tags(feat)
for feat in old_features:
    feat['name'] = name_from_tags(feat)

cur_gdf = gpd.GeoDataFrame(current_features, crs="EPSG:4326")
old_gdf = gpd.GeoDataFrame(old_features, crs="EPSG:4326")

cur_gdf = cur_gdf.dropna(subset=['geometry'])
old_gdf = old_gdf.dropna(subset=['geometry'])

cur_gdf['geometry'] = cur_gdf['geometry'].apply(
    lambda g: make_valid(g) if not g.is_valid else g
)
old_gdf['geometry'] = old_gdf['geometry'].apply(
    lambda g: make_valid(g) if not g.is_valid else g
)

print(f"  Current: {len(cur_gdf)} districts")
print(f"  Historical 2014: {len(old_gdf)} districts")

if len(cur_gdf) < 5 or len(old_gdf) < 5:
    print("ERROR: Too few districts fetched.")
    sys.exit(1)

# ── 4. Normalise names ──────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 4: Normalising district names")
print("=" * 60)

def norm(name):
    if not isinstance(name, str):
        return ""
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)       # 01. → 1.
    s = re.sub(r'bezirk$', 'bezirk', s)
    s = re.sub(r'\s+', ' ', s).strip()
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)    # 1.  → 1
    s = re.sub(r'^wien[-\s]', '', s)
    s = re.sub(r'[-\s]wien$', '', s)
    s = s.replace('st.', 'sankt').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor').replace('dr ', 'doktor ')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s

cur_gdf['name_norm'] = cur_gdf['name'].apply(norm)
old_gdf['name_norm'] = old_gdf['name'].apply(norm)

print("  Current names:", sorted(cur_gdf['name_norm'].tolist()))
print("  Historical names:", sorted(old_gdf['name_norm'].tolist()))

# ── 5. Match districts ───────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 5: Matching districts")
print("=" * 60)

old_lookup = {row['name_norm']: row['geometry'] for _, row in old_gdf.iterrows()}

matches = []
unmatched_cur = []
unmatched_old = set(old_lookup.keys())

for _, row in cur_gdf.iterrows():
    nm = row['name_norm']
    if nm in old_lookup:
        matches.append((row['geometry'], old_lookup[nm], nm))
        unmatched_old.discard(nm)
    else:
        unmatched_cur.append(row)

# Fuzzy by district number
if unmatched_cur or unmatched_old:
    for cur_row in list(unmatched_cur):
        nm = cur_row['name_norm']
        m = re.match(r'^(\d+)', nm)
        if m:
            n = m.group(1)
            for on in list(unmatched_old):
                m2 = re.match(r'^(\d+)', on)
                if m2 and m2.group(1) == n:
                    matches.append((cur_row['geometry'], old_lookup[on], nm))
                    unmatched_cur = [r for r in unmatched_cur if r['name_norm'] != nm]
                    unmatched_old.discard(on)
                    break

print(f"  Matched: {len(matches)}")
print(f"  Unmatched current: {len(unmatched_cur)} ({[r['name_norm'] for r in unmatched_cur]})")
print(f"  Unmatched old: {len(unmatched_old)} ({list(unmatched_old)})")

# ── 6. Compute differences ───────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 6: Computing geometry differences")
print("=" * 60)

fragments = []

for cur_geom, old_geom, nm in matches:
    if not cur_geom.is_valid:
        cur_geom = make_valid(cur_geom)
    if not old_geom.is_valid:
        old_geom = make_valid(old_geom)

    if isinstance(cur_geom, Polygon):
        cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon):
        old_geom = MultiPolygon([old_geom])

    try:
        added = cur_geom.difference(old_geom)
    except Exception:
        added = MultiPolygon()
    try:
        removed = old_geom.difference(cur_geom)
    except Exception:
        removed = MultiPolygon()
    try:
        unchanged = cur_geom.intersection(old_geom)
    except Exception:
        unchanged = MultiPolygon()

    for g, ct in [(added, 'added_since_2014'), (removed, 'removed_since_2014'), (unchanged, 'unchanged')]:
        if g.is_empty:
            continue
        if isinstance(g, Polygon):
            g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon) or g.is_empty:
            continue
        fragments.append({'geometry': g, 'change_type': ct, 'district_name': nm, 'touches_changed': False})

for row in unmatched_cur:
    g = row['geometry']
    if g.is_empty:
        continue
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon):
        continue
    fragments.append({'geometry': g, 'change_type': 'added_since_2014', 'district_name': row['name_norm'], 'touches_changed': True})

for on in unmatched_old:
    g = old_lookup[on]
    if g.is_empty:
        continue
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon):
        continue
    fragments.append({'geometry': g, 'change_type': 'removed_since_2014', 'district_name': on, 'touches_changed': True})

print(f"  Fragments: {len(fragments)}")

# ── 7. Compute touches_changed ───────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 7: Computing touches_changed")
print("=" * 60)

# Build union of all changed areas
changed_polys = [f['geometry'] for f in fragments if f['change_type'] in ('added_since_2014', 'removed_since_2014')]
changed_union = unary_union(changed_polys) if changed_polys else MultiPolygon()
print(f"  Changed area polygons: {len(changed_polys)}")

# Build per-district union
district_polys = {}
for cur_geom, old_geom, nm in matches:
    district_polys.setdefault(nm, []).append(cur_geom)
    district_polys.setdefault(nm, []).append(old_geom)
for row in unmatched_cur:
    district_polys.setdefault(row['name_norm'], []).append(row['geometry'])
for on in unmatched_old:
    district_polys.setdefault(on, []).append(old_lookup[on])

for f in fragments:
    dn = f['district_name']
    if dn in district_polys and not changed_union.is_empty:
        du = unary_union([g for g in district_polys[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                touches = du.touches(changed_union)
                f['touches_changed'] = touches or (f['change_type'] in ('added_since_2014', 'removed_since_2014'))
            except Exception:
                f['touches_changed'] = f['change_type'] in ('added_since_2014', 'removed_since_2014')
        else:
            f['touches_changed'] = False
    else:
        f['touches_changed'] = f['change_type'] in ('added_since_2014', 'removed_since_2014')

tc_true = sum(1 for f in fragments if f['touches_changed'])
print(f"  touches_changed = True: {tc_true}")

# ── 8. Output ────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 8: Saving to GeoJSON")
print("=" * 60)

if not fragments:
    print("ERROR: No fragments. Creating empty output.")
    out = gpd.GeoDataFrame({'change_type': pd.Series(dtype='str'),
                            'district_name': pd.Series(dtype='str'),
                            'touches_changed': pd.Series(dtype='bool'),
                            'geometry': pd.Series(dtype='object')},
                           crs="EPSG:4326")
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326")
    out = out.set_geometry('geometry')
    out['geometry'] = out['geometry'].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )
    # Force MultiPolygon type
    out['geometry'] = out['geometry'].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Wrote {len(out)} features")

# Summary
if len(out):
    print(f"  change_type distribution:")
    print(f"    {out['change_type'].value_counts().to_dict()}")
    print(f"  District names: {sorted(out['district_name'].unique())}")
    print(f"  touches_changed True count: {out['touches_changed'].sum()}")

print()
print("Done!")

#!/usr/bin/env python3
"""
Fetch Vienna district boundaries (current and 2014-01-01) from Overpass,
classify boundary changes, output vienna_boundary_changes.geojson.

Each feature: MultiPolygon, with change_type, district_name (normalised), touches_changed.
"""
import re, sys, time, json
import requests
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import MultiPolygon, Polygon, shape
from shapely.ops import unary_union
from shapely.validation import make_valid

UA = "GIS-Analyst-Agent/1.0"

def overpass_post(query, endpoint="https://overpass-api.de/api/interpreter",
                  timeout=120):
    """POST an Overpass QL query, return parsed JSON."""
    for attempt in range(3):
        try:
            r = requests.post(endpoint, data={"data": query},
                              headers={"User-Agent": UA},
                              timeout=timeout)
            r.raise_for_status()
            return r.json()
        except Exception as e:
            print(f"  Overpass attempt {attempt+1}: {e!r}", file=sys.stderr)
            time.sleep(2 * (attempt + 1))
    raise RuntimeError("Overpass failed after 3 retries")


def geom_from_members(members):
    """Build MultiPolygon from Overpass relation members (role=outer, .geom)."""
    parts = []
    for m in members:
        if m.get("role") == "outer" and "geometry" in m:
            coords = [(g["lon"], g["lat"]) for g in m["geometry"]]
            if len(coords) >= 3:
                parts.append(Polygon(coords))
    if not parts:
        return None
    merged = unary_union(parts)
    if isinstance(merged, Polygon):
        merged = MultiPolygon([merged])
    if merged.is_empty or not isinstance(merged, MultiPolygon):
        return None
    return merged


def norm_name(name):
    """Normalise district name for cross-snapshot matching."""
    if not isinstance(name, str):
        return ""
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)           # 01 → 1
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)       # 1.  → 1
    s = re.sub(r'bezirk$', 'bezirk', s)
    s = re.sub(r'^wien[-\s]', '', s)
    s = re.sub(r'[-\s]wien$', '', s)
    s = s.replace('st.', 'sankt').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s


# ── Helper: fetch relation features (bulk) ──────────────────────────────────

def fetch_districts(date_clause=""):
    """Fetch Vienna admin_level=9 relations, optionally with [date:...]."""
    q = f"""
    [out:json][timeout:180]{date_clause};
    area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
    rel(area.a)["admin_level"="9"]["boundary"="administrative"];
    out geom;
    """
    print(f"    Querying{' (date: 2014-01-01)' if date_clause else ''}...")
    j = overpass_post(q)
    feats = []
    for el in j.get("elements", []):
        if el["type"] != "relation":
            continue
        tags = el.get("tags", {})
        geom = geom_from_members(el.get("members", []))
        if geom is None:
            continue
        feats.append({"relation_id": el["id"], "geometry": geom, **tags})
    return feats


# ═══════════════════════════════════════════════════════════════════════════════
print("=" * 60)
print("STEP 1: Fetching current Vienna districts")
print("=" * 60)
current = fetch_districts("")

# Fallback if bulk returns nothing (try individual ways approach)
if len(current) < 10:
    print("  Bulk returned few results, trying alternative query...")
    q2 = """
    [out:json][timeout:180];
    area["name"="Wien"]["admin_level"="8"]["boundary"="administrative"]->.a;
    (rel(area.a)["admin_level"="9"]["boundary"="administrative"];
     way(r)["admin_level"="9"]["boundary"="administrative"];);
    out geom;
    """
    j2 = overpass_post(q2)
    for el in j2.get("elements", []):
        if el["type"] != "relation":
            continue
        tags = el.get("tags", {})
        geom = geom_from_members(el.get("members", []))
        if geom is None:
            continue
        current.append({"relation_id": el["id"], "geometry": geom, **tags})

print(f"  → {len(current)} current districts")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 2: Fetching 2014-01-01 Vienna districts")
print("=" * 60)
old = fetch_districts('[date:"2014-01-01T00:00:00Z"]')
print(f"  → {len(old)} historical districts")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 3: Building GeoDataFrames & matching")
print("=" * 60)

def name(tags):
    for k in ("name", "name:de", "official_name"):
        v = tags.get(k)
        if v:
            return v
    return "unknown"

for f in current:
    f["name"] = name(f)
for f in old:
    f["name"] = name(f)

cur_gdf = gpd.GeoDataFrame(current, crs="EPSG:4326").dropna(subset=["geometry"])
old_gdf = gpd.GeoDataFrame(old, crs="EPSG:4326").dropna(subset=["geometry"])

# Fix invalid geometries
cur_gdf["geometry"] = cur_gdf["geometry"].apply(
    lambda g: make_valid(g) if not g.is_valid else g
)
old_gdf["geometry"] = old_gdf["geometry"].apply(
    lambda g: make_valid(g) if not g.is_valid else g
)

cur_gdf["name_norm"] = cur_gdf["name"].apply(norm_name)
old_gdf["name_norm"] = old_gdf["name"].apply(norm_name)

print(f"  Current ({len(cur_gdf)}):  {sorted(cur_gdf['name_norm'].tolist())}")
print(f"  Old 2014 ({len(old_gdf)}): {sorted(old_gdf['name_norm'].tolist())}")

# Match
old_by_name = {row["name_norm"]: row["geometry"] for _, row in old_gdf.iterrows()}
matches = []
unmatched_cur = []
unmatched_old = set(old_by_name.keys())

for _, row in cur_gdf.iterrows():
    nm = row["name_norm"]
    if nm in old_by_name:
        matches.append((row["geometry"], old_by_name[nm], nm))
        unmatched_old.discard(nm)
    else:
        unmatched_cur.append(row)

# Fuzzy: match by district number prefix
for cur_row in list(unmatched_cur):
    nm = cur_row["name_norm"]
    m = re.match(r"^(\d+)", nm)
    if not m:
        continue
    n = m.group(1)
    for on in list(unmatched_old):
        m2 = re.match(r"^(\d+)", on)
        if m2 and m2.group(1) == n:
            matches.append((cur_row["geometry"], old_by_name[on], nm))
            unmatched_cur = [r for r in unmatched_cur if r["name_norm"] != nm]
            unmatched_old.discard(on)
            break

print(f"  Matched: {len(matches)}")
if unmatched_cur:
    print(f"  Unmatched current: {[r['name_norm'] for r in unmatched_cur]}")
if unmatched_old:
    print(f"  Unmatched old: {list(unmatched_old)}")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 4: Computing geometry differences")
print("=" * 60)

fragments = []

for cur_geom, old_geom, nm in matches:
    if not cur_geom.is_valid:
        cur_geom = make_valid(cur_geom)
    if not old_geom.is_valid:
        old_geom = make_valid(old_geom)
    if isinstance(cur_geom, Polygon):
        cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon):
        old_geom = MultiPolygon([old_geom])

    try:
        added = cur_geom.difference(old_geom)
    except Exception:
        added = MultiPolygon()
    try:
        removed = old_geom.difference(cur_geom)
    except Exception:
        removed = MultiPolygon()
    try:
        unchanged = cur_geom.intersection(old_geom)
    except Exception:
        unchanged = MultiPolygon()

    for g, ct in [
        (added, "added_since_2014"),
        (removed, "removed_since_2014"),
        (unchanged, "unchanged"),
    ]:
        if g.is_empty:
            continue
        if isinstance(g, Polygon):
            g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon) or g.is_empty:
            continue
        fragments.append({
            "geometry": g,
            "change_type": ct,
            "district_name": nm,
            "touches_changed": False,
        })

for row in unmatched_cur:
    g = row["geometry"]
    if g.is_empty:
        continue
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon):
        continue
    fragments.append({
        "geometry": g,
        "change_type": "added_since_2014",
        "district_name": row["name_norm"],
        "touches_changed": True,
    })

for on in unmatched_old:
    g = old_by_name[on]
    if g.is_empty:
        continue
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon):
        continue
    fragments.append({
        "geometry": g,
        "change_type": "removed_since_2014",
        "district_name": on,
        "touches_changed": True,
    })

print(f"  Fragment count: {len(fragments)}")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 5: Computing touches_changed")
print("=" * 60)

# Union of all changed areas
changed_polys = [
    f["geometry"]
    for f in fragments
    if f["change_type"] in ("added_since_2014", "removed_since_2014")
]
changed_union = unary_union(changed_polys) if changed_polys else MultiPolygon()
print(f"  Changed area pieces: {len(changed_polys)}")

# Per-district geometry union
dist_geoms = {}
for cur_geom, old_geom, nm in matches:
    dist_geoms.setdefault(nm, []).append(cur_geom)
    dist_geoms.setdefault(nm, []).append(old_geom)
for row in unmatched_cur:
    dist_geoms.setdefault(row["name_norm"], []).append(row["geometry"])
for on in unmatched_old:
    dist_geoms.setdefault(on, []).append(old_by_name[on])

for f in fragments:
    dn = f["district_name"]
    if dn in dist_geoms and not changed_union.is_empty:
        du = unary_union([g for g in dist_geoms[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                touches = du.touches(changed_union)
                f["touches_changed"] = touches or (
                    f["change_type"] in ("added_since_2014", "removed_since_2014")
                )
            except Exception:
                f["touches_changed"] = f["change_type"] in (
                    "added_since_2014", "removed_since_2014"
                )
        else:
            f["touches_changed"] = False
    else:
        f["touches_changed"] = f["change_type"] in (
            "added_since_2014", "removed_since_2014"
        )

tc_true = sum(1 for f in fragments if f["touches_changed"])
print(f"  touches_changed=True: {tc_true} / {len(fragments)}")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 6: Writing GeoJSON")
print("=" * 60)

if not fragments:
    out = gpd.GeoDataFrame(
        {
            "change_type": pd.Series(dtype="str"),
            "district_name": pd.Series(dtype="str"),
            "touches_changed": pd.Series(dtype="bool"),
            "geometry": pd.Series(dtype="object"),
        },
        crs="EPSG:4326",
    )
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326").set_geometry("geometry")
    # Enforce MultiPolygon
    out["geometry"] = out["geometry"].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Written {len(out)} features")

if len(out):
    print(f"  change_type: {out['change_type'].value_counts().to_dict()}")
    print(f"  districts: {sorted(out['district_name'].unique())}")
    print(f"  touches_changed True: {out['touches_changed'].sum()}")

print()
print("Done!")

#!/usr/bin/env python3
"""
Fetch Vienna district boundaries from Overpass by known relation IDs.
Output vienna_boundary_changes.geojson.
"""
import re, sys, time, json
import requests
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid

UA = "GIS-Analyst-Agent/1.0"

def overpass_get(query, timeout=60, endpoint="https://overpass-api.de/api/interpreter",
                 retries=5, delay=3):
    for attempt in range(retries):
        try:
            r = requests.post(endpoint, data={"data": query},
                              headers={"User-Agent": UA}, timeout=timeout)
            if r.status_code == 429:
                print(f"  Rate limited, waiting {delay * (attempt+1)}s...")
                time.sleep(delay * (attempt + 1))
                continue
            if r.status_code == 504:
                print(f"  504 Gateway Timeout, retrying...")
                time.sleep(delay * (attempt + 1))
                continue
            r.raise_for_status()
            return r.json()
        except (requests.ConnectionError, requests.Timeout) as e:
            print(f"  Network error attempt {attempt+1}: {e!r}")
            time.sleep(delay * (attempt + 1))
        except requests.HTTPError as e:
            if e.response is not None and e.response.status_code in (429, 504):
                time.sleep(delay * (attempt + 1))
                continue
            raise
    raise RuntimeError("Overpass failed after all retries")


def geom_from_members(members):
    parts = []
    for m in members:
        if m.get("role") == "outer" and "geometry" in m:
            coords = [(g["lon"], g["lat"]) for g in m["geometry"]]
            if len(coords) >= 3:
                parts.append(Polygon(coords))
    if not parts:
        return None
    merged = unary_union(parts)
    if isinstance(merged, Polygon):
        merged = MultiPolygon([merged])
    if merged.is_empty or not isinstance(merged, MultiPolygon):
        return None
    return merged


def norm_name(name):
    if not isinstance(name, str):
        return ""
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    s = re.sub(r'bezirk$', 'bezirk', s)
    s = re.sub(r'^wien[-\s]', '', s)
    s = re.sub(r'[-\s]wien$', '', s)
    s = s.replace('st.', 'sankt').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s


# ── Known Vienna district relation IDs (from the ids query) ──────────────────
rel_ids = [
    1990590, 1990591, 1990592, 1990593, 1990594, 1990595,
    1990596, 1990597, 1990598, 1990599, 1990600,
    1991416, 1991433, 1991434, 1991435, 1991436, 1991437,
    1991438, 1991439, 1991440, 1991441, 1991442, 1991443,
]

print("=" * 60)
print("STEP 1: Fetching current districts (individual relations)")
print("=" * 60)

current_features = []
for rid in rel_ids:
    q = f"""
    [out:json][timeout:60];
    rel({rid});
    (._;>;);
    out geom;
    """
    try:
        j = overpass_get(q, timeout=60, retries=3, delay=3)
        for el in j.get("elements", []):
            if el["type"] != "relation":
                continue
            tags = el.get("tags", {})
            geom = geom_from_members(el.get("members", []))
            if geom is not None:
                current_features.append({"relation_id": rid, "geometry": geom, **tags})
                print(f"  ✓ rel {rid}: {tags.get('name', '?')}")
            break  # only first relation
    except Exception as e:
        print(f"  ✗ rel {rid}: {e!r}")
    time.sleep(0.5)

print(f"\n  → {len(current_features)} current districts")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 2: Fetching 2014-01-01 districts")
print("=" * 60)

old_features = []
for rid in rel_ids:
    q = f"""
    [out:json][timeout:60][date:"2014-01-01T00:00:00Z"];
    rel({rid});
    out geom;
    """
    try:
        j = overpass_get(q, timeout=60, retries=3, delay=3)
        for el in j.get("elements", []):
            if el["type"] != "relation":
                continue
            tags = el.get("tags", {})
            geom = geom_from_members(el.get("members", []))
            if geom is not None:
                old_features.append({"relation_id": rid, "geometry": geom, **tags})
                print(f"  ✓ old rel {rid}: {tags.get('name', '?')}")
            break
    except Exception as e:
        print(f"  ✗ old rel {rid}: {e!r}")
    time.sleep(0.5)

print(f"\n  → {len(old_features)} historical districts")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 3: Building GeoDataFrames & matching")
print("=" * 60)

for f in current_features:
    f['name'] = f.get('name') or f.get('name:de') or f.get('official_name') or 'unknown'
for f in old_features:
    f['name'] = f.get('name') or f.get('name:de') or f.get('official_name') or 'unknown'

cur_gdf = gpd.GeoDataFrame(current_features, crs="EPSG:4326").dropna(subset=["geometry"])
old_gdf = gpd.GeoDataFrame(old_features, crs="EPSG:4326").dropna(subset=["geometry"])

cur_gdf["geometry"] = cur_gdf["geometry"].apply(lambda g: make_valid(g) if not g.is_valid else g)
old_gdf["geometry"] = old_gdf["geometry"].apply(lambda g: make_valid(g) if not g.is_valid else g)

cur_gdf["name_norm"] = cur_gdf["name"].apply(norm_name)
old_gdf["name_norm"] = old_gdf["name"].apply(norm_name)

print(f"  Current ({len(cur_gdf)}): {sorted(cur_gdf['name_norm'].tolist())}")
print(f"  Old ({len(old_gdf)}): {sorted(old_gdf['name_norm'].tolist())}")

old_by_name = {row["name_norm"]: row["geometry"] for _, row in old_gdf.iterrows()}
matches = []
unmatched_cur = []
unmatched_old = set(old_by_name.keys())

for _, row in cur_gdf.iterrows():
    nm = row["name_norm"]
    if nm in old_by_name:
        matches.append((row["geometry"], old_by_name[nm], nm))
        unmatched_old.discard(nm)
    else:
        unmatched_cur.append(row)

# Fuzzy match by number
for cur_row in list(unmatched_cur):
    nm = cur_row["name_norm"]
    m = re.match(r"^(\d+)", nm)
    if not m:
        continue
    n = m.group(1)
    for on in list(unmatched_old):
        m2 = re.match(r"^(\d+)", on)
        if m2 and m2.group(1) == n:
            matches.append((cur_row["geometry"], old_by_name[on], nm))
            unmatched_cur = [r for r in unmatched_cur if r["name_norm"] != nm]
            unmatched_old.discard(on)
            break

print(f"  Matched: {len(matches)}")
if unmatched_cur: print(f"  Unmatched cur: {[r['name_norm'] for r in unmatched_cur]}")
if unmatched_old: print(f"  Unmatched old: {list(unmatched_old)}")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 4: Computing geometry differences")
print("=" * 60)

fragments = []

for cur_geom, old_geom, nm in matches:
    if not cur_geom.is_valid: cur_geom = make_valid(cur_geom)
    if not old_geom.is_valid: old_geom = make_valid(old_geom)
    if isinstance(cur_geom, Polygon): cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon): old_geom = MultiPolygon([old_geom])

    try:
        added = cur_geom.difference(old_geom)
    except:
        added = MultiPolygon()
    try:
        removed = old_geom.difference(cur_geom)
    except:
        removed = MultiPolygon()
    try:
        unchanged = cur_geom.intersection(old_geom)
    except:
        unchanged = MultiPolygon()

    for g, ct in [('added_since_2014', added), ('removed_since_2014', removed), ('unchanged', unchanged)]:
        g = ct
        ct_type = g
        g = ct  # fix reversed
    # Oops, fixing loop
    for g, ct in [(added, 'added_since_2014'), (removed, 'removed_since_2014'), (unchanged, 'unchanged')]:
        if g.is_empty: continue
        if isinstance(g, Polygon): g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon) or g.is_empty: continue
        fragments.append({"geometry": g, "change_type": ct, "district_name": nm, "touches_changed": False})

for row in unmatched_cur:
    g = row["geometry"]
    if g.is_empty: continue
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon): continue
    fragments.append({"geometry": g, "change_type": "added_since_2014", "district_name": row["name_norm"], "touches_changed": True})

for on in unmatched_old:
    g = old_by_name[on]
    if g.is_empty: continue
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon): continue
    fragments.append({"geometry": g, "change_type": "removed_since_2014", "district_name": on, "touches_changed": True})

print(f"  Fragment count: {len(fragments)}")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 5: touches_changed")
print("=" * 60)

changed_polys = [f["geometry"] for f in fragments if f["change_type"] in ("added_since_2014", "removed_since_2014")]
changed_union = unary_union(changed_polys) if changed_polys else MultiPolygon()
print(f"  Changed pieces: {len(changed_polys)}")

dist_geoms = {}
for cur_geom, old_geom, nm in matches:
    dist_geoms.setdefault(nm, []).extend([cur_geom, old_geom])
for row in unmatched_cur:
    dist_geoms.setdefault(row["name_norm"], []).append(row["geometry"])
for on in unmatched_old:
    dist_geoms.setdefault(on, []).append(old_by_name[on])

for f in fragments:
    dn = f["district_name"]
    if dn in dist_geoms and not changed_union.is_empty:
        du = unary_union([g for g in dist_geoms[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                t = du.touches(changed_union)
                f["touches_changed"] = t or (f["change_type"] in ("added_since_2014", "removed_since_2014"))
            except:
                f["touches_changed"] = f["change_type"] in ("added_since_2014", "removed_since_2014")
        else:
            f["touches_changed"] = False
    else:
        f["touches_changed"] = f["change_type"] in ("added_since_2014", "removed_since_2014")

tc_true = sum(1 for f in fragments if f["touches_changed"])
print(f"  touches_changed=True: {tc_true}/{len(fragments)}")

# ─────────────────────────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 6: Writing GeoJSON")
print("=" * 60)

if not fragments:
    out = gpd.GeoDataFrame({"change_type": pd.Series(dtype="str"),
                            "district_name": pd.Series(dtype="str"),
                            "touches_changed": pd.Series(dtype="bool"),
                            "geometry": pd.Series(dtype="object")}, crs="EPSG:4326")
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326").set_geometry("geometry")
    out["geometry"] = out["geometry"].apply(lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g)

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Written {len(out)} features")
if len(out):
    print(f"  change_type: {out['change_type'].value_counts().to_dict()}")
    print(f"  districts: {sorted(out['district_name'].unique())}")
    print(f"  touches_changed True: {out['touches_changed'].sum()}")
print("\nDone!")

#!/usr/bin/env python3
"""
Fetch Vienna district boundaries from Overpass by known relation IDs.
Output vienna_boundary_changes.geojson.
"""
import re, sys, time, json
import requests
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid

UA = "GIS-Analyst-Agent/1.0"

def overpass_get(query, timeout=60, endpoint="https://overpass-api.de/api/interpreter",
                 retries=5, delay=3):
    for attempt in range(retries):
        try:
            r = requests.post(endpoint, data={"data": query},
                              headers={"User-Agent": UA}, timeout=timeout)
            if r.status_code == 429:
                print(f"  Rate limited, waiting {delay * (attempt+1)}s...")
                time.sleep(delay * (attempt + 1))
                continue
            if r.status_code == 504:
                print(f"  504, retrying...")
                time.sleep(delay * (attempt + 1))
                continue
            r.raise_for_status()
            return r.json()
        except (requests.ConnectionError, requests.Timeout) as e:
            print(f"  Network error attempt {attempt+1}: {e!r}")
            time.sleep(delay * (attempt + 1))
        except requests.HTTPError as e:
            if e.response is not None and e.response.status_code in (429, 504):
                time.sleep(delay * (attempt + 1))
                continue
            raise
    raise RuntimeError("Overpass failed after all retries")


def geom_from_members(members):
    parts = []
    for m in members:
        if m.get("role") == "outer" and "geometry" in m:
            coords = [(g["lon"], g["lat"]) for g in m["geometry"]]
            if len(coords) >= 3:
                parts.append(Polygon(coords))
    if not parts:
        return None
    merged = unary_union(parts)
    if isinstance(merged, Polygon):
        merged = MultiPolygon([merged])
    if merged.is_empty or not isinstance(merged, MultiPolygon):
        return None
    return merged


def norm_name(name):
    if not isinstance(name, str):
        return ""
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    s = re.sub(r'bezirk$', 'bezirk', s)
    s = re.sub(r'^wien[-\s]', '', s)
    s = re.sub(r'[-\s]wien$', '', s)
    s = s.replace('st.', 'sankt').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s


# ── Known Vienna district relation IDs ───────────────────────────────────────
rel_ids = [
    1990590, 1990591, 1990592, 1990593, 1990594, 1990595,
    1990596, 1990597, 1990598, 1990599, 1990600,
    1991416, 1991433, 1991434, 1991435, 1991436, 1991437,
    1991438, 1991439, 1991440, 1991441, 1991442, 1991443,
]

# ── 1. Fetch current districts ───────────────────────────────────────────────
print("=" * 60)
print("STEP 1: Fetching current districts")
print("=" * 60)

current_features = []
for rid in rel_ids:
    q = f"""[out:json][timeout:60];rel({rid});(._;>;);out geom;"""
    try:
        j = overpass_get(q, timeout=60, retries=3, delay=3)
        for el in j.get("elements", []):
            if el["type"] != "relation":
                continue
            tags = el.get("tags", {})
            geom = geom_from_members(el.get("members", []))
            if geom is not None:
                current_features.append({"relation_id": rid, "geometry": geom, **tags})
                print(f"  ✓ rel {rid}: {tags.get('name', '?')}")
            break
    except Exception as e:
        print(f"  ✗ rel {rid}: {e!r}")
    time.sleep(0.3)

print(f"  → {len(current_features)} current districts")

# ── 2. Fetch 2014-01-01 districts ────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 2: Fetching 2014-01-01 districts")
print("=" * 60)

old_features = []
for rid in rel_ids:
    q = f"""[out:json][timeout:60][date:"2014-01-01T00:00:00Z"];rel({rid});out geom;"""
    try:
        j = overpass_get(q, timeout=60, retries=3, delay=3)
        for el in j.get("elements", []):
            if el["type"] != "relation":
                continue
            tags = el.get("tags", {})
            geom = geom_from_members(el.get("members", []))
            if geom is not None:
                old_features.append({"relation_id": rid, "geometry": geom, **tags})
                print(f"  ✓ old rel {rid}: {tags.get('name', '?')}")
            break
    except Exception as e:
        print(f"  ✗ old rel {rid}: {e!r}")
    time.sleep(0.3)

print(f"  → {len(old_features)} historical districts")

# ── 3. Build GeoDataFrames ───────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 3: Building GeoDataFrames & matching")
print("=" * 60)

for f in current_features:
    f['name'] = f.get('name') or f.get('name:de') or f.get('official_name') or 'unknown'
for f in old_features:
    f['name'] = f.get('name') or f.get('name:de') or f.get('official_name') or 'unknown'

cur_gdf = gpd.GeoDataFrame(current_features, crs="EPSG:4326").dropna(subset=["geometry"])
old_gdf = gpd.GeoDataFrame(old_features, crs="EPSG:4326").dropna(subset=["geometry"])

cur_gdf["geometry"] = cur_gdf["geometry"].apply(lambda g: make_valid(g) if not g.is_valid else g)
old_gdf["geometry"] = old_gdf["geometry"].apply(lambda g: make_valid(g) if not g.is_valid else g)

cur_gdf["name_norm"] = cur_gdf["name"].apply(norm_name)
old_gdf["name_norm"] = old_gdf["name"].apply(norm_name)

print(f"  Current ({len(cur_gdf)}): {sorted(cur_gdf['name_norm'].tolist())}")
print(f"  Old ({len(old_gdf)}): {sorted(old_gdf['name_norm'].tolist())}")

old_by_name = {row["name_norm"]: row["geometry"] for _, row in old_gdf.iterrows()}
matches = []
unmatched_cur = []
unmatched_old = set(old_by_name.keys())

for _, row in cur_gdf.iterrows():
    nm = row["name_norm"]
    if nm in old_by_name:
        matches.append((row["geometry"], old_by_name[nm], nm))
        unmatched_old.discard(nm)
    else:
        unmatched_cur.append(row)

# Fuzzy match by number prefix
for cur_row in list(unmatched_cur):
    nm = cur_row["name_norm"]
    m = re.match(r"^(\d+)", nm)
    if not m:
        continue
    n = m.group(1)
    for on in list(unmatched_old):
        m2 = re.match(r"^(\d+)", on)
        if m2 and m2.group(1) == n:
            matches.append((cur_row["geometry"], old_by_name[on], nm))
            unmatched_cur = [r for r in unmatched_cur if r["name_norm"] != nm]
            unmatched_old.discard(on)
            break

print(f"  Matched: {len(matches)}")
if unmatched_cur:
    print(f"  Unmatched current: {[r['name_norm'] for r in unmatched_cur]}")
if unmatched_old:
    print(f"  Unmatched old: {list(unmatched_old)}")

# ── 4. Compute differences ───────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 4: Computing geometry differences")
print("=" * 60)

fragments = []

for cur_geom, old_geom, nm in matches:
    if not cur_geom.is_valid:
        cur_geom = make_valid(cur_geom)
    if not old_geom.is_valid:
        old_geom = make_valid(old_geom)
    if isinstance(cur_geom, Polygon):
        cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon):
        old_geom = MultiPolygon([old_geom])

    try:
        added = cur_geom.difference(old_geom)
    except Exception:
        added = MultiPolygon()
    try:
        removed = old_geom.difference(cur_geom)
    except Exception:
        removed = MultiPolygon()
    try:
        unchanged = cur_geom.intersection(old_geom)
    except Exception:
        unchanged = MultiPolygon()

    for g, ct in [(added, "added_since_2014"), (removed, "removed_since_2014"), (unchanged, "unchanged")]:
        if g.is_empty:
            continue
        if isinstance(g, Polygon):
            g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon) or g.is_empty:
            continue
        fragments.append({"geometry": g, "change_type": ct, "district_name": nm, "touches_changed": False})

for row in unmatched_cur:
    g = row["geometry"]
    if g.is_empty:
        continue
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon):
        continue
    fragments.append({"geometry": g, "change_type": "added_since_2014", "district_name": row["name_norm"], "touches_changed": True})

for on in unmatched_old:
    g = old_by_name[on]
    if g.is_empty:
        continue
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon):
        continue
    fragments.append({"geometry": g, "change_type": "removed_since_2014", "district_name": on, "touches_changed": True})

print(f"  Total fragments: {len(fragments)}")

# ── 5. Compute touches_changed ───────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 5: Computing touches_changed")
print("=" * 60)

changed_polys = [f["geometry"] for f in fragments if f["change_type"] in ("added_since_2014", "removed_since_2014")]
changed_union = unary_union(changed_polys) if changed_polys else MultiPolygon()
print(f"  Changed area pieces: {len(changed_polys)}")

dist_geoms = {}
for cur_geom, old_geom, nm in matches:
    dist_geoms.setdefault(nm, []).extend([cur_geom, old_geom])
for row in unmatched_cur:
    dist_geoms.setdefault(row["name_norm"], []).append(row["geometry"])
for on in unmatched_old:
    dist_geoms.setdefault(on, []).append(old_by_name[on])

for f in fragments:
    dn = f["district_name"]
    if dn in dist_geoms and not changed_union.is_empty:
        du = unary_union([g for g in dist_geoms[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                t = du.touches(changed_union)
                f["touches_changed"] = t or (f["change_type"] in ("added_since_2014", "removed_since_2014"))
            except Exception:
                f["touches_changed"] = f["change_type"] in ("added_since_2014", "removed_since_2014")
        else:
            f["touches_changed"] = False
    else:
        f["touches_changed"] = f["change_type"] in ("added_since_2014", "removed_since_2014")

tc = sum(1 for f in fragments if f["touches_changed"])
print(f"  touches_changed=True: {tc}/{len(fragments)}")

# ── 6. Write output ──────────────────────────────────────────────────────────
print()
print("=" * 60)
print("STEP 6: Writing GeoJSON")
print("=" * 60)

if not fragments:
    out = gpd.GeoDataFrame({
        "change_type": pd.Series(dtype="str"),
        "district_name": pd.Series(dtype="str"),
        "touches_changed": pd.Series(dtype="bool"),
        "geometry": pd.Series(dtype="object"),
    }, crs="EPSG:4326")
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326").set_geometry("geometry")
    out["geometry"] = out["geometry"].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Written {len(out)} features")
if len(out):
    print(f"  change_type: {out['change_type'].value_counts().to_dict()}")
    print(f"  districts: {sorted(out['district_name'].unique())}")
    print(f"  touches_changed=True: {out['touches_changed'].sum()}")

print("\nDone!")

#!/usr/bin/env python3
"""
Build Vienna boundary changes from cached Overpass data.
Output: vienna_boundary_changes.geojson
"""
import re, sys, json, pickle

import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid

# ── Load cached data ─────────────────────────────────────────────────────────
print("Loading cached data...")

with open("/work/_rel_info.pkl", "rb") as f:
    cache = pickle.load(f)
current_info = cache["current"]          # dict: rel_id -> {ways, tags, members}
outer_ways_all = cache["outer_ways"]

with open("/work/_rel_2014.pkl", "rb") as f:
    rel_2014_info = pickle.load(f)       # dict: rel_id -> {ways, tags}

with open("/work/_ways_2014.pkl", "rb") as f:
    ways_2014 = pickle.load(f)           # dict: way_id -> [(lon,lat),...]

print(f"Current relations: {len(current_info)}")
print(f"2014 relations: {len(rel_2014_info)}")
print(f"2014 way geometries: {len(ways_2014)}")

# ── Build current district geometries from relation members ──────────────────
print("\nBuilding current district geometries...")

def build_geom_from_members(members):
    """Build MultiPolygon from relation members list (with geometry)."""
    parts = []
    for m in members:
        if m.get("role") == "outer" and "geometry" in m:
            coords = [(g["lon"], g["lat"]) for g in m["geometry"]]
            if len(coords) >= 3:
                parts.append(Polygon(coords))
    if not parts:
        return None
    merged = unary_union(parts)
    if isinstance(merged, Polygon):
        merged = MultiPolygon([merged])
    if merged.is_empty or not isinstance(merged, MultiPolygon):
        return None
    return merged

current_features = []
for rid, info in current_info.items():
    geom = build_geom_from_members(info.get("members", []))
    if geom is not None:
        current_features.append({
            "relation_id": rid,
            "geometry": geom,
            **info["tags"]
        })
print(f"  Built {len(current_features)} current district geometries")

# ── Build 2014 district geometries from way geometries ───────────────────────
print("\nBuilding 2014 district geometries...")

def build_geom_from_way_ids(way_ids, way_geom_dict):
    """Build MultiPolygon from a list of outer way IDs."""
    parts = []
    for wid in way_ids:
        coords = way_geom_dict.get(wid)
        if coords and len(coords) >= 3:
            parts.append(Polygon(coords))
    if not parts:
        return None
    merged = unary_union(parts)
    if isinstance(merged, Polygon):
        merged = MultiPolygon([merged])
    if merged.is_empty or not isinstance(merged, MultiPolygon):
        return None
    return merged

old_features = []
for rid, info in rel_2014_info.items():
    geom = build_geom_from_way_ids(info["ways"], ways_2014)
    if geom is not None:
        old_features.append({
            "relation_id": rid,
            "geometry": geom,
            **info["tags"]
        })
print(f"  Built {len(old_features)} 2014 district geometries")

# ── Normalise names ──────────────────────────────────────────────────────────
print("\nNormalising district names...")

def norm_name(name):
    if not isinstance(name, str):
        return ""
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    s = re.sub(r'bezirk$', 'bezirk', s)
    s = re.sub(r'^wien[-\s]', '', s)
    s = re.sub(r'[-\s]wien$', '', s)
    s = s.replace('st.', 'sankt').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s

for f in current_features:
    f['name'] = f.get('name') or f.get('name:de') or f.get('official_name') or 'unknown'
for f in old_features:
    f['name'] = f.get('name') or f.get('name:de') or f.get('official_name') or 'unknown'

# Build GeoDataFrames
cur_gdf = gpd.GeoDataFrame(current_features, crs="EPSG:4326").dropna(subset=["geometry"])
old_gdf = gpd.GeoDataFrame(old_features, crs="EPSG:4326").dropna(subset=["geometry"])

cur_gdf["geometry"] = cur_gdf["geometry"].apply(lambda g: make_valid(g) if not g.is_valid else g)
old_gdf["geometry"] = old_gdf["geometry"].apply(lambda g: make_valid(g) if not g.is_valid else g)

cur_gdf["name_norm"] = cur_gdf["name"].apply(norm_name)
old_gdf["name_norm"] = old_gdf["name"].apply(norm_name)

print(f"  Current names ({len(cur_gdf)}): {sorted(cur_gdf['name_norm'].tolist())}")
print(f"  Old names ({len(old_gdf)}): {sorted(old_gdf['name_norm'].tolist())}")

# ── Match districts ──────────────────────────────────────────────────────────
print("\nMatching districts...")

# Filter 2014 to only Vienna districts (admin_level=9 within Vienna)
# The extra 2014 relations are nearby towns that aren't Vienna districts
# We match by name to a current Vienna district
old_by_name = {row["name_norm"]: row["geometry"] for _, row in old_gdf.iterrows()}
matches = []
unmatched_cur = []
unmatched_old = set(old_by_name.keys())

for _, row in cur_gdf.iterrows():
    nm = row["name_norm"]
    if nm in old_by_name:
        matches.append((row["geometry"], old_by_name[nm], nm))
        unmatched_old.discard(nm)
    else:
        unmatched_cur.append(row)

# Fuzzy match by district number
for cur_row in list(unmatched_cur):
    nm = cur_row["name_norm"]
    m = re.match(r"^(\d+)", nm)
    if not m:
        continue
    n = m.group(1)
    for on in list(unmatched_old):
        m2 = re.match(r"^(\d+)", on)
        if m2 and m2.group(1) == n:
            matches.append((cur_row["geometry"], old_by_name[on], nm))
            unmatched_cur = [r for r in unmatched_cur if r["name_norm"] != nm]
            unmatched_old.discard(on)
            break

print(f"  Matched: {len(matches)}")
if unmatched_cur:
    print(f"  Unmatched current: {[r['name_norm'] for r in unmatched_cur]}")
if unmatched_old:
    print(f"  Unmatched old: {list(unmatched_old)}")

# ── Compute geometry differences ─────────────────────────────────────────────
print("\nComputing geometry differences...")

fragments = []

for cur_geom, old_geom, nm in matches:
    if not cur_geom.is_valid: cur_geom = make_valid(cur_geom)
    if not old_geom.is_valid: old_geom = make_valid(old_geom)
    if isinstance(cur_geom, Polygon): cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon): old_geom = MultiPolygon([old_geom])

    try:
        added = cur_geom.difference(old_geom)
    except Exception:
        added = MultiPolygon()
    try:
        removed = old_geom.difference(cur_geom)
    except Exception:
        removed = MultiPolygon()
    try:
        unchanged = cur_geom.intersection(old_geom)
    except Exception:
        unchanged = MultiPolygon()

    for g, ct in [(added, "added_since_2014"), (removed, "removed_since_2014"), (unchanged, "unchanged")]:
        if g.is_empty: continue
        if isinstance(g, Polygon): g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon) or g.is_empty: continue
        fragments.append({"geometry": g, "change_type": ct, "district_name": nm, "touches_changed": False})

for row in unmatched_cur:
    g = row["geometry"]
    if g.is_empty: continue
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon): continue
    fragments.append({"geometry": g, "change_type": "added_since_2014", "district_name": row["name_norm"], "touches_changed": True})

for on in unmatched_old:
    g = old_by_name[on]
    if g.is_empty: continue
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon): continue
    fragments.append({"geometry": g, "change_type": "removed_since_2014", "district_name": on, "touches_changed": True})

print(f"  Total fragments: {len(fragments)}")

# ── Compute touches_changed ──────────────────────────────────────────────────
print("\nComputing touches_changed...")

changed_polys = [f["geometry"] for f in fragments if f["change_type"] in ("added_since_2014", "removed_since_2014")]
changed_union = unary_union(changed_polys) if changed_polys else MultiPolygon()
print(f"  Changed area pieces: {len(changed_polys)}")

# Per-district geometry union
dist_geoms = {}
for cur_geom, old_geom, nm in matches:
    dist_geoms.setdefault(nm, []).extend([cur_geom, old_geom])
for row in unmatched_cur:
    dist_geoms.setdefault(row["name_norm"], []).append(row["geometry"])
for on in unmatched_old:
    dist_geoms.setdefault(on, []).append(old_by_name[on])

for f in fragments:
    dn = f["district_name"]
    if dn in dist_geoms and not changed_union.is_empty:
        du = unary_union([g for g in dist_geoms[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                t = du.touches(changed_union)
                f["touches_changed"] = t or (f["change_type"] in ("added_since_2014", "removed_since_2014"))
            except Exception:
                f["touches_changed"] = f["change_type"] in ("added_since_2014", "removed_since_2014")
        else:
            f["touches_changed"] = False
    else:
        f["touches_changed"] = f["change_type"] in ("added_since_2014", "removed_since_2014")

tc = sum(1 for f in fragments if f["touches_changed"])
print(f"  touches_changed=True: {tc}/{len(fragments)}")

# ── Write output ─────────────────────────────────────────────────────────────
print("\nWriting GeoJSON...")

if not fragments:
    out = gpd.GeoDataFrame({
        "change_type": pd.Series(dtype="str"),
        "district_name": pd.Series(dtype="str"),
        "touches_changed": pd.Series(dtype="bool"),
        "geometry": pd.Series(dtype="object"),
    }, crs="EPSG:4326")
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326").set_geometry("geometry")
    out["geometry"] = out["geometry"].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Written {len(out)} features")
if len(out):
    print(f"  change_type: {out['change_type'].value_counts().to_dict()}")
    print(f"  districts: {sorted(out['district_name'].unique())}")
    print(f"  touches_changed True: {out['touches_changed'].sum()}")

print("\nDone!")

#!/usr/bin/env python3
"""
Build Vienna boundary changes from cached Overpass data.
Output: vienna_boundary_changes.geojson
"""
import re, sys, json, pickle

import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid

# Load cached data
print("Loading cached data...")
with open("/work/_rel_info.pkl", "rb") as f:
    cache = pickle.load(f)
current_info = cache["current"]
with open("/work/_rel_2014.pkl", "rb") as f:
    rel_2014_info = pickle.load(f)
with open("/work/_ways_2014.pkl", "rb") as f:
    ways_2014 = pickle.load(f)

print(f"Current relations: {len(current_info)}")
print(f"2014 relations: {len(rel_2014_info)}")
print(f"2014 way geometries: {len(ways_2014)}")

# ── Helper: build valid MultiPolygon from list of polygon coords ─────────────
def build_multipolygon(polygons):
    """Union a list of Shapely Polygons into a valid MultiPolygon."""
    valid = []
    for p in polygons:
        if p.is_valid and not p.is_empty:
            valid.append(p)
        else:
            bp = make_valid(p)
            if not bp.is_empty:
                valid.append(bp)
    if not valid:
        return None
    try:
        merged = unary_union(valid)
    except Exception:
        merged = unary_union([p.buffer(0) for p in valid])
    if isinstance(merged, Polygon):
        merged = MultiPolygon([merged])
    if merged.is_empty or not isinstance(merged, MultiPolygon):
        return None
    if not merged.is_valid:
        merged = make_valid(merged)
    return merged


# ── Build current district geometries ────────────────────────────────────────
print("\nBuilding current district geometries...")
current_features = []
for rid, info in current_info.items():
    parts = []
    for m in info.get("members", []):
        if m.get("role") == "outer" and "geometry" in m:
            coords = [(g["lon"], g["lat"]) for g in m["geometry"]]
            if len(coords) >= 3:
                parts.append(Polygon(coords))
    geom = build_multipolygon(parts) if parts else None
    if geom is not None:
        current_features.append({
            "relation_id": rid, "geometry": geom, **info["tags"]
        })
print(f"  → {len(current_features)} districts")

# ── Build 2014 district geometries (filter to Vienna only) ───────────────────
print("\nBuilding 2014 district geometries (Vienna-only)...")

# Vienna district relation IDs (the 23 known Vienna districts)
vienna_rel_ids = set(current_info.keys())

old_features = []
for rid, info in rel_2014_info.items():
    # Only keep relations that match a current Vienna district OR are a known Vienna district
    # The extra 2014-only relations (Klein-Engersdorf etc.) are not Vienna districts
    if rid not in vienna_rel_ids:
        nm = info["tags"].get("name", "")
        print(f"  Skipping non-Vienna: rel {rid} '{nm}'")
        continue
    
    parts = []
    for wid in info["ways"]:
        coords = ways_2014.get(wid)
        if coords and len(coords) >= 3:
            parts.append(Polygon(coords))
    geom = build_multipolygon(parts) if parts else None
    if geom is not None:
        old_features.append({
            "relation_id": rid, "geometry": geom, **info["tags"]
        })
print(f"  → {len(old_features)} districts")

# ── Normalise names ──────────────────────────────────────────────────────────
print("\nNormalising names...")

def norm_name(name):
    if not isinstance(name, str):
        return ""
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    s = re.sub(r'bezirk$', 'bezirk', s)
    s = re.sub(r'^wien[-\s]', '', s)
    s = re.sub(r'[-\s]wien$', '', s)
    s = s.replace('st.', 'sankt').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s

for f in current_features:
    f['name'] = f.get('name') or f.get('name:de') or f.get('official_name') or 'unknown'
for f in old_features:
    f['name'] = f.get('name') or f.get('name:de') or f.get('official_name') or 'unknown'

cur_gdf = gpd.GeoDataFrame(current_features, crs="EPSG:4326").dropna(subset=["geometry"])
old_gdf = gpd.GeoDataFrame(old_features, crs="EPSG:4326").dropna(subset=["geometry"])

cur_gdf["name_norm"] = cur_gdf["name"].apply(norm_name)
old_gdf["name_norm"] = old_gdf["name"].apply(norm_name)

print(f"  Current ({len(cur_gdf)}): {sorted(cur_gdf['name_norm'].tolist())}")
print(f"  Old ({len(old_gdf)}): {sorted(old_gdf['name_norm'].tolist())}")

# ── Match districts ──────────────────────────────────────────────────────────
print("\nMatching districts...")

old_by_name = {row["name_norm"]: row["geometry"] for _, row in old_gdf.iterrows()}
matches = []
unmatched_cur = []
unmatched_old_names = set(old_by_name.keys())

for _, row in cur_gdf.iterrows():
    nm = row["name_norm"]
    if nm in old_by_name:
        matches.append((row["geometry"], old_by_name[nm], nm))
        unmatched_old_names.discard(nm)
    else:
        unmatched_cur.append(row)

# Fuzzy by district number
for cur_row in list(unmatched_cur):
    nm = cur_row["name_norm"]
    m = re.match(r"^(\d+)", nm)
    if not m:
        continue
    n = m.group(1)
    for on in list(unmatched_old_names):
        m2 = re.match(r"^(\d+)", on)
        if m2 and m2.group(1) == n:
            matches.append((cur_row["geometry"], old_by_name[on], nm))
            unmatched_cur = [r for r in unmatched_cur if r["name_norm"] != nm]
            unmatched_old_names.discard(on)
            break

print(f"  Matched: {len(matches)}")
if unmatched_cur:
    print(f"  Unmatched current: {[r['name_norm'] for r in unmatched_cur]}")
if unmatched_old_names:
    print(f"  Unmatched old: {list(unmatched_old_names)}")

# ── Compute geometry differences (with small buffer tolerance) ───────────────
print("\nComputing geometry differences (tolerance ~1m)...")

# Use a small buffer to absorb coordinate precision differences
# 1 meter in EPSG:4326 ≈ 0.000009 degrees
TOLERANCE = 0.00001  # ~1 meter

fragments = []

for cur_geom, old_geom, nm in matches:
    if not cur_geom.is_valid: cur_geom = make_valid(cur_geom)
    if not old_geom.is_valid: old_geom = make_valid(old_geom)
    if isinstance(cur_geom, Polygon): cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon): old_geom = MultiPolygon([old_geom])

    # Buffer old slightly, then compute difference
    old_buf = old_geom.buffer(TOLERANCE, join_style=2, mitre_limit=5.0)
    cur_buf = cur_geom.buffer(TOLERANCE, join_style=2, mitre_limit=5.0)

    try:
        added = cur_geom.difference(old_buf)
    except Exception:
        added = MultiPolygon()
    try:
        removed = old_geom.difference(cur_buf)
    except Exception:
        removed = MultiPolygon()
    try:
        unchanged = cur_geom.intersection(old_buf)
    except Exception:
        unchanged = MultiPolygon()

    # Also include any area that's in both with the tolerance
    # If added and removed are tiny, treat as unchanged
    added_area = added.area if not added.is_empty else 0
    removed_area = removed.area if not removed.is_empty else 0
    unchanged_area = unchanged.area if not unchanged.is_empty else 0
    total_area = cur_geom.area

    # If changes are < 0.1% of the district area, treat as unchanged
    if added_area < total_area * 0.001 and removed_area < total_area * 0.001:
        # Boundary didn't meaningfully change
        pass  # Keep as-is, but small fragments may exist
    
    for g, ct in [(added, "added_since_2014"), (removed, "removed_since_2014"), (unchanged, "unchanged")]:
        if g.is_empty:
            continue
        if isinstance(g, Polygon):
            g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon) or g.is_empty:
            continue
        fragments.append({
            "geometry": g,
            "change_type": ct,
            "district_name": nm,
            "touches_changed": False,
        })

for row in unmatched_cur:
    g = row["geometry"]
    if g.is_empty: continue
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon): continue
    fragments.append({"geometry": g, "change_type": "added_since_2014", "district_name": row["name_norm"], "touches_changed": True})

for on in unmatched_old_names:
    g = old_by_name[on]
    if g.is_empty: continue
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon): continue
    fragments.append({"geometry": g, "change_type": "removed_since_2014", "district_name": on, "touches_changed": True})

print(f"  Total fragments: {len(fragments)}")
ct_counts = {}
for f in fragments:
    ct_counts[f["change_type"]] = ct_counts.get(f["change_type"], 0) + 1
print(f"  By type: {ct_counts}")

# ── Compute touches_changed ──────────────────────────────────────────────────
print("\nComputing touches_changed...")

changed_polys = [f["geometry"] for f in fragments if f["change_type"] in ("added_since_2014", "removed_since_2014")]
changed_union = unary_union(changed_polys) if changed_polys else MultiPolygon()
print(f"  Changed area pieces: {len(changed_polys)}")

# Per-district geometry union (combine current + old for each matched district)
dist_geoms = {}
for cur_geom, old_geom, nm in matches:
    dist_geoms.setdefault(nm, []).extend([cur_geom, old_geom])
for row in unmatched_cur:
    dist_geoms.setdefault(row["name_norm"], []).append(row["geometry"])
for on in unmatched_old_names:
    dist_geoms.setdefault(on, []).append(old_by_name[on])

for f in fragments:
    dn = f["district_name"]
    # If this fragment is itself a change, touches_changed = True
    if f["change_type"] in ("added_since_2014", "removed_since_2014"):
        f["touches_changed"] = True
        continue
    
    # For unchanged fragments, check if the district touches any changed area
    if dn in dist_geoms and not changed_union.is_empty:
        du = unary_union([g for g in dist_geoms[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                t = du.touches(changed_union)
                f["touches_changed"] = t
            except Exception:
                f["touches_changed"] = False
        else:
            f["touches_changed"] = False
    else:
        f["touches_changed"] = False

tc = sum(1 for f in fragments if f["touches_changed"])
print(f"  touches_changed=True: {tc}/{len(fragments)}")

# ── Write output ─────────────────────────────────────────────────────────────
print("\nWriting GeoJSON...")

if not fragments:
    out = gpd.GeoDataFrame({
        "change_type": pd.Series(dtype="str"),
        "district_name": pd.Series(dtype="str"),
        "touches_changed": pd.Series(dtype="bool"),
        "geometry": pd.Series(dtype="object"),
    }, crs="EPSG:4326")
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326").set_geometry("geometry")
    out["geometry"] = out["geometry"].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Written {len(out)} features")
if len(out):
    print(f"  change_type: {out['change_type'].value_counts().to_dict()}")
    districts = sorted(out['district_name'].unique())
    print(f"  districts ({len(districts)}): {districts}")
    print(f"  touches_changed True: {out['touches_changed'].sum()}")

print("\nDone!")

#!/usr/bin/env python3
"""
Fetch Vienna district boundaries from Overpass, classify boundary changes,
output vienna_boundary_changes.geojson.

Approach:
1. Get current relation IDs and their outer way IDs
2. Fetch all current ways (full geometry, no bbox)
3. Build current district polygons
4. Fetch 2014 ways and build 2014 district polygons
5. Compute differences and classify
"""
import re, sys, time, json, pickle, os

import requests
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid

UA = "GIS-Analyst-Agent/1.0"
CACHE_DIR = "/work/_cache"

def overpass_query(query, timeout=90, retries=5, delay=5):
    """POST an Overpass QL query, return parsed JSON."""
    endpoint = "https://overpass-api.de/api/interpreter"
    for attempt in range(retries):
        try:
            r = requests.post(endpoint, data={"data": query},
                              headers={"User-Agent": UA}, timeout=timeout)
            if r.status_code == 429:
                print(f"  Rate limited, waiting {delay*(attempt+1)}s...")
                time.sleep(delay * (attempt + 1))
                continue
            if r.status_code == 504:
                print(f"  504, waiting {delay*(attempt+1)}s...")
                time.sleep(delay * (attempt + 1))
                continue
            r.raise_for_status()
            return json.loads(r.text)
        except (requests.ConnectionError, requests.Timeout) as e:
            print(f"  Network err attempt {attempt+1}: {e!r}")
            time.sleep(delay * (attempt + 1))
        except requests.HTTPError as e:
            if e.response is not None and e.response.status_code in (429, 504):
                time.sleep(delay * (attempt + 1))
                continue
            raise
    raise RuntimeError("Overpass failed after all retries")


def build_polygon_from_way_coords(way_id, way_geoms):
    """Return a valid Polygon from way coordinates, or None."""
    coords = way_geoms.get(way_id)
    if not coords or len(coords) < 3:
        return None
    p = Polygon(coords)
    if not p.is_valid:
        p = make_valid(p)
    if p.is_valid and not p.is_empty:
        return p.buffer(0)  # heal any self-intersections
    return None


def build_multipolygon_from_ways(way_ids, way_geoms):
    """Union a list of way IDs into a multi-polygon."""
    parts = []
    for wid in way_ids:
        p = build_polygon_from_way_coords(wid, way_geoms)
        if p is not None:
            parts.append(p)
    if not parts:
        return None
    try:
        merged = unary_union(parts)
    except Exception:
        merged = unary_union([p.buffer(0) for p in parts])
    if isinstance(merged, Polygon):
        merged = MultiPolygon([merged])
    if merged.is_empty or not isinstance(merged, MultiPolygon):
        return None
    return merged


def norm_name(name):
    if not isinstance(name, str):
        return ""
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    s = re.sub(r'bezirk$', 'bezirk', s)
    s = re.sub(r'^wien[-\s]', '', s)
    s = re.sub(r'[-\s]wien$', '', s)
    s = s.replace('st.', 'sankt').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s


# ═══════════════════════════════════════════════════════════════════════════════
# STEP 0: Load or fetch data
# ═══════════════════════════════════════════════════════════════════════════════
os.makedirs(CACHE_DIR, exist_ok=True)

# ── 0a. Get current relation IDs and outer way membership ────────────────────
rel_info_file = f"{CACHE_DIR}/rel_info.pkl"
if os.path.exists(rel_info_file):
    print("Loading cached relation info...")
    with open(rel_info_file, "rb") as f:
        current_rel_info = pickle.load(f)
else:
    print("Fetching current relation IDs...")
    # Get relation IDs using bbox (just for IDs, not geometry)
    q = """[out:json][timeout:30];
    rel["admin_level"="9"]["boundary"="administrative"](48.1,16.1,48.35,16.6);
    out body;"""
    j = overpass_query(q, timeout=30)
    current_rel_info = {}
    for el in j.get("elements", []):
        if el["type"] != "relation":
            continue
        rid = el["id"]
        tags = el.get("tags", {})
        outer_ways = [m["ref"] for m in el.get("members", [])
                      if m["type"] == "way" and m.get("role") == "outer"]
        current_rel_info[rid] = {
            "tags": tags,
            "outer_way_ids": outer_ways,
        }
        print(f"  rel {rid}: {tags.get('name','?')} ({len(outer_ways)} outer ways)")
    with open(rel_info_file, "wb") as f:
        pickle.dump(current_rel_info, f)
    print(f"  Saved {len(current_rel_info)} relations")

print(f"Current relations: {len(current_rel_info)}")

# ── 0b. Fetch current way geometries (full, no bbox) ────────────────────────
cur_ways_file = f"{CACHE_DIR}/cur_ways.pkl"
all_outer_ways = set()
for info in current_rel_info.values():
    all_outer_ways.update(info["outer_way_ids"])

if os.path.exists(cur_ways_file):
    print("Loading cached current way geometries...")
    with open(cur_ways_file, "rb") as f:
        cur_way_geoms = pickle.load(f)
    missing = [w for w in all_outer_ways if w not in cur_way_geoms]
else:
    cur_way_geoms = {}
    missing = list(all_outer_ways)

print(f"Total outer ways: {len(all_outer_ways)}, cached: {len(cur_way_geoms)}, missing: {len(missing)}")

if missing:
    # Fetch in batches of 100
    for i in range(0, len(missing), 100):
        batch = missing[i:i+100]
        ids_str = ",".join(str(w) for w in batch)
        q = f"[out:json][timeout:60];way(id:{ids_str});out geom;"
        print(f"  Fetching ways batch {i//100+1} ({len(batch)} ways)...")
        time.sleep(3)
        j = overpass_query(q, timeout=60)
        for el in j.get("elements", []):
            if el["type"] == "way":
                coords = [(g["lon"], g["lat"]) for g in el.get("geometry", [])]
                cur_way_geoms[el["id"]] = coords
        print(f"    Got {sum(1 for e in j.get('elements',[]) if e['type']=='way')} ways")
    with open(cur_ways_file, "wb") as f:
        pickle.dump(cur_way_geoms, f)
    print(f"  Saved {len(cur_way_geoms)} way geometries")

# ── 0c. Fetch 2014 way geometries ──────────────────────────────────────────
old_ways_file = f"{CACHE_DIR}/old_ways.pkl"
if os.path.exists(old_ways_file):
    print("Loading cached 2014 way geometries...")
    with open(old_ways_file, "rb") as f:
        old_way_geoms = pickle.load(f)
    missing_old = [w for w in all_outer_ways if w not in old_way_geoms]
else:
    old_way_geoms = {}
    missing_old = list(all_outer_ways)

print(f"2014 ways cached: {len(old_way_geoms)}, missing: {len(missing_old)}")

if missing_old:
    for i in range(0, len(missing_old), 100):
        batch = missing_old[i:i+100]
        ids_str = ",".join(str(w) for w in batch)
        q = f"""[out:json][timeout:60][date:"2014-01-01T00:00:00Z"];way(id:{ids_str});out geom;"""
        print(f"  Fetching 2014 ways batch {i//100+1} ({len(batch)} ways)...")
        time.sleep(5)
        try:
            j = overpass_query(q, timeout=60)
            for el in j.get("elements", []):
                if el["type"] == "way":
                    coords = [(g["lon"], g["lat"]) for g in el.get("geometry", [])]
                    old_way_geoms[el["id"]] = coords
            print(f"    Got {sum(1 for e in j.get('elements',[]) if e['type']=='way')} ways")
        except Exception as e:
            print(f"    Failed: {e}")
            # Try each way individually
            for wid in batch:
                time.sleep(2)
                q2 = f"""[out:json][timeout:30][date:"2014-01-01T00:00:00Z"];way({wid});out geom;"""
                try:
                    j2 = overpass_query(q2, timeout=30)
                    for el2 in j2.get("elements", []):
                        if el2["type"] == "way":
                            coords = [(g["lon"], g["lat"]) for g in el2.get("geometry", [])]
                            old_way_geoms[el2["id"]] = coords
                except:
                    pass
    with open(old_ways_file, "wb") as f:
        pickle.dump(old_way_geoms, f)
    print(f"  Saved {len(old_way_geoms)} 2014 way geometries")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 1: Build district geometries
# ═══════════════════════════════════════════════════════════════════════════════
print("\nBuilding district geometries...")

current_features = []
old_features = []

for rid, info in current_rel_info.items():
    tags = info["tags"]
    name = tags.get("name") or tags.get("name:de") or tags.get("official_name") or "unknown"
    outer_ways = info["outer_way_ids"]
    
    cur_geom = build_multipolygon_from_ways(outer_ways, cur_way_geoms)
    if cur_geom is None:
        print(f"  ✗ Current {name}: no valid geometry")
        continue
    
    old_geom = build_multipolygon_from_ways(outer_ways, old_way_geoms)
    if old_geom is None:
        print(f"  ✗ Old {name}: no valid geometry (may not exist in 2014)")
        continue
    
    current_features.append({"relation_id": rid, "geometry": cur_geom, "name": name, **tags})
    old_features.append({"relation_id": rid, "geometry": old_geom, "name": name, **tags})

print(f"  Built {len(current_features)} current and {len(old_features)} old district geometries")

# Check areas
for f in current_features:
    gdf = gpd.GeoDataFrame({"geometry": [f["geometry"]]}, crs="EPSG:4326")
    gdf_m = gdf.to_crs("EPSG:31256")
    area_km2 = gdf_m.geometry.area[0] / 1e6
    print(f"  {f['name']}: {area_km2:.2f} km²")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 2: Build GeoDataFrames & match
# ═══════════════════════════════════════════════════════════════════════════════
print("\nBuilding GeoDataFrames & matching...")

cur_gdf = gpd.GeoDataFrame(current_features, crs="EPSG:4326").dropna(subset=["geometry"])
old_gdf = gpd.GeoDataFrame(old_features, crs="EPSG:4326").dropna(subset=["geometry"])

cur_gdf["name_norm"] = cur_gdf["name"].apply(norm_name)
old_gdf["name_norm"] = old_gdf["name"].apply(norm_name)

print(f"  Current ({len(cur_gdf)}): {sorted(cur_gdf['name_norm'].tolist())}")
print(f"  Old ({len(old_gdf)}): {sorted(old_gdf['name_norm'].tolist())}")

old_by_name = {row["name_norm"]: row["geometry"] for _, row in old_gdf.iterrows()}
matches = []
unmatched_cur = []
unmatched_old_names = set(old_by_name.keys())

for _, row in cur_gdf.iterrows():
    nm = row["name_norm"]
    if nm in old_by_name:
        matches.append((row["geometry"], old_by_name[nm], nm))
        unmatched_old_names.discard(nm)
    else:
        unmatched_cur.append(row)

for cur_row in list(unmatched_cur):
    nm = cur_row["name_norm"]
    m = re.match(r"^(\d+)", nm)
    if not m:
        continue
    n = m.group(1)
    for on in list(unmatched_old_names):
        m2 = re.match(r"^(\d+)", on)
        if m2 and m2.group(1) == n:
            matches.append((cur_row["geometry"], old_by_name[on], nm))
            unmatched_cur = [r for r in unmatched_cur if r["name_norm"] != nm]
            unmatched_old_names.discard(on)
            break

print(f"  Matched: {len(matches)}")
if unmatched_cur:
    print(f"  Unmatched current: {[r['name_norm'] for r in unmatched_cur]}")
if unmatched_old_names:
    print(f"  Unmatched old: {list(unmatched_old_names)}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 3: Compute differences (with small tolerance for coordinate noise)
# ═══════════════════════════════════════════════════════════════════════════════
print("\nComputing geometry differences...")

TOL = 0.00001  # ~1m tolerance

fragments = []

for cur_geom, old_geom, nm in matches:
    if not cur_geom.is_valid:
        cur_geom = make_valid(cur_geom).buffer(0)
    if not old_geom.is_valid:
        old_geom = make_valid(old_geom).buffer(0)
    if isinstance(cur_geom, Polygon):
        cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon):
        old_geom = MultiPolygon([old_geom])

    # Buffer to absorb coordinate noise
    old_buf = old_geom.buffer(TOL, join_style=2)
    cur_buf = cur_geom.buffer(TOL, join_style=2)

    try:
        added = cur_geom.difference(old_buf)
    except Exception:
        added = MultiPolygon()
    try:
        removed = old_geom.difference(cur_buf)
    except Exception:
        removed = MultiPolygon()
    try:
        unchanged = cur_geom.intersection(old_buf)
    except Exception:
        unchanged = MultiPolygon()

    for g, ct in [
        (added, "added_since_2014"),
        (removed, "removed_since_2014"),
        (unchanged, "unchanged"),
    ]:
        if g.is_empty:
            continue
        if isinstance(g, Polygon):
            g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon) or g.is_empty:
            continue
        fragments.append({
            "geometry": g,
            "change_type": ct,
            "district_name": nm,
            "touches_changed": False,
        })

for row in unmatched_cur:
    g = row["geometry"]
    if g.is_empty or not isinstance(g, (Polygon, MultiPolygon)):
        continue
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    fragments.append({
        "geometry": g,
        "change_type": "added_since_2014",
        "district_name": row["name_norm"],
        "touches_changed": True,
    })

for on in unmatched_old_names:
    g = old_by_name[on]
    if g.is_empty or not isinstance(g, (Polygon, MultiPolygon)):
        continue
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    fragments.append({
        "geometry": g,
        "change_type": "removed_since_2014",
        "district_name": on,
        "touches_changed": True,
    })

ct_counts = {}
for f in fragments:
    ct_counts[f["change_type"]] = ct_counts.get(f["change_type"], 0) + 1
print(f"  Fragments: {len(fragments)} types: {ct_counts}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 4: Compute touches_changed
# ═══════════════════════════════════════════════════════════════════════════════
print("\nComputing touches_changed...")

changed_polys = [f["geometry"] for f in fragments if f["change_type"] in ("added_since_2014", "removed_since_2014")]
changed_union = unary_union(changed_polys) if changed_polys else MultiPolygon()
print(f"  Changed pieces: {len(changed_polys)}")

# Per-district geometry union (current + old)
dist_geoms = {}
for cur_geom, old_geom, nm in matches:
    dist_geoms.setdefault(nm, []).extend([cur_geom, old_geom])
for row in unmatched_cur:
    dist_geoms.setdefault(row["name_norm"], []).append(row["geometry"])
for on in unmatched_old_names:
    dist_geoms.setdefault(on, []).append(old_by_name[on])

for f in fragments:
    dn = f["district_name"]
    if f["change_type"] in ("added_since_2014", "removed_since_2014"):
        f["touches_changed"] = True
        continue
    # Unchanged: check if district touches changed areas
    if dn in dist_geoms and not changed_union.is_empty:
        du = unary_union([g for g in dist_geoms[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                f["touches_changed"] = du.touches(changed_union)
            except Exception:
                f["touches_changed"] = False
        else:
            f["touches_changed"] = False
    else:
        f["touches_changed"] = False

tc = sum(1 for f in fragments if f["touches_changed"])
print(f"  touches_changed=True: {tc}/{len(fragments)}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 5: Write output
# ═══════════════════════════════════════════════════════════════════════════════
print("\nWriting GeoJSON...")

if not fragments:
    out = gpd.GeoDataFrame({
        "change_type": pd.Series(dtype="str"),
        "district_name": pd.Series(dtype="str"),
        "touches_changed": pd.Series(dtype="bool"),
        "geometry": pd.Series(dtype="object"),
    }, crs="EPSG:4326")
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326").set_geometry("geometry")
    out["geometry"] = out["geometry"].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Written {len(out)} features")
if len(out):
    print(f"  change_type: {out['change_type'].value_counts().to_dict()}")
    print(f"  districts: {sorted(out['district_name'].unique())}")
    print(f"  touches_changed True: {out['touches_changed'].sum()}")

print("\nDone!")

#!/usr/bin/env python3
"""
Fetch Vienna district boundary changes and output vienna_boundary_changes.geojson.

Strategy: use rel(id);out geom; (no bbox, no recurse) for both current and 2014.
Member objects in the response contain full boundary geometries.
"""
import re, sys, time, json, pickle, os

import requests
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid

UA = "GIS-Analyst-Agent/1.0"
CACHE = "/work/_cache2"

def op_query(query, timeout=60, label="", retries=5):
    for att in range(retries):
        try:
            r = requests.post("https://overpass-api.de/api/interpreter",
                              data={"data": query},
                              headers={"User-Agent": UA}, timeout=timeout)
            if r.status_code in (429, 504):
                w = 5 * (att + 1)
                if label: print(f"  {label}: {r.status_code}, wait {w}s")
                time.sleep(w)
                continue
            r.raise_for_status()
            return json.loads(r.text)
        except (requests.ConnectionError, requests.Timeout) as e:
            if label: print(f"  {label}: net err {e!r}")
            time.sleep(5 * (att + 1))
        except requests.HTTPError as e:
            if e.response is not None and e.response.status_code in (429, 504):
                time.sleep(5 * (att + 1))
                continue
            raise
    raise RuntimeError(f"Overpass failed for {label}")

def norm_name(name):
    if not isinstance(name, str): return ""
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    s = re.sub(r'^wien[- ]', '', s)
    s = re.sub(r'[- ]wien$', '', s)
    s = s.replace('st.', 'sankt ').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor ')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s

def members_to_mpoly(members):
    """Extract Polygon from member objects that have .geometry (out geom)."""
    parts = []
    for m in members:
        if m.get("role") == "outer" and "geometry" in m:
            coords = [(g["lon"], g["lat"]) for g in m["geometry"]]
            if len(coords) >= 3:
                p = Polygon(coords)
                if not p.is_valid:
                    p = make_valid(p)
                if p.is_valid and not p.is_empty:
                    parts.append(p.buffer(0))
    if not parts:
        return None
    try:
        u = unary_union(parts)
    except:
        u = unary_union([p.buffer(0) for p in parts])
    if isinstance(u, Polygon):
        u = MultiPolygon([u])
    if u.is_empty or not isinstance(u, MultiPolygon):
        return None
    return u

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 0: Get relation IDs
# ═══════════════════════════════════════════════════════════════════════════════
os.makedirs(CACHE, exist_ok=True)

ids_file = f"{CACHE}/rel_ids.json"
if os.path.exists(ids_file):
    with open(ids_file) as f:
        vienna_rel_ids = json.load(f)
    print(f"Loaded {len(vienna_rel_ids)} relation IDs from cache")
else:
    print("Fetching Vienna district relation IDs...")
    q = """[out:json][timeout:30];
    rel["admin_level"="9"]["boundary"="administrative"](48.1,16.1,48.35,16.6);
    out body;"""
    j = op_query(q, timeout=30)
    vienna_rel_ids = []
    for el in j.get("elements", []):
        if el["type"] == "relation":
            vienna_rel_ids.append(el["id"])
    print(f"  Found {len(vienna_rel_ids)} relation IDs")
    with open(ids_file, "w") as f:
        json.dump(vienna_rel_ids, f)

print(f"Relation IDs: {vienna_rel_ids}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 1: Fetch current districts (full geom, no bbox)
# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- STEP 1: Fetching current districts ---")

cur_file = f"{CACHE}/current.pkl"
if os.path.exists(cur_file):
    print("Loading from cache...")
    with open(cur_file, "rb") as f:
        current_features = pickle.load(f)
else:
    current_features = []
    for rid in vienna_rel_ids:
        q = f"[out:json][timeout:30];rel({rid});out geom;"
        try:
            j = op_query(q, timeout=30, label=f"cur {rid}")
            for el in j.get("elements", []):
                if el["type"] == "relation":
                    tags = el.get("tags", {})
                    geom = members_to_mpoly(el.get("members", []))
                    if geom is not None:
                        name = tags.get("name") or tags.get("name:de") or "unknown"
                        current_features.append({
                            "relation_id": rid, "geometry": geom, "name": name, **tags
                        })
                        print(f"  ✓ {rid}: {name}")
                    break
        except Exception as e:
            print(f"  ✗ {rid}: {e!r}")
        time.sleep(0.2)
    with open(cur_file, "wb") as f:
        pickle.dump(current_features, f)
    print(f"  Saved {len(current_features)} current features")

print(f"Current districts: {len(current_features)}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 2: Fetch 2014 districts
# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- STEP 2: Fetching 2014-01-01 districts ---")

old_file = f"{CACHE}/old.pkl"
if os.path.exists(old_file):
    print("Loading from cache...")
    with open(old_file, "rb") as f:
        old_features = pickle.load(f)
else:
    old_features = []
    for rid in vienna_rel_ids:
        q = f"""[out:json][timeout:30][date:"2014-01-01T00:00:00Z"];rel({rid});out geom;"""
        try:
            j = op_query(q, timeout=30, label=f"old {rid}")
            for el in j.get("elements", []):
                if el["type"] == "relation":
                    tags = el.get("tags", {})
                    geom = members_to_mpoly(el.get("members", []))
                    if geom is not None:
                        name = tags.get("name") or tags.get("name:de") or "unknown"
                        old_features.append({
                            "relation_id": rid, "geometry": geom, "name": name, **tags
                        })
                        print(f"  ✓ {rid}: {name}")
                    break
        except Exception as e:
            print(f"  ✗ {rid}: {e!r}")
        time.sleep(0.5)
    with open(old_file, "wb") as f:
        pickle.dump(old_features, f)
    print(f"  Saved {len(old_features)} old features")

print(f"Old districts: {len(old_features)}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 3: Build & match
# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- STEP 3: Matching districts ---")

cur_gdf = gpd.GeoDataFrame(current_features, crs="EPSG:4326").dropna(subset=["geometry"])
old_gdf = gpd.GeoDataFrame(old_features, crs="EPSG:4326").dropna(subset=["geometry"])

cur_gdf["name_norm"] = cur_gdf["name"].apply(norm_name)
old_gdf["name_norm"] = old_gdf["name"].apply(norm_name)

print(f"  Current ({len(cur_gdf)}): {sorted(cur_gdf['name_norm'].tolist())}")
print(f"  Old ({len(old_gdf)}): {sorted(old_gdf['name_norm'].tolist())}")

# Verify areas
for _, r in cur_gdf.iterrows():
    g = gpd.GeoDataFrame({"geometry": [r["geometry"]]}, crs="EPSG:4326")
    gm = g.to_crs("EPSG:31256")
    a = gm.geometry.area[0] / 1e6
    print(f"  {r['name']}: {a:.2f} km²")

old_by_name = {row["name_norm"]: row["geometry"] for _, row in old_gdf.iterrows()}
matches = []
unmatched_cur = []
unmatched_old = set(old_by_name.keys())

for _, row in cur_gdf.iterrows():
    nm = row["name_norm"]
    if nm in old_by_name:
        matches.append((row["geometry"], old_by_name[nm], nm))
        unmatched_old.discard(nm)
    else:
        unmatched_cur.append(row)

for cr in list(unmatched_cur):
    nm = cr["name_norm"]
    m = re.match(r"^(\d+)", nm)
    if m:
        n = m.group(1)
        for on in list(unmatched_old):
            m2 = re.match(r"^(\d+)", on)
            if m2 and m2.group(1) == n:
                matches.append((cr["geometry"], old_by_name[on], nm))
                unmatched_cur = [r for r in unmatched_cur if r["name_norm"] != nm]
                unmatched_old.discard(on)
                break

print(f"  Matched: {len(matches)}")
if unmatched_cur:
    print(f"  Unmatched cur: {[r['name_norm'] for r in unmatched_cur]}")
if unmatched_old:
    print(f"  Unmatched old: {list(unmatched_old)}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 4: Differences
# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- STEP 4: Computing differences ---")

TOL = 0.00001  # ~1m
fragments = []

for cur_geom, old_geom, nm in matches:
    if isinstance(cur_geom, Polygon):
        cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon):
        old_geom = MultiPolygon([old_geom])

    old_buf = old_geom.buffer(TOL, join_style=2)
    cur_buf = cur_geom.buffer(TOL, join_style=2)

    try:
        added = cur_geom.difference(old_buf)
    except:
        added = MultiPolygon()
    try:
        removed = old_geom.difference(cur_buf)
    except:
        removed = MultiPolygon()
    try:
        unchanged = cur_geom.intersection(old_buf)
    except:
        unchanged = MultiPolygon()

    for g, ct in [(added, "added_since_2014"), (removed, "removed_since_2014"), (unchanged, "unchanged")]:
        if g.is_empty:
            continue
        if isinstance(g, Polygon):
            g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon):
            continue
        fragments.append({"geometry": g, "change_type": ct, "district_name": nm, "touches_changed": False})

for cr in unmatched_cur:
    g = cr["geometry"]
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon):
        continue
    fragments.append({"geometry": g, "change_type": "added_since_2014", "district_name": cr["name_norm"], "touches_changed": True})

for on in unmatched_old:
    g = old_by_name[on]
    if isinstance(g, Polygon):
        g = MultiPolygon([g])
    if not isinstance(g, MultiPolygon):
        continue
    fragments.append({"geometry": g, "change_type": "removed_since_2014", "district_name": on, "touches_changed": True})

print(f"  Total fragments: {len(fragments)}")
ctc = {}
for f in fragments:
    ctc[f["change_type"]] = ctc.get(f["change_type"], 0) + 1
print(f"  Types: {ctc}")

# Area check
for f in fragments:
    gdf_a = gpd.GeoDataFrame({"geometry": [f["geometry"]]}, crs="EPSG:4326")
    gdf_am = gdf_a.to_crs("EPSG:31256")
    a_km2 = gdf_am.geometry.area[0] / 1e6
    print(f"  {f['district_name']} {f['change_type']}: {a_km2:.6f} km²")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 5: touches_changed
# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- STEP 5: touches_changed ---")

changed_polys = [f["geometry"] for f in fragments if f["change_type"] in ("added_since_2014", "removed_since_2014")]
changed_union = unary_union(changed_polys) if changed_polys else MultiPolygon()

dist_geoms = {}
for cur_geom, old_geom, nm in matches:
    dist_geoms.setdefault(nm, []).extend([cur_geom, old_geom])
for cr in unmatched_cur:
    dist_geoms.setdefault(cr["name_norm"], []).append(cr["geometry"])
for on in unmatched_old:
    dist_geoms.setdefault(on, []).append(old_by_name[on])

for f in fragments:
    dn = f["district_name"]
    if f["change_type"] in ("added_since_2014", "removed_since_2014"):
        f["touches_changed"] = True
        continue
    if dn in dist_geoms and not changed_union.is_empty:
        du = unary_union([g for g in dist_geoms[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                f["touches_changed"] = du.touches(changed_union)
            except:
                f["touches_changed"] = False
        else:
            f["touches_changed"] = False
    else:
        f["touches_changed"] = False

tc = sum(1 for f in fragments if f["touches_changed"])
print(f"  touches_changed=True: {tc}/{len(fragments)}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 6: Output
# ═══════════════════════════════════════════════════════════════════════════════
print("\n--- STEP 6: Writing GeoJSON ---")

if not fragments:
    out = gpd.GeoDataFrame({
        "change_type": pd.Series(dtype="str"),
        "district_name": pd.Series(dtype="str"),
        "touches_changed": pd.Series(dtype="bool"),
        "geometry": pd.Series(dtype="object"),
    }, crs="EPSG:4326")
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326").set_geometry("geometry")
    out["geometry"] = out["geometry"].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Written {len(out)} features")
if len(out):
    print(f"  change_type: {out['change_type'].value_counts().to_dict()}")
    print(f"  districts: {sorted(out['district_name'].unique())}")
    print(f"  touches_changed True: {out['touches_changed'].sum()}")

print("\nDone!")

#!/usr/bin/env python3
"""
Vienna district boundary changes pipeline.
Fast approach: fetch all outer ways in bulk, then build geometries.
"""
import re, sys, time, json, pickle, os
import requests
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid

UA = "GIS-Analyst-Agent/1.0"
CACHE = "/work/_cache3"
os.makedirs(CACHE, exist_ok=True)

def op_q(query, timeout=60, retries=5):
    for att in range(retries):
        try:
            r = requests.post("https://overpass-api.de/api/interpreter",
                              data={"data": query}, headers={"User-Agent": UA},
                              timeout=timeout)
            if r.status_code in (429, 504):
                time.sleep(6 * (att + 1))
                continue
            r.raise_for_status()
            return json.loads(r.text)
        except (requests.ConnectionError, requests.Timeout):
            time.sleep(6 * (att + 1))
        except requests.HTTPError as e:
            if e.response is not None and e.response.status_code in (429, 504):
                time.sleep(6 * (att + 1))
                continue
            raise
    raise RuntimeError("Overpass failed")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 1: Get relation outer-way membership (body, no geom needed)
# ═══════════════════════════════════════════════════════════════════════════════
print("STEP 1: Fetching relation way membership...")
f_rel = f"{CACHE}/rel_members.pkl"
if os.path.exists(f_rel):
    with open(f_rel, "rb") as f:
        rids, rel_info, all_ways = pickle.load(f)
    print(f"  Cached: {len(rids)} relations, {len(all_ways)} outer ways")
else:
    q = """[out:json][timeout:30];
    rel["admin_level"="9"]["boundary"="administrative"](48.1,16.1,48.35,16.6);
    out body;"""
    j = op_q(q, timeout=30)
    rids = []
    rel_info = {}
    all_ways = set()
    for el in j.get("elements", []):
        if el["type"] != "relation":
            continue
        rids.append(el["id"])
        tags = el.get("tags", {})
        outer = [m["ref"] for m in el.get("members", [])
                 if m["type"] == "way" and m.get("role") == "outer"]
        # Also store members list for building geometries from relation out geom
        rel_info[el["id"]] = {"tags": tags, "outer_ways": outer}
        all_ways.update(outer)
    all_ways = list(all_ways)
    with open(f_rel, "wb") as f:
        pickle.dump((rids, rel_info, all_ways), f)
    print(f"  {len(rids)} relations, {len(all_ways)} outer ways")

print(f"  Relations: {[rel_info[r]['tags'].get('name','?') for r in rids]}")
print(f"  Total outer ways: {len(all_ways)}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 2: Fetch current way geometries in bulk
# ═══════════════════════════════════════════════════════════════════════════════
print("\nSTEP 2: Fetching current way geometries...")
f_wcur = f"{CACHE}/ways_cur.pkl"
if os.path.exists(f_wcur):
    with open(f_wcur, "rb") as f:
        ways_cur = pickle.load(f)
    print(f"  Cached: {len(ways_cur)} ways")
else:
    ways_cur = {}
    # Split into batches of 200
    for i in range(0, len(all_ways), 200):
        batch = all_ways[i:i+200]
        ids = ",".join(str(w) for w in batch)
        q = f"[out:json][timeout:60];way(id:{ids});out geom;"
        print(f"  Batch {i//200+1}: {len(batch)} ways...")
        time.sleep(2)
        j = op_q(q, timeout=60)
        for el in j.get("elements", []):
            if el["type"] == "way":
                geom = [(g["lon"], g["lat"]) for g in el.get("geometry", [])]
                ways_cur[el["id"]] = geom
        print(f"    → got {len(j.get('elements',[]))} ways")
    with open(f_wcur, "wb") as f:
        pickle.dump(ways_cur, f)
    print(f"  Saved {len(ways_cur)} ways")

print(f"  Have {len(ways_cur)} / {len(all_ways)} ways")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 3: Fetch 2014 way geometries in bulk
# ═══════════════════════════════════════════════════════════════════════════════
print("\nSTEP 3: Fetching 2014 way geometries...")
f_wold = f"{CACHE}/ways_old.pkl"
if os.path.exists(f_wold):
    with open(f_wold, "rb") as f:
        ways_old = pickle.load(f)
    print(f"  Cached: {len(ways_old)} ways")
else:
    ways_old = {}
    for i in range(0, len(all_ways), 150):  # smaller batches for slower 2014 queries
        batch = all_ways[i:i+150]
        ids = ",".join(str(w) for w in batch)
        q = f"""[out:json][timeout:90][date:"2014-01-01T00:00:00Z"];way(id:{ids});out geom;"""
        print(f"  Batch {i//150+1}: {len(batch)} ways...")
        time.sleep(4)
        try:
            j = op_q(q, timeout=90)
            for el in j.get("elements", []):
                if el["type"] == "way":
                    geom = [(g["lon"], g["lat"]) for g in el.get("geometry", [])]
                    ways_old[el["id"]] = geom
            print(f"    → got {sum(1 for e in j.get('elements',[]) if e['type']=='way')} ways")
        except Exception as e:
            print(f"    → batch failed: {e}, trying individually...")
            for wid in batch:
                time.sleep(2)
                try:
                    q2 = f"""[out:json][timeout:30][date:"2014-01-01T00:00:00Z"];way({wid});out geom;"""
                    j2 = op_q(q2, timeout=30)
                    for el2 in j2.get("elements", []):
                        if el2["type"] == "way":
                            geom = [(g["lon"], g["lat"]) for g in el2.get("geometry", [])]
                            ways_old[el2["id"]] = geom
                except:
                    pass
    with open(f_wold, "wb") as f:
        pickle.dump(ways_old, f)
    print(f"  Saved {len(ways_old)} ways")

print(f"  Have {len(ways_old)} / {len(all_ways)} ways")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 4: Build district geometries
# ═══════════════════════════════════════════════════════════════════════════════
print("\nSTEP 4: Building district geometries...")

def build_district(way_ids, way_dict):
    parts = []
    for wid in way_ids:
        coords = way_dict.get(wid)
        if coords and len(coords) >= 3:
            p = Polygon(coords)
            if not p.is_valid:
                p = make_valid(p)
            if p.is_valid and not p.is_empty:
                parts.append(p.buffer(0))
    if not parts:
        return None
    try:
        u = unary_union(parts)
    except:
        u = unary_union([p.buffer(0) for p in parts])
    if isinstance(u, Polygon):
        u = MultiPolygon([u])
    if u.is_empty or not isinstance(u, MultiPolygon):
        return None
    return u

cur_feats = []
old_feats = []
for rid in rids:
    info = rel_info[rid]
    tags = info["tags"]
    name = tags.get("name") or tags.get("name:de") or "unknown"
    ways = info["outer_ways"]
    
    cg = build_district(ways, ways_cur)
    og = build_district(ways, ways_old)
    
    if cg is not None:
        cur_feats.append({"relation_id": rid, "geometry": cg, "name": name})
    if og is not None:
        old_feats.append({"relation_id": rid, "geometry": og, "name": name})
    
    c_area = 0
    if cg:
        gc = gpd.GeoDataFrame({"geometry": [cg]}, crs="EPSG:4326").to_crs("EPSG:31256")
        c_area = gc.geometry.area[0] / 1e6
    o_area = 0
    if og:
        go = gpd.GeoDataFrame({"geometry": [og]}, crs="EPSG:4326").to_crs("EPSG:31256")
        o_area = go.geometry.area[0] / 1e6
    print(f"  {name}: current={c_area:.2f} km², 2014={o_area:.2f} km²")

print(f"  Built: {len(cur_feats)} current, {len(old_feats)} old")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 5: Match and diff
# ═══════════════════════════════════════════════════════════════════════════════
print("\nSTEP 5: Matching & computing differences...")

def norm_name(name):
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    s = re.sub(r'^wien[- ]', '', s)
    s = re.sub(r'[- ]wien$', '', s)
    s = s.replace('st.', 'sankt ').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor ')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s

for f in cur_feats:
    f["name_norm"] = norm_name(f["name"])
for f in old_feats:
    f["name_norm"] = norm_name(f["name"])

obn = {f["name_norm"]: f["geometry"] for f in old_feats}
matches = []
unmatched_c = []
unmatched_o = set(obn.keys())

for f in cur_feats:
    nm = f["name_norm"]
    if nm in obn:
        matches.append((f["geometry"], obn[nm], nm))
        unmatched_o.discard(nm)
    else:
        unmatched_c.append(f)

for f in list(unmatched_c):
    nm = f["name_norm"]
    m = re.match(r"^(\d+)", nm)
    if m:
        n = m.group(1)
        for on in list(unmatched_o):
            m2 = re.match(r"^(\d+)", on)
            if m2 and m2.group(1) == n:
                matches.append((f["geometry"], obn[on], nm))
                unmatched_c = [x for x in unmatched_c if x["name_norm"] != nm]
                unmatched_o.discard(on)
                break

print(f"  Matched: {len(matches)}")
if unmatched_c: print(f"  Unmatched cur: {[x['name_norm'] for x in unmatched_c]}")
if unmatched_o: print(f"  Unmatched old: {list(unmatched_o)}")

TOL = 0.00001
fragments = []
for cur_geom, old_geom, nm in matches:
    if isinstance(cur_geom, Polygon): cur_geom = MultiPolygon([cur_geom])
    if isinstance(old_geom, Polygon): old_geom = MultiPolygon([old_geom])
    ob = old_geom.buffer(TOL, join_style=2)
    cb = cur_geom.buffer(TOL, join_style=2)
    try:
        a = cur_geom.difference(ob)
    except:
        a = MultiPolygon()
    try:
        r = old_geom.difference(cb)
    except:
        r = MultiPolygon()
    try:
        u = cur_geom.intersection(ob)
    except:
        u = MultiPolygon()
    for g, ct in [(a, "added_since_2014"), (r, "removed_since_2014"), (u, "unchanged")]:
        if g.is_empty: continue
        if isinstance(g, Polygon): g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon): continue
        fragments.append({"geometry": g, "change_type": ct, "district_name": nm, "touches_changed": False})

for f in unmatched_c:
    g = f["geometry"]
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if isinstance(g, MultiPolygon):
        fragments.append({"geometry": g, "change_type": "added_since_2014", "district_name": f["name_norm"], "touches_changed": True})

for on in unmatched_o:
    g = obn[on]
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if isinstance(g, MultiPolygon):
        fragments.append({"geometry": g, "change_type": "removed_since_2014", "district_name": on, "touches_changed": True})

ctc = {}
for f in fragments:
    ctc[f["change_type"]] = ctc.get(f["change_type"], 0) + 1
print(f"  Fragments: {len(fragments)} types: {ctc}")

# Area summary
for f in fragments:
    ga = gpd.GeoDataFrame({"geometry": [f["geometry"]]}, crs="EPSG:4326").to_crs("EPSG:31256")
    ak = ga.geometry.area[0] / 1e6
    if f["change_type"] != "unchanged":
        print(f"  {f['district_name']} {f['change_type']}: {ak:.6f} km²")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 6: touches_changed
# ═══════════════════════════════════════════════════════════════════════════════
print("\nSTEP 6: touches_changed...")

cp = [f["geometry"] for f in fragments if f["change_type"] in ("added_since_2014", "removed_since_2014")]
cu = unary_union(cp) if cp else MultiPolygon()

dg = {}
for cg, og, nm in matches:
    dg.setdefault(nm, []).extend([cg, og])
for f in unmatched_c:
    dg.setdefault(f["name_norm"], []).append(f["geometry"])
for on in unmatched_o:
    dg.setdefault(on, []).append(obn[on])

for f in fragments:
    if f["change_type"] in ("added_since_2014", "removed_since_2014"):
        f["touches_changed"] = True
        continue
    dn = f["district_name"]
    if dn in dg and not cu.is_empty:
        du = unary_union([g for g in dg[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                f["touches_changed"] = du.touches(cu)
            except:
                f["touches_changed"] = False
        else:
            f["touches_changed"] = False
    else:
        f["touches_changed"] = False

tc = sum(1 for f in fragments if f["touches_changed"])
print(f"  touches_changed=True: {tc}/{len(fragments)}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 7: Write output
# ═══════════════════════════════════════════════════════════════════════════════
print("\nSTEP 7: Writing GeoJSON...")

if not fragments:
    out = gpd.GeoDataFrame({"change_type": pd.Series(dtype="str"),
                            "district_name": pd.Series(dtype="str"),
                            "touches_changed": pd.Series(dtype="bool"),
                            "geometry": pd.Series(dtype="object")}, crs="EPSG:4326")
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326").set_geometry("geometry")
    out["geometry"] = out["geometry"].apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g
    )

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Written {len(out)} features")
if len(out):
    print(f"  change_type: {out['change_type'].value_counts().to_dict()}")
    print(f"  districts: {sorted(out['district_name'].unique())}")
    print(f"  touches_changed True: {out['touches_changed'].sum()}")

print("\nDone!")

#!/usr/bin/env python3
"""
Vienna boundary change detection.
Uses rel(id);out geom; (no bbox, verified correct) to get full district geometries.
"""
import re, sys, time, json, pickle, os
import requests
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid

UA = "GIS-Analyst-Agent/1.0"
CACHE = "/work/_cache_final"
os.makedirs(CACHE, exist_ok=True)

def op_q(query, timeout=60, retries=5):
    for att in range(retries):
        try:
            r = requests.post("https://overpass-api.de/api/interpreter",
                              data={"data": query}, headers={"User-Agent": UA},
                              timeout=timeout)
            if r.status_code in (429, 504):
                time.sleep(6 * (att + 1))
                continue
            r.raise_for_status()
            return json.loads(r.text)
        except (requests.ConnectionError, requests.Timeout):
            time.sleep(6 * (att + 1))
        except requests.HTTPError as e:
            if e.response is not None and e.response.status_code in (429, 504):
                time.sleep(6 * (att + 1))
                continue
            raise
    raise RuntimeError("Overpass failed")

def members_to_mpoly(members):
    """Build MultiPolygon from relation member objects (out geom format)."""
    parts = []
    for m in members:
        if m.get("role") == "outer" and "geometry" in m:
            coords = [(g["lon"], g["lat"]) for g in m["geometry"]]
            if len(coords) >= 3:
                p = Polygon(coords)
                if not p.is_valid:
                    p = make_valid(p)
                if p.is_valid and not p.is_empty:
                    parts.append(p.buffer(0))
    if not parts:
        return None
    try:
        u = unary_union(parts)
    except:
        u = unary_union([p.buffer(0) for p in parts])
    if isinstance(u, Polygon):
        u = MultiPolygon([u])
    if u.is_empty or not isinstance(u, MultiPolygon):
        return None
    return u

def norm_name(name):
    s = name.lower().strip()
    s = re.sub(r'\b0+(\d)', r'\1', s)
    s = re.sub(r'^(\d+)\.\s+', r'\1 ', s)
    s = re.sub(r'^wien[- ]', '', s)
    s = re.sub(r'[- ]wien$', '', s)
    s = s.replace('st.', 'sankt ').replace('st ', 'sankt ')
    s = s.replace('dr.', 'doktor ')
    s = s.replace('u.', 'und')
    s = re.sub(r'\s+', ' ', s).strip()
    return s

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 1: Get relation IDs
# ═══════════════════════════════════════════════════════════════════════════════
print("=" * 60)
print("STEP 1: Fetching Vienna district relation IDs")
print("=" * 60)

q = """[out:json][timeout:30];
rel["admin_level"="9"]["boundary"="administrative"](48.1,16.1,48.35,16.6);
out body;"""
j = op_q(q, timeout=40)
rids = sorted([el["id"] for el in j.get("elements", []) if el["type"] == "relation"])
print(f"  {len(rids)} relations: {rids}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 2: Fetch current districts using rel(id);out geom;
# ═══════════════════════════════════════════════════════════════════════════════
print()
print("=" * 60)
print("STEP 2: Fetching current districts")
print("=" * 60)

f_cur = f"{CACHE}/cur.pkl"
if os.path.exists(f_cur):
    with open(f_cur, "rb") as f:
        cur_feats = pickle.load(f)
    print(f"  Loaded {len(cur_feats)} from cache")
else:
    cur_feats = []
    for i, rid in enumerate(rids):
        q = f"[out:json][timeout:30];rel({rid});out geom;"
        try:
            j = op_q(q, timeout=40)
            for el in j.get("elements", []):
                if el["type"] == "relation":
                    tags = el.get("tags", {})
                    geom = members_to_mpoly(el.get("members", []))
                    if geom is not None:
                        n = tags.get("name") or tags.get("name:de") or "unknown"
                        cur_feats.append({"geometry": geom, "name": n})
                        # Quick area check
                        ga = gpd.GeoDataFrame({"geometry": [geom]}, crs="EPSG:4326")
                        gm = ga.to_crs("EPSG:31256")
                        ak = gm.geometry.area[0] / 1e6
                        print(f"  [{i+1}/{len(rids)}] ✓ {n}: {ak:.2f} km²")
                    break
        except Exception as e:
            print(f"  [{i+1}/{len(rids)}] ✗ {rid}: {e!r}")
        time.sleep(0.3)
    with open(f_cur, "wb") as f:
        pickle.dump(cur_feats, f)
    print(f"  Saved {len(cur_feats)} current features")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 3: Fetch 2014 districts
# ═══════════════════════════════════════════════════════════════════════════════
print()
print("=" * 60)
print("STEP 3: Fetching 2014-01-01 districts")
print("=" * 60)

f_old = f"{CACHE}/old.pkl"
if os.path.exists(f_old):
    with open(f_old, "rb") as f:
        old_feats = pickle.load(f)
    print(f"  Loaded {len(old_feats)} from cache")
else:
    old_feats = []
    for i, rid in enumerate(rids):
        q = f"""[out:json][timeout:30][date:"2014-01-01T00:00:00Z"];rel({rid});out geom;"""
        time.sleep(0.5)
        try:
            j = op_q(q, timeout=40)
            for el in j.get("elements", []):
                if el["type"] == "relation":
                    tags = el.get("tags", {})
                    geom = members_to_mpoly(el.get("members", []))
                    if geom is not None:
                        n = tags.get("name") or tags.get("name:de") or "unknown"
                        old_feats.append({"geometry": geom, "name": n})
                        ga = gpd.GeoDataFrame({"geometry": [geom]}, crs="EPSG:4326")
                        gm = ga.to_crs("EPSG:31256")
                        ak = gm.geometry.area[0] / 1e6
                        print(f"  [{i+1}/{len(rids)}] ✓ {n}: {ak:.2f} km²")
                    break
        except Exception as e:
            print(f"  [{i+1}/{len(rids)}] ✗ {rid}: {e!r}")
    with open(f_old, "wb") as f:
        pickle.dump(old_feats, f)
    print(f"  Saved {len(old_feats)} old features")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 4: Match & diff
# ═══════════════════════════════════════════════════════════════════════════════
print()
print("=" * 60)
print("STEP 4: Matching and computing differences")
print("=" * 60)

for f in cur_feats:
    f["name_norm"] = norm_name(f["name"])
for f in old_feats:
    f["name_norm"] = norm_name(f["name"])

obn = {f["name_norm"]: f["geometry"] for f in old_feats}
matches = []
unmatched_c = []
unmatched_o = set(obn.keys())

for f in cur_feats:
    nm = f["name_norm"]
    if nm in obn:
        matches.append((f["geometry"], obn[nm], nm))
        unmatched_o.discard(nm)
    else:
        unmatched_c.append(f)

# Fuzzy
for f in list(unmatched_c):
    nm = f["name_norm"]
    m = re.match(r"^(\d+)", nm)
    if m:
        n = m.group(1)
        for on in list(unmatched_o):
            m2 = re.match(r"^(\d+)", on)
            if m2 and m2.group(1) == n:
                matches.append((f["geometry"], obn[on], nm))
                unmatched_c = [x for x in unmatched_c if x["name_norm"] != nm]
                unmatched_o.discard(on)
                break

print(f"  Matched: {len(matches)}")
if unmatched_c: print(f"  Unmatched cur: {[x['name_norm'] for x in unmatched_c]}")
if unmatched_o: print(f"  Unmatched old: {list(unmatched_o)}")

TOL = 0.00001
fragments = []
for cg, og, nm in matches:
    if isinstance(cg, Polygon): cg = MultiPolygon([cg])
    if isinstance(og, Polygon): og = MultiPolygon([og])
    ob = og.buffer(TOL, join_style=2)
    cb = cg.buffer(TOL, join_style=2)
    try:
        a = cg.difference(ob)
    except:
        a = MultiPolygon()
    try:
        r = og.difference(cb)
    except:
        r = MultiPolygon()
    try:
        u = cg.intersection(ob)
    except:
        u = MultiPolygon()
    for g, ct in [(a, "added_since_2014"), (r, "removed_since_2014"), (u, "unchanged")]:
        if g.is_empty: continue
        if isinstance(g, Polygon): g = MultiPolygon([g])
        if not isinstance(g, MultiPolygon): continue
        fragments.append({"geometry": g, "change_type": ct, "district_name": nm, "touches_changed": False})

for f in unmatched_c:
    g = f["geometry"]
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if isinstance(g, MultiPolygon):
        fragments.append({"geometry": g, "change_type": "added_since_2014", "district_name": f["name_norm"], "touches_changed": True})

for on in unmatched_o:
    g = obn[on]
    if isinstance(g, Polygon): g = MultiPolygon([g])
    if isinstance(g, MultiPolygon):
        fragments.append({"geometry": g, "change_type": "removed_since_2014", "district_name": on, "touches_changed": True})

print(f"  Fragments: {len(fragments)}")
ctc = {}
for f in fragments:
    ctc[f["change_type"]] = ctc.get(f["change_type"], 0) + 1
print(f"  Types: {ctc}")

# Area details
for f in fragments:
    if f["change_type"] == "unchanged":
        continue
    ga = gpd.GeoDataFrame({"geometry": [f["geometry"]]}, crs="EPSG:4326").to_crs("EPSG:31256")
    ak = ga.geometry.area[0] / 1e6
    print(f"  {f['district_name']} {f['change_type']}: {ak:.6f} km²")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 5: touches_changed
# ═══════════════════════════════════════════════════════════════════════════════
print()
print("=" * 60)
print("STEP 5: touches_changed")
print("=" * 60)

cp = [f["geometry"] for f in fragments if f["change_type"] in ("added_since_2014", "removed_since_2014")]
cu = unary_union(cp) if cp else MultiPolygon()

dg = {}
for cg, og, nm in matches:
    dg.setdefault(nm, []).extend([cg, og])
for f in unmatched_c:
    dg.setdefault(f["name_norm"], []).append(f["geometry"])
for on in unmatched_o:
    dg.setdefault(on, []).append(obn[on])

for f in fragments:
    if f["change_type"] in ("added_since_2014", "removed_since_2014"):
        f["touches_changed"] = True
        continue
    dn = f["district_name"]
    if dn in dg and not cu.is_empty:
        du = unary_union([g for g in dg[dn] if not g.is_empty])
        if not du.is_empty:
            try:
                f["touches_changed"] = du.touches(cu)
            except:
                f["touches_changed"] = False
    else:
        f["touches_changed"] = False

tc = sum(1 for f in fragments if f["touches_changed"])
print(f"  touches_changed=True: {tc}/{len(fragments)}")

# ═══════════════════════════════════════════════════════════════════════════════
# STEP 6: Output
# ═══════════════════════════════════════════════════════════════════════════════
print()
print("=" * 60)
print("STEP 6: Writing GeoJSON")
print("=" * 60)

if not fragments:
    out = gpd.GeoDataFrame({"change_type": pd.Series(dtype="str"),
                            "district_name": pd.Series(dtype="str"),
                            "touches_changed": pd.Series(dtype="bool"),
                            "geometry": pd.Series(dtype="object")}, crs="EPSG:4326")
else:
    out = gpd.GeoDataFrame(fragments, crs="EPSG:4326").set_geometry("geometry")
    out["geometry"] = out["geometry"].apply(lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g)

out.to_file("/work/vienna_boundary_changes.geojson", driver="GeoJSON", engine="pyogrio")
print(f"  Written {len(out)} features")
if len(out):
    print(f"  change_type: {out['change_type'].value_counts().to_dict()}")
    print(f"  districts: {sorted(out['district_name'].unique())}")
    print(f"  touches_changed True: {out['touches_changed'].sum()}")

print("\nDone!")