dc-l3-vienna-overpass-historical

Rerun

import requests, json

# Test query for current districts
query = """
[out:json];
area["name"="Wien"]->.a;
rel(area.a)[admin_level=9][boundary=administrative];
out body geom;
"""
r = requests.post("https://overpass-api.de/api/interpreter", data={"data": query})
print(f"Status: {r.status_code}")
data = r.json()
print(f"Elements: {len(data.get('elements', []))}")
for el in data.get('elements', [])[:3]:
    print(json.dumps({k:v for k,v in el.items() if k != 'members'}, indent=2)[:500])
    if 'members' in el:
        print(f"  Members count: {len(el['members'])}")
        for m in el['members'][:3]:
            print(f"    role={m['role']} type={m['type']} ref={m['ref']} geom={'geom' in m}")

import requests, json

# Try with different content type headers
query = """
[out:json];
area["name"="Wien"]->.a;
rel(area.a)[admin_level=9][boundary=administrative];
out body geom;
"""
headers = {"Content-Type": "application/x-www-form-urlencoded"}
r = requests.post("https://overpass-api.de/api/interpreter", data={"data": query}, headers=headers)
print(f"Status: {r.status_code}")
print(f"Response preview: {r.text[:500]}")

import requests, json

query = """
[out:json];
area["name"="Wien"]->.a;
rel(area.a)[admin_level=9][boundary=administrative];
out body geom;
"""
headers = {"Content-Type": "application/x-www-form-urlencoded", "Accept": "application/json"}
r = requests.post("https://overpass-api.de/api/interpreter", data={"data": query}, headers=headers)
print(f"Status: {r.status_code}")
print(f"Response preview: {r.text[:500]}")
if r.status_code == 200:
    data = r.json()
    print(f"Elements: {len(data.get('elements', []))}")

import requests, json

query = """
[out:json];
area["name"="Wien"]->.a;
rel(area.a)[admin_level=9][boundary=administrative];
out body geom;
"""
# Try alternative overpass server
urls = [
    "https://overpass-api.de/api/interpreter",
    "https://lz4.overpass-api.de/api/interpreter",
    "https://overpass.kumi.systems/api/interpreter",
]
for url in urls:
    try:
        r = requests.post(url, data={"data": query}, timeout=30)
        print(f"URL: {url}")
        print(f"Status: {r.status_code}")
        if r.status_code == 200:
            data = r.json()
            print(f"Elements: {len(data.get('elements', []))}")
            break
        else:
            print(f"Response: {r.text[:200]}")
    except Exception as e:
        print(f"Error with {url}: {e}")

import urllib.request, urllib.parse, json

query = '[out:json];area["name"="Wien"]->.a;rel(area.a)[admin_level=9][boundary=administrative];out body geom;'
data = urllib.parse.urlencode({'data': query}).encode()
req = urllib.request.Request('https://overpass-api.de/api/interpreter', data=data, headers={'Accept': 'application/json'})
resp = urllib.request.urlopen(req, timeout=120)
print(f"Status: {resp.status}")
text = resp.read().decode()
print(f"Length: {len(text)}")
print(text[:500])
if resp.status == 200:
    obj = json.loads(text)
    print(f"Elements: {len(obj.get('elements',[]))}")

import urllib.request, urllib.parse, json

query = '[out:json];area["name"="Wien"]->.a;rel(area.a)[admin_level=9][boundary=administrative];out body geom;'
data = urllib.parse.urlencode({'data': query}).encode()
req = urllib.request.Request('https://overpass-api.de/api/interpreter', data=data)
resp = urllib.request.urlopen(req, timeout=120)
print(f"Status: {resp.status}")
text = resp.read().decode()
print(f"Length: {len(text)}")
print(text[:500])
if resp.status == 200:
    obj = json.loads(text)
    print(f"Elements: {len(obj.get('elements',[]))}")

import requests, json

query = '[out:json];area["name"="Wien"]->.a;rel(area.a)[admin_level=9][boundary=administrative];out body geom;'
headers = {
    "User-Agent": "Mozilla/5.0 (X11; Linux x86_64) Python/3.14",
    "Accept": "application/json, text/plain, */*",
    "Content-Type": "application/x-www-form-urlencoded; charset=UTF-8",
}
r = requests.post("https://overpass-api.de/api/interpreter", data={"data": query}, headers=headers, timeout=120)
print(f"Status: {r.status_code}")
print(f"Headers: {dict(r.headers)}")
print(f"Response: {r.text[:500]}")

import requests, json

# Simplified query - no area filter first
query = """
[out:json];
rel(3860455);
out body geom;
"""
# Try using osm3s to prepend _[out:json] differently
query2 = "[out:json];rel(3860455);out body geom;"

# Also try getting node
query3 = "[out:json];node(1);out;"

for i, q in enumerate([query, query2, query3], 1):
    try:
        r = requests.post("https://overpass-api.de/api/interpreter", data={"data": q.strip()}, timeout=30)
        print(f"Query {i}: Status={r.status_code}")
        if r.status_code == 200:
            print(f"  Elements: {len(r.json().get('elements',[]))}")
            print(f"  Text: {r.text[:300]}")
        else:
            print(f"  Error: {r.text[:200]}")
    except Exception as e:
        print(f"Query {i}: Exception: {e}")

import json, io, urllib.request, urllib.parse

query = '[out:json];area["name"="Wien"]->.a;rel(area.a)[admin_level=9][boundary=administrative];out body geom;'
data = urllib.parse.urlencode({'data': query}).encode()
req = urllib.request.Request('https://overpass-api.de/api/interpreter', data=data)
# Remove all accept headers by not setting one
resp = urllib.request.urlopen(req, timeout=120)
print(f"Status: {resp.status}")
print(f"Headers: {dict(resp.headers)}")
text = resp.read().decode()
print(f"Length: {len(text)}")
print(text[:500])

import pycurl, io, json, urllib.parse

query = '[out:json];area["name"="Wien"]->.a;rel(area.a)[admin_level=9][boundary=administrative];out body geom;'

buf = io.BytesIO()
c = pycurl.Curl()
c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
c.setopt(c.WRITEDATA, buf)
c.setopt(c.TIMEOUT, 120)
c.perform()
status = c.getinfo(c.RESPONSE_CODE)
body = buf.getvalue().decode('utf-8')
print(f"Status: {status}")
print(f"Length: {len(body)}")
print(body[:500])
if status == 200:
    obj = json.loads(body)
    print(f"Elements: {len(obj.get('elements',[]))}")

import pycurl, io, json, urllib.parse, pickle

query = '[out:json];area["name"="Wien"]->.a;rel(area.a)[admin_level=9][boundary=administrative];out body geom;'

buf = io.BytesIO()
c = pycurl.Curl()
c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
c.setopt(c.WRITEDATA, buf)
c.setopt(c.TIMEOUT, 120)
c.perform()
status = c.getinfo(c.RESPONSE_CODE)
body = buf.getvalue().decode('utf-8')
print(f"Status: {status}, Length: {len(body)}")
obj = json.loads(body)
print(f"Elements: {len(obj.get('elements',[]))}")
with open('/work/current_districts.json', 'w') as f:
    json.dump(obj, f)
print("Saved current_districts.json")

# Show element types and names
for el in obj['elements']:
    tags = el.get('tags', {})
    name = tags.get('name', '?')
    print(f"  id={el['id']} type={el['type']} name={name}")

import pycurl, io, json, urllib.parse

# Use the 'date' parameter in the [out:json] or use a query with date
query = '[out:json][date:"2014-01-01T00:00:00Z"];area["name"="Wien"]->.a;rel(area.a)[admin_level=9][boundary=administrative];out body geom;'

buf = io.BytesIO()
c = pycurl.Curl()
c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
c.setopt(c.WRITEDATA, buf)
c.setopt(c.TIMEOUT, 120)
c.perform()
status = c.getinfo(c.RESPONSE_CODE)
body = buf.getvalue().decode('utf-8')
print(f"Status: {status}, Length: {len(body)}")
if status == 200:
    obj = json.loads(body)
    print(f"Elements: {len(obj.get('elements',[]))}")
    with open('/work/districts_2014.json', 'w') as f:
        json.dump(obj, f)
    print("Saved districts_2014.json")
    for el in obj['elements']:
        tags = el.get('tags', {})
        name = tags.get('name', '?')
        print(f"  id={el['id']} type={el['type']} name={name}")
else:
    print(body[:500])

import pycurl, io, json, urllib.parse

# Method: use the 'out' modifier with 'geom' and specify date in the request
# Actually Overpass supports [date:"YYYY-MM-DDTHH:MM:SSZ"] in the query header
# But maybe the server doesn't have attic data for that date.
# Let's try without the area filter first, get all relations with admin_level=9 in Vienna 

# Try a different syntax - put date right at the beginning
query = '[date:"2014-01-01T00:00:00Z"][out:json];area["name"="Wien"]->.a;rel(area.a)[admin_level=9][boundary=administrative];out body geom;'

buf = io.BytesIO()
c = pycurl.Curl()
c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
c.setopt(c.WRITEDATA, buf)
c.setopt(c.TIMEOUT, 120)
c.setopt(c.NOSIGNAL, 1)
c.perform()
status = c.getinfo(c.RESPONSE_CODE)
body = buf.getvalue().decode('utf-8')
print(f"Status: {status}, Length: {len(body)}")
if status == 200:
    obj = json.loads(body)
    print(f"Elements: {len(obj.get('elements',[]))}")
    with open('/work/districts_2014.json', 'w') as f:
        json.dump(obj, f)
    print("Saved districts_2014.json")
    for el in obj['elements']:
        tags = el.get('tags', {})
        name = tags.get('name', '?')
        print(f"  id={el['id']} type={el['type']} name={name}")
else:
    print(body[:500])

import pycurl, io, json, urllib.parse

# Let's get all Vienna districts by their relation IDs directly
# First get current to find the IDs
# Actually we know they're 1990590-1990600, 1991416, 1991433-1991443
# Let's try with a simpler query first to see if attic data works

# Try getting a single relation at that date
query = '[out:json][date:"2014-01-01T00:00:00Z"];rel(1990590);out body geom;'

buf = io.BytesIO()
c = pycurl.Curl()
c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
c.setopt(c.WRITEDATA, buf)
c.setopt(c.TIMEOUT, 120)
c.setopt(c.NOSIGNAL, 1)
c.perform()
status = c.getinfo(c.RESPONSE_CODE)
body = buf.getvalue().decode('utf-8')
print(f"Status: {status}, Length: {len(body)}")
if status == 200:
    obj = json.loads(body)
    print(f"Elements: {len(obj.get('elements',[]))}")
    for el in obj['elements']:
        tags = el.get('tags', {})
        name = tags.get('name', '?')
        print(f"  id={el['id']} name={name}")
    with open('/work/test_2014.json', 'w') as f:
        json.dump(obj, f)
else:
    print(body[:500])

import pycurl, io, json, urllib.parse

# Get all 23 Vienna district relations from 2014-01-01
ids = [1990590, 1990591, 1990592, 1990593, 1990594, 1990595, 1990596, 
       1990597, 1990598, 1990599, 1990600, 1991416, 1991433, 1991434, 
       1991435, 1991436, 1991437, 1991438, 1991439, 1991440, 1991441, 
       1991442, 1991443]

# Build query for all relations
rel_ids = ','.join(str(i) for i in ids)
query = f'[out:json][date:"2014-01-01T00:00:00Z"];rel({rel_ids});out body geom;'

buf = io.BytesIO()
c = pycurl.Curl()
c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
c.setopt(c.WRITEDATA, buf)
c.setopt(c.TIMEOUT, 300)
c.setopt(c.NOSIGNAL, 1)
c.perform()
status = c.getinfo(c.RESPONSE_CODE)
body = buf.getvalue().decode('utf-8')
print(f"Status: {status}, Length: {len(body)}")
if status == 200:
    obj = json.loads(body)
    print(f"Elements: {len(obj.get('elements',[]))}")
    with open('/work/districts_2014.json', 'w') as f:
        json.dump(obj, f)
    print("Saved districts_2014.json")
    for el in obj['elements']:
        tags = el.get('tags', {})
        name = tags.get('name', '?')
        print(f"  id={el['id']} name={name}")
else:
    print(body[:1000])

import pycurl, io, json, urllib.parse, time

ids = [1990590, 1990591, 1990592, 1990593, 1990594, 1990595, 1990596, 
       1990597, 1990598, 1990599, 1990600, 1991416, 1991433, 1991434, 
       1991435, 1991436, 1991437, 1991438, 1991439, 1991440, 1991441, 
       1991442, 1991443]

all_elements = []

for rid in ids:
    query = f'[out:json][date:"2014-01-01T00:00:00Z"];rel({rid});out body geom;'
    buf = io.BytesIO()
    c = pycurl.Curl()
    c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
    c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
    c.setopt(c.WRITEDATA, buf)
    c.setopt(c.TIMEOUT, 120)
    c.setopt(c.NOSIGNAL, 1)
    c.perform()
    status = c.getinfo(c.RESPONSE_CODE)
    body = buf.getvalue().decode('utf-8')
    if status == 200:
        obj = json.loads(body)
        for el in obj['elements']:
            all_elements.append(el)
            tags = el.get('tags', {})
            name = tags.get('name', '?')
            print(f"  id={el['id']} name={name}")
    else:
        print(f"  id={rid} status={status} error={body[:200]}")
    time.sleep(2)  # be nice to the server

print(f"\nTotal elements: {len(all_elements)}")
with open('/work/districts_2014.json', 'w') as f:
    json.dump({"version": 0.6, "generator": "Overpass", "elements": all_elements}, f)
print("Saved districts_2014.json")

import pycurl, io, json, urllib.parse, time

ids = [1991442, 1991443]
all_elements = []

for rid in ids:
    query = f'[out:json][date:"2014-01-01T00:00:00Z"];rel({rid});out body geom;'
    buf = io.BytesIO()
    c = pycurl.Curl()
    c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
    c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
    c.setopt(c.WRITEDATA, buf)
    c.setopt(c.TIMEOUT, 120)
    c.setopt(c.NOSIGNAL, 1)
    c.setopt(c.FRESH_CONNECT, 1)
    c.perform()
    status = c.getinfo(c.RESPONSE_CODE)
    body = buf.getvalue().decode('utf-8')
    if status == 200:
        obj = json.loads(body)
        for el in obj['elements']:
            all_elements.append(el)
            tags = el.get('tags', {})
            name = tags.get('name', '?')
            print(f"  id={el['id']} name={name}")
    else:
        print(f"  id={rid} status={status}")
        print(f"  {body[:200]}")
    time.sleep(5)

with open('/work/districts_2014.json', 'r') as f:
    data = json.load(f)
data['elements'].extend(all_elements)
with open('/work/districts_2014.json', 'w') as f:
    json.dump(data, f)
print(f"Total: {len(data['elements'])}")

import json

with open('/work/current_districts.json') as f:
    cur = json.load(f)
with open('/work/districts_2014.json') as f:
    old = json.load(f)

# Examine a relation structure
for el in cur['elements']:
    if el['id'] == 1990592:  # Innere Stadt
        print("=== Current Innere Stadt ===")
        print(f"Tags: {json.dumps(el.get('tags',{}), indent=2)}")
        mems = el.get('members', [])
        print(f"Members: {len(mems)}")
        for m in mems[:5]:
            print(f"  role={m['role']} type={m['type']} ref={m['ref']} geom={'geom' in m}")
        # Check if any member has full geometry
        for m in mems:
            if 'geom' in m:
                print(f"  Member with geom: role={m['role']} type={m['type']}")
                print(f"    Geometry points: {len(m['geom'])}")
                print(f"    First 3: {m['geom'][:3]}")
                break
        break

# Same for 2014
for el in old['elements']:
    if el['id'] == 1990592:
        print("\n=== 2014 Innere Stadt ===")
        print(f"Tags: {json.dumps(el.get('tags',{}), indent=2)}")
        mems = el.get('members', [])
        print(f"Members: {len(mems)}")
        for m in mems[:5]:
            print(f"  role={m['role']} type={m['type']} ref={m['ref']} geom={'geom' in m}")
        for m in mems:
            if 'geom' in m:
                print(f"  Member with geom: role={m['role']} type={m['type']}")
                print(f"    Geometry points: {len(m['geom'])}")
                print(f"    First 3: {m['geom'][:3]}")
                break
        break

import json

with open('/work/current_districts.json') as f:
    cur = json.load(f)

# Check if any element has geometry directly
for el in cur['elements']:
    if el['type'] == 'relation' and el['id'] == 1990592:
        # Check if 'geometry' key exists
        print(f"Keys: {list(el.keys())}")
        # The members
        for m in el.get('members', []):
            if m['role'] == 'outer':
                print(f"Way ref={m['ref']}, has_geom={'geom' in m}")
        break
    
# Check the whole response for any ways
ways = [el for el in cur['elements'] if el['type'] == 'way']
print(f"\nTotal ways in response: {len(ways)}")
if ways:
    print(f"First way keys: {list(ways[0].keys())}")
    print(f"First way nodes: {ways[0].get('nodes', [])[:5]}")
    if 'geometry' in ways[0]:
        print(f"Geometry: {ways[0]['geometry'][:3]}")
    # Check if there's a 'geom' key (Overpass often calls it 'geometry')
    if 'geom' in ways[0]:
        print(f"geom: {ways[0]['geom'][:3]}")

import pycurl, io, json, urllib.parse

# Try different query styles to get geometry for relation members

# Style 1: out geom
query1 = '[out:json];rel(1990592);out geom;'

# Style 2: recurse then out
query2 = '[out:json];rel(1990592);(._;>;);out geom;'

for label, query in [("out geom", query1), ("recurse", query2)]:
    buf = io.BytesIO()
    c = pycurl.Curl()
    c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
    c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
    c.setopt(c.WRITEDATA, buf)
    c.setopt(c.TIMEOUT, 30)
    c.setopt(c.NOSIGNAL, 1)
    c.perform()
    status = c.getinfo(c.RESPONSE_CODE)
    body = buf.getvalue().decode('utf-8')
    print(f"\n=== {label} (status={status}) ===")
    if status == 200:
        obj = json.loads(body)
        print(f"Elements: {len(obj.get('elements',[]))}")
        for el in obj['elements']:
            if el['type'] == 'relation':
                print(f"  Relation id={el['id']}")
                for m in el.get('members', []):
                    print(f"    Member: role={m['role']} type={m['type']} ref={m['ref']} geom={'geom' in m}")
            elif el['type'] == 'way':
                print(f"  Way id={el['id']} nodes={len(el.get('nodes',[]))} geom={'geom' in el or 'geometry' in el}")
                if 'geometry' in el:
                    print(f"    Geometry points: {len(el['geometry'])} first: {el['geometry'][:2]}")
                elif 'geom' in el:
                    print(f"    geom points: {len(el['geom'])}")
            elif el['type'] == 'node':
                if len(obj['elements']) < 20:
                    print(f"  Node id={el['id']} lat={el.get('lat')} lon={el.get('lon')}")
    else:
        print(body[:300])

import pycurl, io, json, urllib.parse

ids = [1990590, 1990591, 1990592, 1990593, 1990594, 1990595, 1990596, 
       1990597, 1990598, 1990599, 1990600, 1991416, 1991433, 1991434, 
       1991435, 1991436, 1991437, 1991438, 1991439, 1991440, 1991441, 
       1991442, 1991443]

rel_ids = ','.join(str(i) for i in ids)
# Query: get relations, recurse to ways with geometry
query = f'[out:json];rel({rel_ids});(._;>;);out geom;'

buf = io.BytesIO()
c = pycurl.Curl()
c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
c.setopt(c.WRITEDATA, buf)
c.setopt(c.TIMEOUT, 300)
c.setopt(c.NOSIGNAL, 1)
c.perform()
status = c.getinfo(c.RESPONSE_CODE)
body = buf.getvalue().decode('utf-8')
print(f"Status: {status}, Length: {len(body)}")
if status == 200:
    obj = json.loads(body)
    print(f"Elements: {len(obj.get('elements',[]))}")
    with open('/work/current_full.json', 'w') as f:
        json.dump(obj, f)
    print("Saved current_full.json")
    
    # Count types
    types = {}
    for el in obj['elements']:
        t = el['type']
        types[t] = types.get(t, 0) + 1
    print(f"Types: {types}")
else:
    print(body[:500])

import pycurl, io, json, urllib.parse, time

ids = [1990590, 1990591, 1990592, 1990593, 1990594, 1990595, 1990596, 
       1990597, 1990598, 1990599, 1990600, 1991416, 1991433, 1991434, 
       1991435, 1991436, 1991437, 1991438, 1991439, 1991440, 1991441, 
       1991442, 1991443]

def fetch_batch(rids, label):
    rel_ids = ','.join(str(i) for i in rids)
    query = f'[out:json];rel({rel_ids});(._;>;);out geom;'
    buf = io.BytesIO()
    c = pycurl.Curl()
    c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
    c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
    c.setopt(c.WRITEDATA, buf)
    c.setopt(c.TIMEOUT, 120)
    c.setopt(c.NOSIGNAL, 1)
    c.perform()
    status = c.getinfo(c.RESPONSE_CODE)
    body = buf.getvalue().decode('utf-8')
    print(f"  {label}: status={status}, length={len(body)}")
    if status == 200:
        return json.loads(body)
    else:
        print(f"  Error: {body[:300]}")
        return None

all_data = {"version": 0.6, "generator": "Overpass", "elements": []}

# Try batches of 5
batch_size = 5
for i in range(0, len(ids), batch_size):
    batch = ids[i:i+batch_size]
    label = f"ids {batch[0]}-{batch[-1]}"
    print(f"Fetching {label}...")
    data = fetch_batch(batch, label)
    if data:
        all_data['elements'].extend(data['elements'])
    time.sleep(3)

print(f"\nTotal elements collected: {len(all_data['elements'])}")
with open('/work/current_full.json', 'w') as f:
    json.dump(all_data, f)
print("Saved current_full.json")

import pycurl, io, json, urllib.parse, time

ids = [1990590, 1990591, 1990592, 1990593, 1990594, 1990595, 1990596, 
       1990597, 1990598, 1990599, 1990600, 1991416, 1991433, 1991434, 
       1991435, 1991436, 1991437, 1991438, 1991439, 1991440, 1991441, 
       1991442, 1991443]

all_elements = []

for i, rid in enumerate(ids):
    query = f'[out:json];rel({rid});(._;>;);out geom;'
    buf = io.BytesIO()
    c = pycurl.Curl()
    c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
    c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
    c.setopt(c.WRITEDATA, buf)
    c.setopt(c.TIMEOUT, 60)
    c.setopt(c.NOSIGNAL, 1)
    c.perform()
    status = c.getinfo(c.RESPONSE_CODE)
    body = buf.getvalue().decode('utf-8')
    
    tags = {}
    if status == 200:
        data = json.loads(body)
        for el in data['elements']:
            all_elements.append(el)
            if el['type'] == 'relation':
                tags = el.get('tags', {})
        name = tags.get('name', '?')
        print(f"  [{i+1}/{len(ids)}] id={rid} name={name} status=200 elements={len(data['elements'])}")
    else:
        print(f"  [{i+1}/{len(ids)}] id={rid} status={status}")
        if status == 429:
            print("  Rate limited, waiting 10s...")
            time.sleep(10)
            # Retry
            c.perform()
            status2 = c.getinfo(c.RESPONSE_CODE)
            body2 = buf.getvalue().decode('utf-8')
            if status2 == 200:
                data = json.loads(body2)
                for el in data['elements']:
                    all_elements.append(el)
                    if el['type'] == 'relation':
                        tags = el.get('tags', {})
                name = tags.get('name', '?')
                print(f"  Retry: id={rid} name={name} status=200 elements={len(data['elements'])}")
    
    # Be nice to server
    time.sleep(5)

print(f"\nTotal elements: {len(all_elements)}")
with open('/work/current_full.json', 'w') as f:
    json.dump({"version": 0.6, "generator": "Overpass", "elements": all_elements}, f)
print("Saved current_full.json")

import pycurl, io, json, urllib.parse

# Check full error for batch of 5
ids = [1990590, 1990591, 1990592, 1990593, 1990594]
rel_ids = ','.join(str(i) for i in ids)
query = f'[out:json];rel({rel_ids});(._;>;);out geom;'
print(f"Query: {query}")

buf = io.BytesIO()
c = pycurl.Curl()
c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
c.setopt(c.WRITEDATA, buf)
c.setopt(c.TIMEOUT, 60)
c.setopt(c.NOSIGNAL, 1)
c.perform()
status = c.getinfo(c.RESPONSE_CODE)
body = buf.getvalue().decode('utf-8')
print(f"Status: {status}")
print(f"Full body:\n{body}")

import pycurl, io, json, urllib.parse, time

query = '[out:json];rel(1990590);(._;>;);out geom;'
print(f"Query: {query}")

buf = io.BytesIO()
c = pycurl.Curl()
c.setopt(c.URL, 'https://overpass-api.de/api/interpreter')
c.setopt(c.POSTFIELDS, 'data=' + urllib.parse.quote(query))
c.setopt(c.WRITEDATA, buf)
c.setopt(c.TIMEOUT, 60)
c.setopt(c.NOSIGNAL, 1)
c.setopt(c.FRESH_CONNECT, 1)
c.setopt(c.FORBID_REUSE, 1)
c.perform()
status = c.getinfo(c.RESPONSE_CODE)
body = buf.getvalue().decode('utf-8')
print(f"Status: {status}, Length: {len(body)}")
if status == 200:
    data = json.loads(body)
    print(f"Elements: {len(data['elements'])}")
    # Count types
    types = {}
    for el in data['elements']:
        t = el['type']
        types[t] = types.get(t, 0) + 1
    print(f"Types: {types}")
else:
    print(body[:500])

import json

# Check what geometry data we already have
for fname, label in [('/work/current_districts.json', 'current'), ('work/districts_2014.json', '2014')]:
    try:
        with open(fname) as f:
            data = json.load(f)
        print(f"\n=== {label} ===")
        print(f"Elements: {len(data['elements'])}")
        for el in data['elements']:
            if el['type'] == 'relation':
                tags = el.get('tags', {})
                name = tags.get('name', '?')
                mems = el.get('members', [])
                outer_count = sum(1 for m in mems if m['role'] == 'outer')
                has_geom_members = any('geom' in m for m in mems)
                print(f"  id={el['id']} name={name} outer_ways={outer_count} has_member_geom={has_geom_members}")
                if has_geom_members:
                    for m in mems:
                        if 'geom' in m:
                            print(f"    First member with geom: role={m['role']} pts={len(m['geom'])}")
                            print(f"    First pt: {m['geom'][0]}")
                            break
    except FileNotFoundError:
        print(f"\nFile not found: {fname}")

import json

# Check if members have 'geom' key in the data we fetched
with open('/work/current_districts.json') as f:
    cur = json.load(f)
with open('/work/districts_2014.json') as f:
    old = json.load(f)

# Check current
print("=== Current ===")
for el in cur['elements']:
    if el['type'] == 'relation':
        name = el.get('tags', {}).get('name', '?')
        mems = el.get('members', [])
        for m in mems:
            if 'geom' in m:
                print(f"  id={el['id']} name={name}: member type={m['type']} ref={m['ref']} geom_pts={len(m['geom'])}")
                break
        else:
            print(f"  id={el['id']} name={name}: NO geometry on members")
            # Check if 'geometry' key exists
            for m in mems:
                if 'geometry' in m:
                    print(f"    -> Found 'geometry' on member!")
                    break
            else:
                print(f"    -> Member keys sample: {list(mems[0].keys()) if mems else 'no members'}")
                # Check if relation itself has geometry
                if 'geometry' in el:
                    print(f"    -> Relation has geometry key!")
                if 'geom' in el:
                    print(f"    -> Relation has geom key!")

print("\n=== 2014 ===")
for el in old['elements']:
    if el['type'] == 'relation':
        name = el.get('tags', {}).get('name', '?')
        mems = el.get('members', [])
        for m in mems:
            if 'geom' in m:
                print(f"  id={el['id']} name={name}: member type={m['type']} ref={m['ref']} geom_pts={len(m['geom'])}")
                break
        else:
            print(f"  id={el['id']} name={name}: NO geometry on members")

import json

# Check geometry details
with open('/work/current_districts.json') as f:
    cur = json.load(f)
with open('/work/districts_2014.json') as f:
    old = json.load(f)

# Check first member of first relation in current  
el = cur['elements'][0]
print(f"Current first element: id={el['id']} type={el['type']}")
mems = el.get('members', [])
if mems:
    m = mems[0]
    print(f"  First member keys: {list(m.keys())}")
    print(f"  type={m['type']} ref={m['ref']} role={m['role']}")
    if 'geometry' in m:
        print(f"  geometry len: {len(m['geometry'])}")
        print(f"  first pt: {m['geometry'][0]}")

# Check 2014
el = old['elements'][0]
print(f"\n2014 first element: id={el['id']} type={el['type']}")
mems = el.get('members', [])
if mems:
    m = mems[0]
    print(f"  First member keys: {list(m.keys())}")
    print(f"  type={m['type']} ref={m['ref']} role={m['role']}")

# Let's check all 2014 members for a single district
for el in old['elements']:
    if el['id'] == 1990590:
        print(f"\n2014 Alsergrund members:")
        for m in el.get('members', []):
            print(f"  type={m['type']} ref={m['ref']} role={m['role']} keys={list(m.keys())}")
        break

# And for current
for el in cur['elements']:
    if el['id'] == 1990590:
        print(f"\nCurrent Alsergrund members:")
        for m in el.get('members', []):
            print(f"  type={m['type']} ref={m['ref']} role={m['role']} keys={list(m.keys())}")
        break

import json

with open('/work/current_districts.json') as f:
    cur = json.load(f)
with open('/work/districts_2014.json') as f:
    old = json.load(f)

# Check 2014 geometry data
for el in old['elements']:
    if el['id'] == 1990590:
        for m in el.get('members', []):
            if 'geometry' in m:
                print(f"2014: type={m['type']} ref={m['ref']} role={m['role']} geometry_pts={len(m['geometry'])}")
                print(f"  first pt: {m['geometry'][0]}")
                break
        break

# Check if 2014 data has geometry for all
for el in old['elements']:
    if el['type'] == 'relation':
        name = el.get('tags', {}).get('name', '?')
        mems = el.get('members', [])
        has_geom = any('geometry' in m for m in mems if m['role'] == 'outer')
        if not has_geom:
            print(f"2014: {name} has NO geometry on outer ways!")
        else:
            print(f"2014: {name} has geometry OK")

import json

with open('/work/current_districts.json') as f:
    cur = json.load(f)
with open('/work/districts_2014.json') as f:
    old = json.load(f)

print("=== Current names ===")
for el in sorted(cur['elements'], key=lambda x: int(x.get('tags',{}).get('ref', '99'))):
    tags = el.get('tags', {})
    print(f"  ref={tags.get('ref','?')} name='{tags.get('name','?')}'")

print("\n=== 2014 names ===")
for el in sorted(old['elements'], key=lambda x: int(x.get('tags',{}).get('ref', '99'))):
    tags = el.get('tags', {})
    print(f"  ref={tags.get('ref','?')} name='{tags.get('name','?')}'")

import json
from shapely.geometry import Polygon, MultiPolygon, shape, mapping
from shapely.ops import unary_union
import geopandas as gpd

def build_polygon_from_relation(rel):
    """Build a Polygon from a relation's outer way geometries."""
    outer_rings = []
    for m in rel.get('members', []):
        if m['role'] == 'outer' and 'geometry' in m:
            coords = [(pt['lon'], pt['lat']) for pt in m['geometry']]
            if len(coords) >= 3:
                outer_rings.append(coords)
    
    if not outer_rings:
        return None
    
    # Try to merge into a single polygon (most districts have 1 outer ring)
    # For districts with multiple outer ways, we need to assemble them
    # First, collect all rings as linestrings
    from shapely.geometry import LineString
    from shapely.ops import polygonize
    
    lines = [LineString(ring) for ring in outer_rings]
    
    # Try simple approach: if the first and last points match, it's a closed ring
    # Use polygonize to merge the lines
    merged = unary_union(lines)
    
    if merged.geom_type == 'LineString':
        # Single line - try to close it
        if merged.is_ring:
            return Polygon(merged.coords)
        else:
            print(f"  WARNING: merged line is not a ring for rel {rel['id']}")
            return None
    elif merged.geom_type == 'MultiLineString':
        # Try polygonizing
        polygons = list(polygonize(merged.geoms))
        if len(polygons) == 1:
            return polygons[0]
        elif len(polygons) > 1:
            return MultiPolygon(polygons)
        else:
            print(f"  WARNING: cannot polygonize for rel {rel['id']}")
            return None
    else:
        return None


def build_polygon_simple(rel):
    """
    Build polygon from ways in order. 
    Overpass returns geometry arrays which are sequential segments.
    We need to connect them properly.
    """
    from shapely.geometry import LineString, Polygon
    
    outer_ways = []
    for m in rel.get('members', []):
        if m['role'] == 'outer' and 'geometry' in m:
            coords = [(pt['lon'], pt['lat']) for pt in m['geometry']]
            if len(coords) >= 2:
                outer_ways.append(coords)
    
    if not outer_ways:
        return None
    
    # For Vienna districts, outer ways form closed rings
    # Let's try to connect them by matching endpoints
    # First, let's see if there's just one closed ring
    if len(outer_ways) == 1 and len(outer_ways[0]) >= 3:
        coords = outer_ways[0]
        if coords[0] == coords[-1]:
            return Polygon(coords)
        else:
            # Try to close it
            return Polygon(coords + [coords[0]])
    
    # Multiple ways - sort by connecting endpoints
    # Build a dict mapping endpoint tuples to line segments
    segments = []
    for coords in outer_ways:
        seg = LineString(coords)
        segments.append(seg)
    
    merged = unary_union(segments)
    
    if merged.geom_type == 'LineString':
        if merged.is_ring:
            return Polygon(merged.coords)
        else:
            # Try to close
            coords = list(merged.coords)
            if coords[0] != coords[-1]:
                coords.append(coords[0])
            return Polygon(coords)
    elif merged.geom_type == 'MultiLineString':
        polygs = list(polygonize(merged.geoms))
        if len(polygs) == 1:
            return polygs[0]
        elif len(polygs) > 1:
            # Check for holes - inner rings
            # For Vienna districts, there might be holes (enclaves)
            # Sort by area, largest is outer
            polygs_sorted = sorted(polygs, key=lambda p: p.area, reverse=True)
            outer = polygs_sorted[0]
            holes = [p.exterior.coords for p in polygs_sorted[1:] if p.area > 0.000001]
            if holes:
                return Polygon(outer.exterior.coords, holes)
            return outer
        else:
            print(f"  WARNING: cannot polygonize for rel {rel['id']}")
            return None
    else:
        return None


def build_polygon_simple2(rel):
    """
    Simpler approach: just collect all outer way coordinates,
    close them, and build polygon via unary_union of LineStrings.
    """
    from shapely.geometry import LineString
    from shapely.ops import polygonize
    
    lines = []
    for m in rel.get('members', []):
        if m['role'] == 'outer' and 'geometry' in m:
            coords = [(pt['lon'], pt['lat']) for pt in m['geometry']]
            if len(coords) >= 2:
                lines.append(LineString(coords))
    
    if not lines:
        return None
    
    merged = unary_union(lines)
    
    if merged.geom_type == 'GeometryCollection':
        geoms = [g for g in merged.geoms if g.geom_type in ('LineString', 'MultiLineString')]
        if not geoms:
            return None
        merged = unary_union(geoms)
    
    if merged.geom_type in ('LineString',):
        if merged.is_ring:
            return Polygon(merged)
        else:
            coords = list(merged.coords)
            if len(coords) >= 3:
                if coords[0] != coords[-1]:
                    coords.append(coords[0])
                return Polygon(coords)
            return None
    elif merged.geom_type == 'MultiLineString':
        polys = list(polygonize(merged.geoms))
        if not polys:
            return None
        if len(polys) == 1:
            return polys[0]
        # Sort by area
        polys.sort(key=lambda p: p.area, reverse=True)
        outer_poly = polys[0]
        holes = []
        for p in polys[1:]:
            if p.representative_point().within(outer_poly) and not p.touches(outer_poly):
                holes.append(p.exterior.coords)
        if holes:
            return Polygon(outer_poly.exterior.coords, holes)
        return outer_poly
    else:
        return None


def build_polygon_robust(rel):
    """
    Build a polygon from the relation's outer way geometries.
    This handles Vienna's district boundaries which are made of multiple ways.
    """
    from shapely.geometry import LineString, Polygon, MultiLineString
    from shapely.ops import polygonize, linemerge
    
    # Collect all outer way geometries as LineStrings
    lines = []
    for m in rel.get('members', []):
        if m['role'] == 'outer' and 'geometry' in m:
            coords = [(pt['lon'], pt['lat']) for pt in m['geometry']]
            if len(coords) >= 2:
                lines.append(LineString(coords))
    
    if not lines:
        return None
    
    # Merge connected lines
    merged = linemerge(lines)
    
    if merged.geom_type == 'LineString':
        coords = list(merged.coords)
        if coords[0] != coords[-1]:
            coords.append(coords[0])
        return Polygon(coords)
    elif merged.geom_type == 'MultiLineString':
        # Try to polygonize
        polys = list(polygonize(merged.geoms))
        if not polys:
            return None
        
        if len(polys) == 1:
            return polys[0]
        
        # Multiple polygons - check for holes
        polys_sorted = sorted(polys, key=lambda p: p.area, reverse=True)
        outer = polys_sorted[0]
        holes = []
        for p in polys_sorted[1:]:
            # Check if this polygon is a hole (inner ring)
            if p.area > 1e-10:
                # Check if it's inside the outer polygon
                centroid = p.representative_point()
                if outer.contains(centroid) or outer.touches(p):
                    holes.append(p.exterior.coords)
        if holes:
            return Polygon(outer.exterior.coords, holes)
        return outer
    else:
        return None


# Load data
print("Loading data...")
with open('/work/current_districts.json') as f:
    current_data = json.load(f)
with open('/work/districts_2014.json') as f:
    old_data = json.load(f)

# Build district polygons
print("Building current district polygons...")
current_polys = {}
for el in current_data['elements']:
    if el['type'] == 'relation':
        ref = el.get('tags', {}).get('ref')
        name = el.get('tags', {}).get('name', '')
        poly = build_polygon_robust(el)
        if poly:
            current_polys[ref] = {'name': name, 'polygon': poly, 'relation': el}
        else:
            print(f"  FAILED to build polygon for current: {name} (ref={ref})")

print("Building 2014 district polygons...")
old_polys = {}
for el in old_data['elements']:
    if el['type'] == 'relation':
        ref = el.get('tags', {}).get('ref')
        name = el.get('tags', {}).get('name', '')
        poly = build_polygon_robust(el)
        if poly:
            old_polys[ref] = {'name': name, 'polygon': poly, 'relation': el}
        else:
            print(f"  FAILED to build polygon for 2014: {name} (ref={ref})")

print(f"\nCurrent districts built: {len(current_polys)}")
print(f"2014 districts built: {len(old_polys)}")

# Check which refs match
print(f"\nCurrent refs: {sorted(current_polys.keys())}")
print(f"2014 refs: {sorted(old_polys.keys())}")

# Build complete Vienna boundaries for both time periods
print("\nComputing Vienna boundary differences...")
current_vienna = unary_union([current_polys[r]['polygon'] for r in current_polys])
old_vienna = unary_union([old_polys[r]['polygon'] for r in old_polys])

# Compute symmetric difference
# added_since_2014 = in current but not in 2014
added = current_vienna.difference(old_vienna)
# removed_since_2014 = in 2014 but not in current
removed = old_vienna.difference(current_vienna)
# unchanged = intersection
unchanged = current_vienna.intersection(old_vienna)

print(f"  Added area: {added.area:.6f} (type={added.geom_type})")
print(f"  Removed area: {removed.area:.6f} (type={removed.geom_type})")
print(f"  Unchanged area: {unchanged.area:.6f}")

# Build features
features = []

def geom_to_feature(geom, change_type, district_name, touches_changed):
    """Convert a shapely geometry to a GeoJSON feature."""
    if geom.is_empty:
        return None
    if geom.geom_type == 'GeometryCollection':
        # Flatten
        result = []
        for g in geom.geoms:
            f = geom_to_feature(g, change_type, district_name, touches_changed)
            if f:
                result.append(f)
        return result
    
    # Ensure MultiPolygon
    if geom.geom_type == 'Polygon':
        geom = MultiPolygon([geom])
    elif geom.geom_type != 'MultiPolygon':
        return None
    
    return {
        "type": "Feature",
        "geometry": mapping(geom),
        "properties": {
            "change_type": change_type,
            "district_name": district_name,
            "touches_changed": touches_changed
        }
    }

# For added and removed, the district_name is not clearly one district
# We'll classify by which districts are involved
# For added: which current districts overlap with the added area
# For removed: which old districts overlap with the removed area

# Compute per-district changes
features_list = []

# Unchanged - per district
for ref in sorted(current_polys.keys(), key=int):
    name = current_polys[ref]['name']
    poly = current_polys[ref]['polygon']
    
    district_unchanged = poly.intersection(unchanged)
    district_added = poly.intersection(added)
    district_removed = poly.intersection(removed)
    
    touches = not district_added.is_empty or not district_removed.is_empty
    
    if not district_unchanged.is_empty:
        f = geom_to_feature(district_unchanged, 'unchanged', name, touches)
        if f:
            if isinstance(f, list):
                features_list.extend(f)
            else:
                features_list.append(f)
    
    if not district_added.is_empty:
        f = geom_to_feature(district_added, 'added_since_2014', name, True)
        if f:
            if isinstance(f, list):
                features_list.extend(f)
            else:
                features_list.append(f)
    
    if not district_removed.is_empty:
        f = geom_to_feature(district_removed, 'removed_since_2014', name, True)
        if f:
            if isinstance(f, list):
                features_list.extend(f)
            else:
                features_list.append(f)

print(f"\nTotal features: {len(features_list)}")

# Build GeoJSON
geojson = {
    "type": "FeatureCollection",
    "features": features_list
}

with open('/work/vienna_boundary_changes.geojson', 'w') as f:
    json.dump(geojson, f, indent=2)

print("\nSaved vienna_boundary_changes.geojson")
print(f"Feature count: {len(features_list)}")

# Summary stats
stats = {}
for feat in features_list:
    ct = feat['properties']['change_type']
    stats[ct] = stats.get(ct, 0) + 1
print(f"By type: {stats}")

import json

with open('/work/vienna_boundary_changes.geojson') as f:
    data = json.load(f)

print(f"Total features: {len(data['features'])}")

# Check added features
for feat in data['features']:
    props = feat['properties']
    if props['change_type'] != 'unchanged':
        geom = feat['geometry']
        from shapely.geometry import shape
        s = shape(geom)
        print(f"  change_type={props['change_type']} district={props['district_name']} touches={props['touches_changed']} area={s.area:.8f}")

print("\nAll district names in output:")
names = set()
for feat in data['features']:
    names.add(feat['properties']['district_name'])
for n in sorted(names):
    print(f"  {n}")

import json
from shapely.geometry import shape, Polygon, MultiPolygon
from shapely.ops import unary_union

with open('/work/current_districts.json') as f:
    cur = json.load(f)

# Check a polygon
for el in cur['elements']:
    if el['tags'].get('ref') == '1':  # Innere Stadt
        print(f"Innere Stadt members:")
        for m in el['members']:
            if m['role'] == 'outer':
                coords = [(pt['lon'], pt['lat']) for pt in m['geometry']]
                print(f"  way {m['ref']}: {len(coords)} pts, first={coords[0]}, last={coords[-1]}")
        break

# Also check a district with changes, like Favoriten
for el in cur['elements']:
    if el['tags'].get('ref') == '10':  # Favoriten
        print(f"\nFavoriten (current) members:")
        for m in el['members']:
            if m['role'] == 'outer':
                coords = [(pt['lon'], pt['lat']) for pt in m['geometry']]
                print(f"  way {m['ref']}: {len(coords)} pts, first={coords[0]}, last={coords[-1]}")
        break

# Check 2014 Favoriten
with open('/work/districts_2014.json') as f:
    old = json.load(f)
for el in old['elements']:
    if el['tags'].get('ref') == '10':  # Favoriten
        print(f"\nFavoriten (2014) members:")
        for m in el['members']:
            if m['role'] == 'outer':
                coords = [(pt['lon'], pt['lat']) for pt in m['geometry']]
                print(f"  way {m['ref']}: {len(coords)} pts, first={coords[0]}, last={coords[-1]}")
        break

import json
from shapely.geometry import Polygon, MultiPolygon, LineString, shape, mapping
from shapely.ops import unary_union, polygonize, linemerge
import geopandas as gpd

def build_polygon_robust(rel):
    """
    Build a polygon from the relation's outer way geometries.
    Vienna districts are made of multiple ways forming closed rings.
    """
    lines = []
    for m in rel.get('members', []):
        if m['role'] == 'outer' and 'geometry' in m:
            coords = [(pt['lon'], pt['lat']) for pt in m['geometry']]
            if len(coords) >= 2:
                lines.append(LineString(coords))
    
    if not lines:
        return None
    
    merged = linemerge(lines)
    
    if merged is None or merged.is_empty:
        return None
    
    if merged.geom_type == 'LineString':
        coords = list(merged.coords)
        if coords[0] != coords[-1]:
            coords.append(coords[0])
        return Polygon(coords)
    elif merged.geom_type == 'MultiLineString':
        polys = list(polygonize(merged.geoms))
        if not polys:
            return None
        
        if len(polys) == 1:
            return polys[0]
        
        # Sort by area descending
        polys_sorted = sorted(polys, key=lambda p: p.area, reverse=True)
        outer = polys_sorted[0]
        holes = []
        for p in polys_sorted[1:]:
            if p.area > 1e-10:
                centroid = p.representative_point()
                if outer.contains(centroid) or outer.touches(p):
                    holes.append(p.exterior.coords)
        if holes:
            return Polygon(outer.exterior.coords, holes)
        return outer
    else:
        return None


def ensure_multipolygon(geom):
    """Convert to MultiPolygon if needed."""
    if geom.is_empty:
        return None
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    if geom.geom_type == 'MultiPolygon':
        return geom
    if geom.geom_type == 'GeometryCollection':
        polys = [g for g in geom.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
        if not polys:
            return None
        merged = unary_union(polys)
        return ensure_multipolygon(merged)
    return None


# Load data
print("Loading data...")
with open('/work/current_districts.json') as f:
    current_data = json.load(f)
with open('/work/districts_2014.json') as f:
    old_data = json.load(f)

# Build polygons for each district, keyed by ref number
print("Building district polygons...")
current_districts = {}  # ref -> {name, polygon}
old_districts = {}

for el in current_data['elements']:
    if el['type'] == 'relation':
        ref = el.get('tags', {}).get('ref')
        name = el.get('tags', {}).get('name', '')
        poly = build_polygon_robust(el)
        if poly and not poly.is_empty:
            current_districts[ref] = {'name': name, 'polygon': poly}
        else:
            print(f"  FAILED: current {name} (ref={ref})")

for el in old_data['elements']:
    if el['type'] == 'relation':
        ref = el.get('tags', {}).get('ref')
        name = el.get('tags', {}).get('name', '')
        poly = build_polygon_robust(el)
        if poly and not poly.is_empty:
            old_districts[ref] = {'name': name, 'polygon': poly}
        else:
            print(f"  FAILED: 2014 {name} (ref={ref})")

print(f"Current districts: {len(current_districts)}")
print(f"2014 districts: {len(old_districts)}")

# Build complete Vienna outlines
current_vienna = unary_union([current_districts[r]['polygon'] for r in current_districts])
old_vienna = unary_union([old_districts[r]['polygon'] for r in old_districts])

# Buffer slightly to handle numerical precision issues
BUF = 1e-10

# Compute the symmetric difference areas
added_total = current_vienna.difference(old_vienna)
removed_total = old_vienna.difference(current_vienna)
unchanged_total = current_vienna.intersection(old_vienna)

# Clean up small slivers
def clean_geom(geom, min_area=1e-10):
    if geom.is_empty:
        return None
    if geom.geom_type == 'MultiPolygon':
        polys = [p for p in geom.geoms if p.area >= min_area]
        if not polys:
            return None
        if len(polys) == 1:
            g = polys[0]
        else:
            g = MultiPolygon(polys)
    else:
        g = geom
    return ensure_multipolygon(g)

added_total = clean_geom(added_total)
removed_total = clean_geom(removed_total)

print(f"\nAdded area: {added_total.area:.10f}" if added_total else "Added: None")
print(f"Removed area: {removed_total.area:.10f}" if removed_total else "Removed: None")
print(f"Unchanged area: {unchanged_total.area:.10f}")

# Compute per-district properties
# For each district (by ref), get its current and 2014 polygon
# Determine which parts are added, removed, unchanged

features = []

# Check all district refs (union of both)
all_refs = sorted(set(list(current_districts.keys()) + list(old_districts.keys())), key=int)

for ref in all_refs:
    cur_info = current_districts.get(ref)
    old_info = old_districts.get(ref)
    
    # Use current name if available, else 2014 name
    if cur_info:
        name = cur_info['name']
        cur_poly = cur_info['polygon']
    else:
        name = old_info['name']
        cur_poly = None
    
    if old_info:
        old_poly = old_info['polygon']
    else:
        old_poly = None
    
    # Compute district-level changes
    district_added = None
    district_removed = None
    district_unchanged = None
    
    if cur_poly:
        # What's in this district's current area
        if added_total and not added_total.is_empty:
            district_added = cur_poly.intersection(added_total)
            district_added = clean_geom(district_added)
        district_unchanged = cur_poly.intersection(unchanged_total)
        district_unchanged = clean_geom(district_unchanged)
    
    if old_poly:
        if removed_total and not removed_total.is_empty:
            district_removed = old_poly.intersection(removed_total)
            district_removed = clean_geom(district_removed)
    
    # Does this district border any changed area?
    # Check if the district polygon (either current or old) touches added or removed areas
    touches_changed = False
    
    if cur_poly:
        if added_total and not added_total.is_empty:
            if cur_poly.touches(added_total) or cur_poly.intersects(added_total.buffer(BUF)):
                touches_changed = True
        if removed_total and not removed_total.is_empty:
            if cur_poly.touches(removed_total) or cur_poly.intersects(removed_total.buffer(BUF)):
                touches_changed = True
    
    if old_poly and not touches_changed:
        if added_total and not added_total.is_empty:
            if old_poly.touches(added_total) or old_poly.intersects(added_total.buffer(BUF)):
                touches_changed = True
        if removed_total and not removed_total.is_empty:
            if old_poly.touches(removed_total) or old_poly.intersects(removed_total.buffer(BUF)):
                touches_changed = True
    
    # Create features
    if district_unchanged:
        f = {
            "type": "Feature",
            "geometry": mapping(district_unchanged),
            "properties": {
                "change_type": "unchanged",
                "district_name": name,
                "touches_changed": touches_changed
            }
        }
        features.append(f)
    
    if district_added:
        f = {
            "type": "Feature",
            "geometry": mapping(district_added),
            "properties": {
                "change_type": "added_since_2014",
                "district_name": name,
                "touches_changed": True
            }
        }
        features.append(f)
    
    if district_removed:
        f = {
            "type": "Feature",
            "geometry": mapping(district_removed),
            "properties": {
                "change_type": "removed_since_2014",
                "district_name": name,
                "touches_changed": True
            }
        }
        features.append(f)

print(f"\nTotal features: {len(features)}")

# Build final GeoJSON
geojson = {
    "type": "FeatureCollection",
    "features": features
}

with open('/work/vienna_boundary_changes.geojson', 'w') as f:
    json.dump(geojson, f, indent=2)

print("Saved vienna_boundary_changes.geojson")

# Summary
stats = {}
for feat in features:
    ct = feat['properties']['change_type']
    stats[ct] = stats.get(ct, 0) + 1
print(f"By type: {stats}")

# Check touches_changed values
touches_stats = {}
for feat in features:
    tc = feat['properties']['touches_changed']
    touches_stats[tc] = touches_stats.get(tc, 0) + 1
print(f"By touches_changed: {touches_stats}")

# Print details
for feat in features:
    p = feat['properties']
    geom = shape(feat['geometry'])
    print(f"  {p['change_type']:20s} | {p['district_name']:25s} | touches={str(p['touches_changed']):5s} | area={geom.area:.10f}")

import json
from shapely.geometry import Polygon, MultiPolygon, LineString, shape, mapping
from shapely.ops import unary_union, polygonize, linemerge

def build_polygon(rel):
    """Build a polygon from the relation's outer way geometries."""
    lines = []
    for m in rel.get('members', []):
        if m['role'] == 'outer' and 'geometry' in m:
            coords = [(pt['lon'], pt['lat']) for pt in m['geometry']]
            if len(coords) >= 2:
                lines.append(LineString(coords))
    if not lines:
        return None
    merged = linemerge(lines)
    if merged is None or merged.is_empty:
        return None
    if merged.geom_type == 'LineString':
        coords = list(merged.coords)
        if coords[0] != coords[-1]:
            coords.append(coords[0])
        return Polygon(coords)
    elif merged.geom_type == 'MultiLineString':
        polys = list(polygonize(merged.geoms))
        if not polys:
            return None
        if len(polys) == 1:
            return polys[0]
        polys_sorted = sorted(polys, key=lambda p: p.area, reverse=True)
        outer = polys_sorted[0]
        holes = []
        for p in polys_sorted[1:]:
            if p.area > 1e-10:
                centroid = p.representative_point()
                if outer.contains(centroid) or outer.touches(p):
                    holes.append(p.exterior.coords)
        if holes:
            return Polygon(outer.exterior.coords, holes)
        return outer
    return None

def to_multipolygon(geom):
    """Convert to MultiPolygon, return None if empty."""
    if geom is None or geom.is_empty:
        return None
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    if geom.geom_type == 'MultiPolygon':
        return geom
    if geom.geom_type == 'GeometryCollection':
        polys = [g for g in geom.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
        if not polys:
            return None
        merged = unary_union(polys)
        return to_multipolygon(merged)
    return None

def clean_slivers(geom, min_area=1e-10):
    """Remove tiny sliver polygons."""
    if geom is None:
        return None
    geom = to_multipolygon(geom)
    if geom is None:
        return None
    polys = [p for p in geom.geoms if p.area >= min_area]
    if not polys:
        return None
    if len(polys) == 1:
        return MultiPolygon([polys[0]])
    return MultiPolygon(polys)

# Load data
print("Loading data...")
with open('/work/current_districts.json') as f:
    current_data = json.load(f)
with open('/work/districts_2014.json') as f:
    old_data = json.load(f)

print("Building district polygons...")
current_districts = {}
old_districts = {}

for el in current_data['elements']:
    if el['type'] == 'relation':
        ref = el.get('tags', {}).get('ref')
        name = el.get('tags', {}).get('name', '')
        poly = build_polygon(el)
        if poly:
            current_districts[ref] = {'name': name, 'polygon': poly}

for el in old_data['elements']:
    if el['type'] == 'relation':
        ref = el.get('tags', {}).get('ref')
        name = el.get('tags', {}).get('name', '')
        poly = build_polygon(el)
        if poly:
            old_districts[ref] = {'name': name, 'polygon': poly}

print(f"Current: {len(current_districts)}, 2014: {len(old_districts)}")

# Collect all refs
all_refs = sorted(set(list(current_districts.keys()) + list(old_districts.keys())), key=int)

# For each district, compute the symmetric difference between current and 2014
features = []
BUF = 1e-10

# Also compute overall changes for boundary detection
# Collect all changed areas (union of all per-district changes)
all_changed_areas = []

for ref in all_refs:
    cur_info = current_districts.get(ref)
    old_info = old_districts.get(ref)
    
    if cur_info:
        cur_name = cur_info['name']
        cur_poly = cur_info['polygon']
    else:
        cur_name = old_info['name']
        cur_poly = None
    
    if old_info:
        old_poly = old_info['polygon']
    else:
        old_poly = None
    
    if cur_poly and old_poly:
        # Normalise names - use the normalised name (current name)
        name = cur_info['name']
        
        # Compute differences
        added = cur_poly.difference(old_poly)
        removed = old_poly.difference(cur_poly)
        unchanged = cur_poly.intersection(old_poly)
        
        added = clean_slivers(added)
        removed = clean_slivers(removed)
        unchanged = clean_slivers(unchanged)
        
        if added:
            all_changed_areas.append(added)
        if removed:
            all_changed_areas.append(removed)
        
        touches = False  # Will compute later
        if added:
            features.append({
                "type": "Feature",
                "geometry": mapping(added),
                "properties": {"change_type": "added_since_2014", "district_name": name, "touches_changed": False}
            })
            touches = True
        if removed:
            features.append({
                "type": "Feature",
                "geometry": mapping(removed),
                "properties": {"change_type": "removed_since_2014", "district_name": name, "touches_changed": False}
            })
            touches = True
        if unchanged:
            features.append({
                "type": "Feature",
                "geometry": mapping(unchanged),
                "properties": {"change_type": "unchanged", "district_name": name, "touches_changed": False}
            })
    elif cur_poly and not old_poly:
        # New district (existed in current but not 2014 - unlikely for Vienna)
        name = cur_info['name']
        features.append({
            "type": "Feature",
            "geometry": mapping(to_multipolygon(cur_poly)),
            "properties": {"change_type": "added_since_2014", "district_name": name, "touches_changed": False}
        })
        all_changed_areas.append(to_multipolygon(cur_poly))
    elif old_poly and not cur_poly:
        # District removed (existed in 2014 but not current)
        name = old_info['name']
        features.append({
            "type": "Feature",
            "geometry": mapping(to_multipolygon(old_poly)),
            "properties": {"change_type": "removed_since_2014", "district_name": name, "touches_changed": False}
        })
        all_changed_areas.append(to_multipolygon(old_poly))

print(f"Features before touches: {len(features)}")

# Compute the union of all changed areas
if all_changed_areas:
    changed_union = unary_union(all_changed_areas)
    changed_union = clean_slivers(changed_union)
else:
    changed_union = None

print(f"Changed union: {changed_union.area if changed_union else 0}")

# Now set touches_changed for each feature
for feat in features:
    p = feat['properties']
    geom = shape(feat['geometry'])
    
    if p['change_type'] in ('added_since_2014', 'removed_since_2014'):
        p['touches_changed'] = True
    else:
        # unchanged - check if this geometry touches any changed area
        if changed_union and not changed_union.is_empty:
            # Buffer slightly for robustness
            touches = geom.buffer(BUF).intersects(changed_union.buffer(BUF))
            p['touches_changed'] = bool(touches)
        else:
            p['touches_changed'] = False

# Build GeoJSON
geojson = {
    "type": "FeatureCollection",
    "features": features
}

with open('/work/vienna_boundary_changes.geojson', 'w') as f:
    json.dump(geojson, f, indent=2)

print("Saved vienna_boundary_changes.geojson")

# Summary
stats = {}
for feat in features:
    ct = feat['properties']['change_type']
    stats[ct] = stats.get(ct, 0) + 1
print(f"By type: {stats}")

touches_stats = {}
for feat in features:
    tc = feat['properties']['touches_changed']
    touches_stats[tc] = touches_stats.get(tc, 0) + 1
print(f"By touches_changed: {touches_stats}")

for feat in features:
    p = feat['properties']
    g = shape(feat['geometry'])
    print(f"  {p['change_type']:20s} | {p['district_name']:25s} | touches={str(p['touches_changed']):5s} | area={g.area:.12f}")