crs-l3-tokyo-jgd-crossings

Rerun

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

1. Fetch ward boundaries from OSM via Overpass.
2. Fetch drivable roads (excl. service) in the same area.
3. Compute crossing points where roads intersect ward boundaries.
4. Buffer crossing points by 50 m (metric CRS).
5. Clip buffers to the intersected ward.
6. Compute crossing density per ward.
Output: tokyo_crossings.gpkg (5 layers).
"""

import itertools
import time
import warnings
from collections import defaultdict

import geopandas as gpd
import numpy as np
import pandas as pd
import requests
from shapely.geometry import Point, LineString, shape, MultiLineString, MultiPolygon
from shapely.ops import unary_union, linemerge
from shapely import wkb

warnings.filterwarnings("ignore")

# ---------------------------------------------------------------------------
# Constants
# ---------------------------------------------------------------------------
OVERPASS_URLS = [
    "https://overpass-api.de/api/interpreter",
    "https://overpass.kumi.systems/api/interpreter",
    "https://lz4.overpass-api.de/api/interpreter",
]

# JGD2011 / Japan Plane Rectangular CS IX (Tokyo)
METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"

OUTPUT = "/work/tokyo_crossings.gpkg"

# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------

def overpass_query(query, max_retries=3):
    """Run an Overpass QL query, return parsed JSON, retrying on failure."""
    for attempt in range(max_retries):
        for url in OVERPASS_URLS:
            try:
                r = requests.get(url, params={"data": query}, timeout=180)
                r.raise_for_status()
                data = r.json()
                if "elements" in data:
                    return data
                # Might be an error element
                if "remark" in data:
                    print(f"  Overpass remark: {data['remark']}")
                return data
            except requests.exceptions.RequestException as e:
                print(f"  Attempt {attempt+1}, URL {url} failed: {e}")
                time.sleep(2)
        time.sleep(5)
    raise RuntimeError("All Overpass retries exhausted.")


def parse_polygon_from_geometry(elem):
    """Convert an Overpass element's geometry to a Shapely geometry."""
    if elem["type"] == "way":
        coords = [(p["lon"], p["lat"]) for p in elem.get("geometry", [])]
        if len(coords) < 3:
            return None
        return LineString(coords) if coords[0] != coords[-1] else None  # not a polygon from just a way
    return None


def geom_from_element(elem):
    """Build a Shapely geometry from an Overpass element that has geometry/coords."""
    if elem["type"] == "node":
        return Point(elem["lon"], elem["lat"])
    if elem["type"] == "way":
        coords = [(p["lon"], p["lat"]) for p in elem.get("geometry", [])]
        if len(coords) < 2:
            return None
        if len(coords) == 2 or coords[0] != coords[-1]:
            return LineString(coords)
        return None
    return None


def geometry_from_element(elem, node_coords=None):
    """Build a geometry from a relation or way element."""
    if elem["type"] == "way":
        coords = [(p["lon"], p["lat"]) for p in elem.get("geometry", [])]
        if len(coords) < 2:
            return None
        if coords[0] == coords[-1] and len(coords) >= 4:
            return None  # It's a polygon ring, handle via relations
        return LineString(coords)
    return None


# ---------------------------------------------------------------------------
# Step 1: Fetch Tokyo 23 special ward boundaries (admin_level=7, boundary=administrative)
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching Tokyo 23 special ward boundaries ===")

# First, let's find the Tokyo-to relation and then the 23 special wards
# The 23 special wards in Tokyo are admin_level=7 with "特別区" in their name
# Let's query all boundary relations in Tokyo area

ward_query = """
[out:json][timeout:180];
area["name:en"="Tokyo"]["admin_level=4"]->.tokyo;
(
  relation["admin_level"="7"]["boundary"="administrative"](area.tokyo);
  relation["admin_level"="8"]["boundary"="administrative"](area.tokyo);
);
out body geom;
"""

print("  Querying Overpass for ward boundaries...")
ward_data = overpass_query(ward_query)
print(f"  Got {len(ward_data.get('elements', []))} elements")

# Filter to the 23 special wards - they have admin_level=7 and are part of Tokyo-to
# Actually, the special wards in Tokyo can be admin_level=7 or 8. Let's filter by their relation tags.
# The 23 special wards have names like "千代田区" (Chiyoda-ku), "中央区" (Chuo-ku), etc.
# Let's first collect all admin_level=7 relations since those are the special wards

relations = []
for elem in ward_data.get("elements", []):
    if elem["type"] != "relation":
        continue
    tags = elem.get("tags", {})
    name = tags.get("name", "")
    name_en = tags.get("name:en", "")
    admin_level = tags.get("admin_level", "")
    
    # Tokyo 23 special wards - they're admin_level=7
    # Also sometimes admin_level=8 for newer divisions
    # Let's filter by name pattern - ending with "区" (ku/ward)
    if admin_level == "7" and name.endswith("区"):
        # Also verify it's in Tokyo by checking the "addr:ward" or parent
        ward_type = tags.get("ward", "")
        if "特別区" in tags.get("name:ja", "") or "特別区" in tags.get("description", "") or ward_type == "特別区":
            relations.append(elem)
        elif name_en and "ku" in name_en.lower() and "Tokyo" in name_en:
            relations.append(elem)
        elif name.endswith("区"):
            # Check if it's one of the 23 special wards by known names
            known_wards = [
                "千代田区", "中央区", "港区", "新宿区", "文京区", "台東区", "墨田区",
                "江東区", "品川区", "目黒区", "大田区", "世田谷区", "渋谷区", "中野区",
                "杉並区", "豊島区", "北区", "荒川区", "板橋区", "練馬区", "足立区",
                "葛飾区", "江戸川区"
            ]
            if name in known_wards:
                relations.append(elem)

print(f"  Found {len(relations)} candidate ward relations")

# Let's get more detail if we don't have enough
if len(relations) < 23:
    print("  Not enough wards found, trying broader query...")
    # Broader query - get all admin_level=7 in Tokyo area with more detail
    ward_query2 = """
    [out:json][timeout:180];
    area["name:en"="Tokyo"]["admin_level=4"]->.tokyo;
    (
      relation["admin_level"="7"]["boundary"="administrative"](area.tokyo);
    );
    out body geom;
    """
    ward_data2 = overpass_query(ward_query2)
    relations2 = []
    for elem in ward_data2.get("elements", []):
        if elem["type"] != "relation":
            continue
        tags = elem.get("tags", {})
        name = tags.get("name", "")
        admin_level = tags.get("admin_level", "")
        if name.endswith("区"):
            relations2.append(elem)
    print(f"  Found {len(relations2)} admin_level=7 relations with '区' suffix")
    
    # If still not enough, try admin_level=8 too
    if len(relations2) < 23:
        ward_query3 = """
        [out:json][timeout:180];
        area["name:en"="Tokyo"]["admin_level=4"]->.tokyo;
        (
          relation["admin_level"]["boundary"="administrative"](area.tokyo);
        );
        out body geom;
        """
        ward_data3 = overpass_query(ward_query3)
        relations3 = []
        for elem in ward_data3.get("elements", []):
            if elem["type"] != "relation":
                continue
            tags = elem.get("tags", {})
            name = tags.get("name", "")
            name_en = tags.get("name:en", "")
            admin_level = tags.get("admin_level", "")
            # Check if it's one of the 23 special wards
            known_wards_jp = [
                "千代田区", "中央区", "港区", "新宿区", "文京区", "台東区", "墨田区",
                "江東区", "品川区", "目黒区", "大田区", "世田谷区", "渋谷区", "中野区",
                "杉並区", "豊島区", "北区", "荒川区", "板橋区", "練馬区", "足立区",
                "葛飾区", "江戸川区"
            ]
            if name in known_wards_jp:
                relations3.append(elem)
        print(f"  Found {len(relations3)} known special ward relations via broader query")
        relations = relations3 if len(relations3) >= len(relations) else relations
    else:
        relations = relations2

# Print found wards
for r in relations:
    tags = r.get("tags", {})
    print(f"  - ID={r['id']}, name={tags.get('name','')}, name:en={tags.get('name:en','')}, level={tags.get('admin_level','')}")

if len(relations) < 23:
    print(f"  WARNING: Only found {len(relations)} of 23 special wards")

# Build geometries from relation members
ward_polygons = []
ward_info = []

for rel in relations:
    tags = rel.get("tags", {})
    name = tags.get("name", "")
    name_en = tags.get("name:en", "")
    rel_id = rel["id"]
    
    # Get the outer boundary polygon
    members = rel.get("members", [])
    outer_rings = []
    inner_rings = []
    
    for mem in members:
        if mem["role"] == "outer" and mem["type"] == "way":
            geom = mem.get("geometry")
            if geom and len(geom) >= 3:
                coords = [(p["lon"], p["lat"]) for p in geom]
                if coords[0] == coords[-1]:  # Closed ring
                    outer_rings.append(coords)
        elif mem["role"] == "inner" and mem["type"] == "way":
            geom = mem.get("geometry")
            if geom and len(geom) >= 3:
                coords = [(p["lon"], p["lat"]) for p in geom]
                if coords[0] == coords[-1]:  # Closed ring
                    inner_rings.append(coords)
    
    if outer_rings:
        from shapely.geometry import Polygon
        outer = Polygon(outer_rings[0])  # Main outer
        for inner in inner_rings:
            try:
                outer = outer.difference(Polygon(inner))
            except Exception:
                pass
        for i in range(1, len(outer_rings)):
            try:
                outer = unary_union([outer, Polygon(outer_rings[i])])
            except Exception:
                pass
        
        if outer.is_valid and not outer.is_empty:
            ward_polygons.append(outer)
            ward_info.append({"ward_id": rel_id, "ward_name": name, "ward_name_en": name_en})
            print(f"  Built geometry for {name} (ID={rel_id})")
    else:
        # Try using the relation geometry directly
        if "geometry" in rel:
            geom_objs = rel["geometry"]
            from shapely.geometry import Polygon
            for g in geom_objs:
                coords = [(p["lon"], p["lat"]) for p in g]
                if len(coords) >= 4 and coords[0] == coords[-1]:
                    try:
                        poly = Polygon(coords)
                        if poly.is_valid:
                            ward_polygons.append(poly)
                            ward_info.append({"ward_id": rel_id, "ward_name": name, "ward_name_en": name_en})
                            print(f"  Built geometry for {name} from rel geom (ID={rel_id})")
                            break
                    except Exception:
                        pass

print(f"\n  Successfully built {len(ward_polygons)} ward polygons")

if len(ward_polygons) == 0:
    print("  ERROR: No ward polygons built. Trying alternative approach...")
    # Fallback: fetch ward boundaries as ways directly
    ward_query_fb = """
    [out:json][timeout:180];
    area["name:en"="Tokyo"]["admin_level=4"]->.tokyo;
    (
      rel["admin_level"="7"]["boundary"="administrative"](area.tokyo);
    );
    out body geom;
    """
    ward_data_fb = overpass_query(ward_query_fb)
    # The above gives relations with geometry in members
    # Let's try a different approach - get the full relation with all members
    for elem in ward_data_fb.get("elements", []):
        if elem["type"] != "relation":
            continue
        tags = elem.get("tags", {})
        name = tags.get("name", "")
        if not name.endswith("区"):
            continue
        
        rel_id = elem["id"]
        # Fetch the full relation with all ways
        detail_query = f"""
        [out:json][timeout:180];
        rel({rel_id});
        out body geom;
        """
        detail_data = overpass_query(detail_query)
        # Process...
        for d_elem in detail_data.get("elements", []):
            if d_elem["type"] != "relation":
                continue
            members = d_elem.get("members", [])
            outer_coords_list = []
            inner_coords_list = []
            for mem in members:
                if mem.get("geometry"):
                    coords = [(p["lon"], p["lat"]) for p in mem["geometry"]]
                    if len(coords) >= 3 and coords[0] == coords[-1]:
                        if mem["role"] == "outer":
                            outer_coords_list.append(coords)
                        elif mem["role"] == "inner":
                            inner_coords_list.append(coords)
            
            if outer_coords_list:
                from shapely.geometry import Polygon
                main_poly = Polygon(outer_coords_list[0])
                for ic in inner_coords_list:
                    try:
                        main_poly = main_poly.difference(Polygon(ic))
                    except Exception:
                        pass
                for oc in outer_coords_list[1:]:
                    try:
                        main_poly = unary_union([main_poly, Polygon(oc)])
                    except Exception:
                        pass
                
                if main_poly.is_valid and not main_poly.is_empty:
                    ward_polygons.append(main_poly)
                    name_en = tags.get("name:en", "")
                    ward_info.append({"ward_id": rel_id, "ward_name": name, "ward_name_en": name_en})
                    print(f"  Built geometry for {name} (ID={rel_id}) from detail fetch")

print(f"\n  Total ward polygons built: {len(ward_polygons)}")

if len(ward_polygons) == 0:
    print("  CRITICAL: Could not fetch ward boundaries. Exiting.")
    exit(1)

# Create GeoDataFrame
wards_gdf = gpd.GeoDataFrame(ward_info, geometry=ward_polygons, crs=WGS84)
print(f"\n  Wards DataFrame: {wards_gdf.shape}")
print(wards_gdf[["ward_id", "ward_name", "ward_name_en"]].to_string())

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Uses OSMnx for data fetching (more reliable than raw Overpass).
Output: tokyo_crossings.gpkg (5 layers).
"""

import warnings
import time
import itertools

import geopandas as gpd
import numpy as np
import pandas as pd
import osmnx as ox
from shapely.geometry import Point, LineString, MultiPoint, MultiLineString, MultiPolygon
from shapely.ops import unary_union, linemerge, split, nearest_points
from shapely import wkb

warnings.filterwarnings("ignore")

# ---------------------------------------------------------------------------
# Constants
# ---------------------------------------------------------------------------
# JGD2011 / Japan Plane Rectangular CS IX (Tokyo) — metric
METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"

OUTPUT = "/work/tokyo_crossings.gpkg"

# The 23 special wards of Tokyo
WARD_NAMES_JP = [
    "千代田区", "中央区", "港区", "新宿区", "文京区", "台東区", "墨田区",
    "江東区", "品川区", "目黒区", "大田区", "世田谷区", "渋谷区", "中野区",
    "杉並区", "豊島区", "北区", "荒川区", "板橋区", "練馬区", "足立区",
    "葛飾区", "江戸川区"
]

WARD_NAMES_EN = [
    "Chiyoda", "Chuo", "Minato", "Shinjuku", "Bunkyo", "Taito", "Sumida",
    "Koto", "Shinagawa", "Meguro", "Ota", "Setagaya", "Shibuya", "Nakano",
    "Suginami", "Toshima", "Kita", "Arakawa", "Itabashi", "Nerima", "Adachi",
    "Katsushika", "Edogawa"
]

# OSM relation IDs for the 23 special wards (annotated from OSM)
WARD_IDS = [
    1759467,  # 千代田区
    1760068,  # 中央区
    1759475,  # 港区
    1759472,  # 新宿区
    1759474,  # 文京区
    1759471,  # 台東区
    1759470,  # 墨田区
    1759469,  # 江東区
    1759477,  # 品川区
    1759478,  # 目黒区
    1759476,  # 大田区
    1759473,  # 世田谷区
    1759479,  # 渋谷区
    1759480,  # 中野区
    1759481,  # 杉並区
    1759484,  # 豊島区
    1761848,  # 北区
    1759482,  # 荒川区
    1761847,  # 板橋区
    1761846,  # 練馬区
    1761845,  # 足立区
    1761844,  # 葛飾区
    1761843,  # 江戸川区
]

# Map jp name -> en name, id
JP_TO_EN = dict(zip(WARD_NAMES_JP, WARD_NAMES_EN))
JP_TO_ID = dict(zip(WARD_NAMES_JP, WARD_IDS))


def retry_ox(max_retries=3, sleep=5):
    """Decorator to retry OSMnx calls with backoff."""
    def deco(func):
        def wrapper(*args, **kwargs):
            for attempt in range(max_retries):
                try:
                    return func(*args, **kwargs)
                except Exception as e:
                    print(f"  Attempt {attempt+1}/{max_retries} failed: {e}")
                    if attempt < max_retries - 1:
                        time.sleep(sleep * (attempt + 1))
            raise RuntimeError(f"All {max_retries} attempts for {func.__name__} failed.")
        return wrapper
    return deco


# ---------------------------------------------------------------------------
# Step 1: Fetch Tokyo 23 special ward boundaries
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching Tokyo 23 special ward boundaries ===")

wards = []
for jp_name, en_name, ward_id in zip(WARD_NAMES_JP, WARD_NAMES_EN, WARD_IDS):
    print(f"  Fetching {jp_name} ({en_name})...")
    try:
        # Use OSMnx to get the boundary
        gdf = ox.features_from_place(f"{jp_name}, Tokyo, Japan", tags={"boundary": "administrative"})
        # Filter to the relation we want
        # Look for the specific relation
        rel_mask = (gdf["element_type"] == "relation") & (gdf["osmid"] == ward_id)
        if rel_mask.any():
            row = gdf[rel_mask].iloc[0]
        else:
            # Try any relation or way that has the right name
            name_mask = gdf.get("name", "") == jp_name
            if name_mask.any() and gdf[name_mask]["element_type"].eq("relation").any():
                row = gdf[name_mask & (gdf["element_type"] == "relation")].iloc[0]
            else:
                # Check for ways with admin_level=7 or 8
                admin_mask = gdf.get("admin_level", "").isin(["7", "8"])
                name_mask = gdf.get("name", "") == jp_name
                combined = admin_mask & name_mask
                if combined.any():
                    row = gdf[combined].iloc[0]
                elif name_mask.any():
                    row = gdf[name_mask].iloc[0]
                else:
                    print(f"    WARNING: Could not find boundary for {jp_name}, trying feature query...")
                    # Broader approach
                    gdf2 = ox.features_from_place(f"{jp_name}, Tokyo, Japan", tags={"admin_level": ["7", "8"]})
                    name_mask2 = gdf2.get("name", "") == jp_name
                    if name_mask2.any():
                        row = gdf2[name_mask2].iloc[0]
                    else:
                        print(f"    WARNING: Still could not find {jp_name}, skipping")
                        continue
        
        geom = row.geometry
        if geom is None or geom.is_empty:
            print(f"    WARNING: Empty geometry for {jp_name}")
            continue
        
        # Ensure it's a polygon/multipolygon (not line from boundary way without role)
        if geom.geom_type not in ("Polygon", "MultiPolygon"):
            print(f"    WARNING: {jp_name} geometry is {geom.geom_type}, skipping")
            continue
        
        wards.append({
            "ward_id": ward_id,
            "ward_name": jp_name,
            "ward_name_en": en_name,
            "geometry": geom
        })
        print(f"    ✓ Got {jp_name}")
    except Exception as e:
        print(f"    ERROR fetching {jp_name}: {e}")
        continue

print(f"\n  Successfully fetched {len(wards)} ward boundaries")

if len(wards) < 23:
    print(f"  Only got {len(wards)} wards. Trying alternative approach for remaining...")
    # Try fetching the full Tokyo area and filtering
    try:
        tokyo_gdf = ox.features_from_place("Tokyo, Japan", tags={"boundary": "administrative", "admin_level": "7"})
        print(f"  Got {len(tokyo_gdf)} admin_level=7 features in Tokyo")
        missing_jp = [n for n in WARD_NAMES_JP if n not in [w["ward_name"] for w in wards]]
        for jp_name in missing_jp:
            mask = tokyo_gdf.get("name", "") == jp_name
            if mask.any():
                row = tokyo_gdf[mask].iloc[0]
                if row.geometry is not None and not row.geometry.is_empty:
                    geom = row.geometry
                    if geom.geom_type in ("Polygon", "MultiPolygon"):
                        ward_id = JP_TO_ID.get(jp_name, 0)
                        en_name = JP_TO_EN.get(jp_name, "")
                        if isinstance(row.get("osmid"), (list, tuple)):
                            osmid_val = row["osmid"][0] if row["osmid"] else ward_id
                        else:
                            osmid_val = row.get("osmid", ward_id)
                        wards.append({
                            "ward_id": osmid_val if osmid_val != 0 else ward_id,
                            "ward_name": jp_name,
                            "ward_name_en": en_name,
                            "geometry": geom
                        })
                        print(f"    ✓ Got {jp_name} from batch query")
    except Exception as e:
        print(f"  Batch query failed: {e}")

print(f"\n  Final ward count: {len(wards)}")

if len(wards) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    exit(1)

wards_gdf = gpd.GeoDataFrame(wards, crs=WGS84)
print(f"  Columns: {wards_gdf.columns.tolist()}")

# ---------------------------------------------------------------------------
# Step 2: Fetch drivable roads for the Tokyo 23 wards area
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads (excluding service) ===")

# Get the union boundary to restrict the road query
union_geom = unary_union(wards_gdf.geometry.values)
union_bounds = union_geom.bounds  # (minx, miny, maxx, maxy)
print(f"  Bounding box: {union_bounds}")

# Use OSMnx graph module to get drivable streets
# Filter out service roads by custom_filter
print("  Fetching drivable road network via OSMnx graph module...")
try:
    # Get drivable network for the combined area
    G = ox.graph_from_polygon(union_geom, network_type="drive", simplify=True, retain_all=True)
    print(f"  Got graph with {len(G.nodes)} nodes, {len(G.edges)} edges")
    
    # Convert edges to GeoDataFrame
    edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
    print(f"  Edges shape: {edges.shape}")
    
    # Filter out service roads
    # OSM highway types considered drivable and non-service:
    drivable_highways = {
        "motorway", "motorway_link", "trunk", "trunk_link", "primary", "primary_link",
        "secondary", "secondary_link", "tertiary", "tertiary_link",
        "unclassified", "residential", "living_street", "road"
    }
    
    # highway column might have a single value or list
    def is_drivable(hw):
        if isinstance(hw, list):
            return any(h in drivable_highways for h in hw)
        return hw in drivable_highways
    
    # Also filter by highway tag - exclude service, track, path, footway, cycleway etc.
    mask = edges["highway"].apply(is_drivable) if "highway" in edges.columns else pd.Series(True, index=edges.index)
    edges_drivable = edges[mask].copy()
    print(f"  After filtering out service/non-drivable: {len(edges_drivable)} edges")
    
    # Remove any geometry issues
    edges_drivable = edges_drivable[edges_drivable.geometry.notna() & ~edges_drivable.geometry.is_empty].copy()
    
    # Ensure all geometries are LineString
    def to_linestring(geom):
        if geom.geom_type == "MultiLineString":
            # Take the longest part
            parts = list(geom.geoms)
            if parts:
                return max(parts, key=lambda x: x.length)
        return geom
    
    edges_drivable["geometry"] = edges_drivable["geometry"].apply(to_linestring)
    edges_drivable = edges_drivable.set_geometry("geometry")
    edges_drivable = edges_drivable[edges_drivable.geometry.geom_type == "LineString"]
    print(f"  Final drivable edges: {len(edges_drivable)}")
    
except Exception as e:
    print(f"  Graph-based approach failed: {e}")
    print("  Falling back to OSMnx features query for roads...")
    # Fallback: get highway features directly
    try:
        roads = ox.features_from_polygon(union_geom, tags={"highway": True})
        print(f"  Got {len(roads)} highway features")
        
        # Filter
        def is_drivable(hw):
            if isinstance(hw, list):
                return any(h in drivable_highways for h in hw)
            return hw in drivable_highways
        
        mask = roads["highway"].apply(is_drivable) if "highway" in roads.columns else pd.Series(True, index=roads.index)
        edges_drivable = roads[mask].copy()
        
        # Filter geometry types
        edges_drivable = edges_drivable[edges_drivable.geometry.notna() & ~edges_drivable.geometry.is_empty].copy()
        mask_geom = edges_drivable.geometry.geom_type.isin(["LineString", "MultiLineString"])
        edges_drivable = edges_drivable[mask_geom].copy()
        
        def to_linestring(geom):
            if geom.geom_type == "MultiLineString":
                parts = list(geom.geoms)
                if parts:
                    return max(parts, key=lambda x: x.length)
            return geom
        
        edges_drivable["geometry"] = edges_drivable["geometry"].apply(to_linestring)
        edges_drivable = edges_drivable[edges_drivable.geometry.geom_type == "LineString"].copy()
        edges_drivable = edges_drivable.set_geometry("geometry")
        print(f"  Drivable edges (non-service): {len(edges_drivable)}")
    except Exception as e2:
        print(f"  Fallback also failed: {e2}")
        exit(1)

# ---------------------------------------------------------------------------
# Step 3: Find road/ward-boundary crossings
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing crossing points ===")

# Project to metric CRS
print("  Projecting to metric CRS...")
wards_metric = wards_gdf.to_crs(METRIC_CRS)
edges_metric = edges_drivable.to_crs(METRIC_CRS)
union_metric = unary_union(wards_metric.geometry.values)

# For each ward, get its boundary line and intersect with nearby roads
# Build a combined boundary line for all wards
print("  Building ward boundary lines...")
ward_boundaries = []
for idx, row in wards_metric.iterrows():
    geom = row.geometry
    boundary_line = geom.boundary
    if boundary_line.geom_type == "MultiLineString":
        # Split into individual lines
        for line in boundary_line.geoms:
            ward_boundaries.append({"ward_id": row["ward_id"], "geometry": line})
    else:
        ward_boundaries.append({"ward_id": row["ward_id"], "geometry": boundary_line})

boundary_gdf = gpd.GeoDataFrame(ward_boundaries, crs=METRIC_CRS, geometry="geometry")

# Spatial index for roads
edges_sindex = edges_metric.sindex

# Find crossing points - where a road line and a ward boundary intersect
print("  Finding road-boundary intersections...")
crossings = []

for idx, boundary_row in boundary_gdf.iterrows():
    boundary_line = boundary_row.geometry
    ward_id = boundary_row["ward_id"]
    
    # Find roads that might intersect this boundary segment
    possible_matches_idx = list(edges_sindex.intersection(boundary_line.bounds))
    if not possible_matches_idx:
        continue
    
    for road_idx in possible_matches_idx:
        road_geom = edges_metric.iloc[road_idx].geometry
        
        # Check intersection
        if boundary_line.intersects(road_geom):
            inter = boundary_line.intersection(road_geom)
            
            if inter.is_empty:
                continue
            
            # Handle various intersection types
            if inter.geom_type == "Point":
                crossings.append({
                    "ward_id": ward_id,
                    "geometry": inter
                })
            elif inter.geom_type == "MultiPoint":
                for pt in inter.geoms:
                    crossings.append({
                        "ward_id": ward_id,
                        "geometry": pt
                    })
            elif inter.geom_type == "GeometryCollection":
                for geom in inter.geoms:
                    if geom.geom_type == "Point":
                        crossings.append({
                            "ward_id": ward_id,
                            "geometry": geom
                        })
            # Also handle LineString intersections (road runs along boundary) - skip if it's a line
            # If we get a line, it means the road follows the boundary, not crosses it

print(f"  Found {len(crossings)} raw crossing points")

if len(crossings) == 0:
    print("  WARNING: No crossing points found!")
    # Create empty crossing layer
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64")},
        geometry=gpd.GeoSeries(dtype="geometry"),
        crs=METRIC_CRS
    )
else:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    
    # Remove duplicate points (within small tolerance)
    print("  Deduplicating crossing points...")
    # Buffer ever so slightly and union to dissolve duplicates
    crossing_gdf["x"] = crossing_gdf.geometry.x
    crossing_gdf["y"] = crossing_gdf.geometry.y
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    # Add ward names
    name_map = {}
    for _, wrow in wards_metric.iterrows():
        name_map[wrow["ward_id"]] = (wrow["ward_name"], wrow["ward_name_en"])
    
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[1])
    
print(f"  Total unique crossing points: {len(crossing_gdf)}")

# ---------------------------------------------------------------------------
# Step 4: Create 50m buffers around crossing points
# ---------------------------------------------------------------------------
print("\n=== Step 4: Creating 50m buffers ===")

if len(crossing_gdf) > 0:
    buffer_gdf = crossing_gdf.copy()
    buffer_gdf["geometry"] = buffer_gdf.geometry.buffer(50.0)
    buffer_gdf = buffer_gdf.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers")
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"),
         "ward_name": pd.Series(dtype="str"),
         "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"),
        crs=METRIC_CRS
    )
    print("  No buffers (no crossing points)")

# ---------------------------------------------------------------------------
# Step 5: Clip buffers to their respective ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to intersected ward ===")

if len(buffer_gdf) > 0:
    # For each buffer, clip to the ward it's associated with
    ward_geom_map = {}
    for _, wrow in wards_metric.iterrows():
        ward_geom_map[wrow["ward_id"]] = wrow.geometry
    
    clipped_geoms = []
    for idx, brow in buffer_gdf.iterrows():
        ward_geom = ward_geom_map.get(brow["ward_id"])
        if ward_geom is not None:
            clipped = brow.geometry.intersection(ward_geom)
            if not clipped.is_empty:
                clipped_geoms.append(clipped)
            else:
                clipped_geoms.append(None)
        else:
            clipped_geoms.append(None)
    
    clipped_gdf = buffer_gdf.copy()
    clipped_gdf["geometry"] = clipped_geoms
    clipped_gdf = clipped_gdf[clipped_gdf.geometry.notna() & ~clipped_gdf.geometry.is_empty].copy()
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    
    # Fix any invalid geometries
    clipped_gdf["geometry"] = clipped_gdf.geometry.buffer(0)
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    
    print(f"  Created {len(clipped_gdf)} clipped buffer polygons")
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"),
         "ward_name": pd.Series(dtype="str"),
         "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"),
        crs=METRIC_CRS
    )
    print("  No clipped buffers (no buffers)")

# ---------------------------------------------------------------------------
# Step 6: Compute ward crossing density (WGS84 output)
# ---------------------------------------------------------------------------
print("\n=== Step 6: Computing crossing density ===")

# Count crossings per ward
if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

# Compute ward area in km² (using metric CRS for real area)
ward_area_metric = wards_metric.copy()
ward_area_metric["ward_area_m2"] = ward_area_metric.geometry.area
ward_area_metric["ward_area_km2"] = ward_area_metric["ward_area_m2"] / 1_000_000

density_df = ward_area_metric[["ward_id", "ward_name", "ward_name_en", "ward_area_km2"]].copy()
density_df = density_df.merge(counts, on="ward_id", how="left")
density_df["crossing_count"] = density_df["crossing_count"].fillna(0).astype(int)
density_df["crossings_per_km2"] = density_df["crossing_count"] / density_df["ward_area_km2"]

# Convert back to WGS84 for the density layer
density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(density_df[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save to GeoPackage
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===")

# Layer 1: wards_jgd (in metric CRS)
wards_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wards_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ Saved wards_jgd ({len(wards_out)} features)")

# Layer 2: crossing_points (in metric CRS)
if len(crossing_gdf) > 0:
    crossing_out = crossing_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
    crossing_out.to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ Saved crossing_points ({len(crossing_out)} features)")
else:
    gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print("  ✓ Saved empty crossing_points")

# Layer 3: crossing_buffers_50m (in metric CRS)
if len(buffer_gdf) > 0:
    buffer_out = buffer_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
    buffer_out.to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ Saved crossing_buffers_50m ({len(buffer_out)} features)")
else:
    gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print("  ✓ Saved empty crossing_buffers_50m")

# Layer 4: buffer_ward_intersection (in metric CRS)
if len(clipped_gdf) > 0:
    clipped_out = clipped_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
    clipped_out.to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ Saved buffer_ward_intersection ({len(clipped_out)} features)")
else:
    gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print("  ✓ Saved empty buffer_ward_intersection")

# Layer 5: ward_crossing_density_wgs84 (in WGS84)
density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ Saved ward_crossing_density_wgs84 ({len(density_gdf)} features)")

print("\n=== DONE ===")
print(f"Output: {OUTPUT}")

# Verify
print("\n=== Verification ===")
for layer_name in ["wards_jgd", "crossing_points", "crossing_buffers_50m", "buffer_ward_intersection", "ward_crossing_density_wgs84"]:
    try:
        gdf = gpd.read_file(OUTPUT, layer=layer_name)
        print(f"  {layer_name}: {len(gdf)} features, CRS={gdf.crs}, cols={gdf.columns.tolist()}")
    except Exception as e:
        print(f"  {layer_name}: ERROR - {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Efficient single-pass approach:
  1. One Overpass query for all ward boundaries
  2. One Overpass query for all drivable roads
  3. Local geometric processing

Output: tokyo_crossings.gpkg (5 layers).
"""

import json
import time
import warnings
import requests
import sys

import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiPoint, MultiLineString, shape
from shapely.ops import unary_union, linemerge
from shapely import wkb, wkt
from pyproj import Transformer

warnings.filterwarnings("ignore")

# ---------------------------------------------------------------------------
# Constants
# ---------------------------------------------------------------------------
METRIC_CRS = "EPSG:6675"   # JGD2011 / Japan Plane Rectangular CS IX (Tokyo)
WGS84 = "EPSG:4326"

OUTPUT = "/work/tokyo_crossings.gpkg"

OVERPASS_URLS = [
    "https://overpass-api.de/api/interpreter",
    "https://overpass.kumi.systems/api/interpreter",
    "https://lz4.overpass-api.de/api/interpreter",
]

# Known 23 special ward OSM relation IDs and names
WARD_DATA = [
    (1759467, "千代田区", "Chiyoda"),
    (1760068, "中央区", "Chuo"),
    (1759475, "港区", "Minato"),
    (1759472, "新宿区", "Shinjuku"),
    (1759474, "文京区", "Bunkyo"),
    (1759471, "台東区", "Taito"),
    (1759470, "墨田区", "Sumida"),
    (1759469, "江東区", "Koto"),
    (1759477, "品川区", "Shinagawa"),
    (1759478, "目黒区", "Meguro"),
    (1759476, "大田区", "Ota"),
    (1759473, "世田谷区", "Setagaya"),
    (1759479, "渋谷区", "Shibuya"),
    (1759480, "中野区", "Nakano"),
    (1759481, "杉並区", "Suginami"),
    (1759484, "豊島区", "Toshima"),
    (1761848, "北区", "Kita"),
    (1759482, "荒川区", "Arakawa"),
    (1761847, "板橋区", "Itabashi"),
    (1761846, "練馬区", "Nerima"),
    (1761845, "足立区", "Adachi"),
    (1761844, "葛飾区", "Katsushika"),
    (1761843, "江戸川区", "Edogawa"),
]

WARD_IDS = [w[0] for w in WARD_DATA]


def run_overpass(query, max_retries=4):
    """Run an Overpass QL query, return parsed JSON."""
    for attempt in range(max_retries):
        for url in OVERPASS_URLS:
            try:
                r = requests.get(url, params={"data": query}, timeout=300)
                r.raise_for_status()
                data = r.json()
                if "elements" in data:
                    return data
                if "remark" in data:
                    print(f"  Overpass remark: {data['remark']}", file=sys.stderr)
                return data
            except requests.exceptions.RequestException as e:
                print(f"  Attempt {attempt+1}, URL {url}: {e}", file=sys.stderr)
                time.sleep(2)
        time.sleep(5)
    raise RuntimeError("All Overpass retries exhausted.")


def build_polygon_from_ring(coords_list):
    """Build a Polygon or MultiPolygon from a list of ring coordinate lists.
    coords_list: list of lists of (lon,lat) tuples. First is outer, rest are inner."""
    if not coords_list:
        return None
    outer = Polygon(coords_list[0])
    inners = coords_list[1:]
    if inners:
        for inner in inners:
            try:
                hole = Polygon(inner)
                if hole.is_valid:
                    outer = outer.difference(hole)
            except Exception:
                pass
    if outer.is_valid and not outer.is_empty:
        return outer
    return None


# ---------------------------------------------------------------------------
# Step 1: Fetch all 23 ward boundaries in one query
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching ward boundaries ===")

# Build a query that fetches each ward relation by ID
rel_ids_str = ",".join(str(wid) for wid in WARD_IDS)
ward_query = f"""
[out:json][timeout:300];
(
  rel({rel_ids_str});
);
out body geom;
"""
print("  Querying Overpass for ward relation geometries...")
ward_data = run_overpass(ward_query)
print(f"  Got {len(ward_data.get('elements',[]))} elements")

# Now fetch the full members (ways) for these relations too
# But the above with 'geom' should have the geometry on the relation members
# Let's see what we have

ward_relations = [e for e in ward_data["elements"] if e["type"] == "relation"]
print(f"  Found {len(ward_relations)} relations")

# For each relation, reconstruct the polygon from its member ways
wards_list = []

for rel in ward_relations:
    rel_id = rel["id"]
    tags = rel.get("tags", {})
    name = tags.get("name", "")
    name_en = tags.get("name:en", "")
    
    # Find match in our known list
    match = next((w for w in WARD_DATA if w[0] == rel_id), None)
    if match:
        name = name or match[1]
        name_en = name_en or match[2]
    
    # Collect outer and inner rings from members
    members = rel.get("members", [])
    outer_rings = []
    inner_rings = []
    
    for mem in members:
        role = mem.get("role", "")
        geom_coords = mem.get("geometry")
        if not geom_coords or len(geom_coords) < 3:
            continue
        
        coords = [(p["lon"], p["lat"]) for p in geom_coords]
        # Check if closed
        if coords[0] == coords[-1]:
            if role == "outer":
                outer_rings.append(coords)
            elif role == "inner":
                inner_rings.append(coords)
    
    if outer_rings:
        # Build the polygon
        poly = build_polygon_from_ring(outer_rings, inner_rings)
        if poly is None:
            # Try with just first outer ring
            poly = Polygon(outer_rings[0])
            for ic in inner_rings:
                try:
                    poly = poly.difference(Polygon(ic))
                except Exception:
                    pass
        
        if poly is not None and poly.is_valid and not poly.is_empty:
            wards_list.append({
                "ward_id": rel_id,
                "ward_name": name,
                "ward_name_en": name_en,
                "geometry": poly
            })
            print(f"  ✓ {name} ({name_en}) — {poly.area:.6f} deg²")
        else:
            print(f"  ✗ {name} — invalid geometry")
    else:
        print(f"  ✗ {name} — no outer rings found")

print(f"\n  Built {len(wards_list)} ward polygons")

# If we don't have 23, try fetching full way geometries
if len(wards_list) < 23:
    print("  Trying alternative approach: full way geometry fetch...")
    missing_ids = [w[0] for w in WARD_DATA if w[0] not in [x["ward_id"] for x in wards_list]]
    
    # Try a different query: get the relations with full way geometries using 'out geom' on the ways
    for rid in missing_ids:
        detail_query = f"""
        [out:json][timeout:300];
        rel({rid});
        (._;>>;);
        out body geom;
        """
        try:
            detail = run_overpass(detail_query)
            # Collect way geometries
            ways = {}
            for elem in detail["elements"]:
                if elem["type"] == "way":
                    coords = [(p["lon"], p["lat"]) for p in elem.get("geometry", [])]
                    ways[elem["id"]] = coords
            
            # Find the relation
            for elem in detail["elements"]:
                if elem["type"] == "relation" and elem["id"] == rid:
                    members = elem.get("members", [])
                    outer_rings = []
                    inner_rings = []
                    
                    for mem in members:
                        way_id = mem.get("ref")
                        role = mem.get("role", "")
                        if way_id in ways and len(ways[way_id]) >= 3:
                            coords = ways[way_id]
                            if coords[0] == coords[-1]:
                                if role == "outer":
                                    outer_rings.append(coords)
                                elif role == "inner":
                                    inner_rings.append(coords)
                    
                    if outer_rings:
                        poly = build_polygon_from_ring(outer_rings, inner_rings)
                        if poly is None:
                            poly = Polygon(outer_rings[0])
                            for ic in inner_rings:
                                try:
                                    poly = poly.difference(Polygon(ic))
                                except Exception:
                                    pass
                        
                        if poly is not None and poly.is_valid and not poly.is_empty:
                            tags = elem.get("tags", {})
                            name = tags.get("name", "")
                            name_en = tags.get("name:en", "")
                            match = next((w for w in WARD_DATA if w[0] == rid), None)
                            if match:
                                name = name or match[1]
                                name_en = name_en or match[2]
                            wards_list.append({
                                "ward_id": rid,
                                "ward_name": name,
                                "ward_name_en": name_en,
                                "geometry": poly
                            })
                            print(f"  ✓ {name} (via detail query)")
                    break
        except Exception as e:
            print(f"  Failed for {rid}: {e}")

print(f"\n  Final ward count: {len(wards_list)}")

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    exit(1)

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
print(f"  Wards columns: {wards_gdf.columns.tolist()}")

# ---------------------------------------------------------------------------
# Step 2: Fetch drivable roads (excluding service)
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads ===")

# Get bounding box from combined wards for the road query
union_bounds_wgs84 = unary_union(wards_gdf.geometry.values).bounds
print(f"  Bounding box: {union_bounds_wgs84}")

# Overpass query for drivable roads (excluding service)
# Get all highways in the area, we'll filter client-side
road_query = f"""
[out:json][timeout:300];
(
  way["highway"]({union_bounds_wgs84[1]},{union_bounds_wgs84[0]},{union_bounds_wgs84[3]},{union_bounds_wgs84[2]});
);
out geom;
"""
print("  Querying Overpass for roads...")
road_data = run_overpass(road_query)
print(f"  Got {len(road_data.get('elements',[]))} road elements")

# Filter and process roads
drivable_highways = {
    "motorway", "motorway_link", "trunk", "trunk_link", "primary", "primary_link",
    "secondary", "secondary_link", "tertiary", "tertiary_link",
    "unclassified", "residential", "living_street", "road"
}

road_lines = []
for elem in road_data.get("elements", []):
    if elem["type"] != "way":
        continue
    
    tags = elem.get("tags", {})
    hw = tags.get("highway", "")
    
    # Skip service roads and non-drivable types
    if hw not in drivable_highways:
        continue
    
    coords = [(p["lon"], p["lat"]) for p in elem.get("geometry", [])]
    if len(coords) < 2:
        continue
    
    road_lines.append(LineString(coords))

print(f"  Found {len(road_lines)} drivable road segments (non-service)")

if len(road_lines) == 0:
    print("CRITICAL: No roads fetched.")
    exit(1)

roads_gdf = gpd.GeoDataFrame(geometry=road_lines, crs=WGS84)
print(f"  Roads columns: {roads_gdf.columns.tolist()}")

# Project to metric
print("\n  Projecting to metric CRS...")
wards_metric = wards_gdf.to_crs(METRIC_CRS)
roads_metric = roads_gdf.to_crs(METRIC_CRS)

# ---------------------------------------------------------------------------
# Step 3: Find crossing points
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing road-boundary crossing points ===")

# Build boundary lines for each ward
boundary_lines = []  # list of (ward_id, LineString)
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_lines.append((ward_id, line))
    else:
        boundary_lines.append((ward_id, boundary))

print(f"  Built {len(boundary_lines)} boundary segments")

# Build spatial index for roads
roads_sindex = roads_metric.sindex

crossing_pts = []  # list of dicts

for ward_id, boundary_line in boundary_lines:
    # Skip very short boundary lines
    if boundary_line.length < 1:
        continue
    
    # Find candidate roads
    bounds = boundary_line.bounds
    candidates = list(roads_sindex.intersection(bounds))
    
    for ridx in candidates:
        road_geom = roads_metric.iloc[ridx].geometry
        
        if not boundary_line.intersects(road_geom):
            continue
        
        inter = boundary_line.intersection(road_geom)
        
        if inter.is_empty:
            continue
        elif inter.geom_type == "Point":
            crossing_pts.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossing_pts.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossing_pts.append({"ward_id": ward_id, "geometry": g})
        # Skip LineString intersections (road runs along boundary)

print(f"  Found {len(crossing_pts)} raw crossing points")

# Deduplicate
if crossing_pts:
    crossing_gdf = gpd.GeoDataFrame(crossing_pts, crs=METRIC_CRS, geometry="geometry")
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    # Add ward names
    name_map = {}
    for _, wrow in wards_metric.iterrows():
        name_map[wrow["ward_id"]] = (wrow["ward_name"], wrow["ward_name_en"])
    
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[1])
    
    print(f"  Deduplicated to {len(crossing_gdf)} crossing points")
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossing points found")

# ---------------------------------------------------------------------------
# Step 4: Create 50m buffers
# ---------------------------------------------------------------------------
print("\n=== Step 4: Creating 50m buffers ===")

if len(crossing_gdf) > 0:
    buffer_gdf = crossing_gdf.copy()
    buffer_gdf["geometry"] = buffer_gdf.geometry.buffer(50.0)
    buffer_gdf = buffer_gdf.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers")
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No buffers (no crossings)")

# ---------------------------------------------------------------------------
# Step 5: Clip buffers to respective ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to ward ===")

if len(buffer_gdf) > 0:
    ward_geom_map = {}
    for _, wrow in wards_metric.iterrows():
        ward_geom_map[wrow["ward_id"]] = wrow.geometry
    
    clipped_geoms = []
    for idx, brow in buffer_gdf.iterrows():
        ward_geom = ward_geom_map.get(brow["ward_id"])
        if ward_geom is not None:
            clipped = brow.geometry.intersection(ward_geom)
            clipped_geoms.append(clipped if not clipped.is_empty else None)
        else:
            clipped_geoms.append(None)
    
    clipped_gdf = buffer_gdf.copy()
    clipped_gdf["geometry"] = clipped_geoms
    clipped_gdf = clipped_gdf[clipped_gdf.geometry.notna() & ~clipped_gdf.geometry.is_empty].copy()
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    
    # Fix invalid geometries
    clipped_gdf["geometry"] = clipped_gdf.geometry.buffer(0)
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    
    print(f"  Created {len(clipped_gdf)} clipped buffer polygons")
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No clipped buffers")

# ---------------------------------------------------------------------------
# Step 6: Crossing density
# ---------------------------------------------------------------------------
print("\n=== Step 6: Crossing density ===")

# Count per ward
if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

# Ward areas in metric CRS
ward_areas = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
ward_areas["ward_area_m2"] = ward_areas.geometry.area
ward_areas["ward_area_km2"] = ward_areas["ward_area_m2"] / 1_000_000

density_df = ward_areas.merge(counts, on="ward_id", how="left")
density_df["crossing_count"] = density_df["crossing_count"].fillna(0).astype(int)
density_df["crossings_per_km2"] = density_df["crossing_count"] / density_df["ward_area_km2"]

# Convert to WGS84 for output
density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(density_df[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print("\nDensity summary:")
print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===")

# Layer 1: wards_jgd (metric)
wards_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wards_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wards_out)} features")

# Layer 2: crossing_points (metric)
if len(crossing_gdf) > 0:
    crossing_out = crossing_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
    crossing_out.to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_out)} features")
else:
    gpd.GeoDataFrame({"ward_id": [], "ward_name": [], "ward_name_en": []},
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print("  ✓ crossing_points: 0 features (empty)")

# Layer 3: crossing_buffers_50m (metric)
if len(buffer_gdf) > 0:
    buffer_out = buffer_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
    buffer_out.to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_out)} features")
else:
    gpd.GeoDataFrame({"ward_id": [], "ward_name": [], "ward_name_en": []},
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print("  ✓ crossing_buffers_50m: 0 features (empty)")

# Layer 4: buffer_ward_intersection (metric)
if len(clipped_gdf) > 0:
    clipped_out = clipped_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
    clipped_out.to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_out)} features")
else:
    gpd.GeoDataFrame({"ward_id": [], "ward_name": [], "ward_name_en": []},
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print("  ✓ buffer_ward_intersection: 0 features (empty)")

# Layer 5: ward_crossing_density_wgs84 (WGS84)
density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)} features")

print("\n=== DONE ===")

# Verify
print("\n=== Verification ===")
for layer_name in ["wards_jgd", "crossing_points", "crossing_buffers_50m", "buffer_ward_intersection", "ward_crossing_density_wgs84"]:
    try:
        gdf = gpd.read_file(OUTPUT, layer=layer_name)
        print(f"  {layer_name}: {len(gdf)} features, CRS={gdf.crs}, cols={gdf.columns.tolist()}")
        if len(gdf) > 0:
            bb = gdf.total_bounds
            print(f"    bounds: {bb}")
    except Exception as e:
        print(f"  {layer_name}: ERROR - {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Uses OSMnx per-ward for robust boundary and road fetching.
Output: tokyo_crossings.gpkg (5 layers).
"""

import time
import warnings

import geopandas as gpd
import numpy as np
import pandas as pd
import osmnx as ox
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiPoint
from shapely.ops import unary_union

warnings.filterwarnings("ignore")

# ---------------------------------------------------------------------------
# Configuration
# ---------------------------------------------------------------------------
METRIC_CRS = "EPSG:6675"   # JGD2011 / Japan Plane Rectangular CS IX (Tokyo)
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

ox.settings.timeout = 300
ox.settings.max_query_area_size = 100 * 1000 * 1000  # 100 km² sub-queries
ox.settings.bulk_download = True
ox.settings.memory = 8192

# Known 23 special ward data
WARD_DATA = [
    (1759467, "千代田区", "Chiyoda"),
    (1760068, "中央区", "Chuo"),
    (1759475, "港区", "Minato"),
    (1759472, "新宿区", "Shinjuku"),
    (1759474, "文京区", "Bunkyo"),
    (1759471, "台東区", "Taito"),
    (1759470, "墨田区", "Sumida"),
    (1759469, "江東区", "Koto"),
    (1759477, "品川区", "Shinagawa"),
    (1759478, "目黒区", "Meguro"),
    (1759476, "大田区", "Ota"),
    (1759473, "世田谷区", "Setagaya"),
    (1759479, "渋谷区", "Shibuya"),
    (1759480, "中野区", "Nakano"),
    (1759481, "杉並区", "Suginami"),
    (1759484, "豊島区", "Toshima"),
    (1761848, "北区", "Kita"),
    (1759482, "荒川区", "Arakawa"),
    (1761847, "板橋区", "Itabashi"),
    (1761846, "練馬区", "Nerima"),
    (1761845, "足立区", "Adachi"),
    (1761844, "葛飾区", "Katsushika"),
    (1761843, "江戸川区", "Edogawa"),
]


def fetch_ward_boundary(jp_name, en_name, rel_id):
    """Fetch one ward boundary polygon via OSMnx."""
    place = f"{jp_name}, Tokyo, Japan"
    
    # Try fetching the boundary with the boundary tag
    for attempt in range(3):
        try:
            # First approach: get admin boundary features
            gdf = ox.features_from_place(place, tags={"boundary": "administrative"})
            
            # The osmid is in the index (might be MultiIndex)
            # Look for a polygon that matches our ward
            candidates = []
            for idx, row in gdf.iterrows():
                geom = row.geometry
                if geom is None or geom.is_empty:
                    continue
                
                # Check if this is a relation with our ward ID
                osm_id = idx[1] if isinstance(idx, tuple) else idx
                
                # Get the name
                name = row.get("name", "")
                name_ja = row.get("name:ja", "")
                admin_level = row.get("admin_level", "")
                
                # Match criteria
                matches = False
                if isinstance(osm_id, (int, np.integer)) and osm_id == rel_id:
                    matches = True
                elif name == jp_name and admin_level in ("7", "8"):
                    matches = True
                elif name_ja == jp_name and admin_level in ("7", "8"):
                    matches = True
                
                if not matches:
                    continue
                
                if geom.geom_type in ("Polygon", "MultiPolygon"):
                    candidates.append((osm_id, geom))
            
            if candidates:
                # Pick the best one (prefer the exact relation ID)
                for osm_id, geom in candidates:
                    if isinstance(osm_id, (int, np.integer)) and osm_id == rel_id:
                        return geom
                # Otherwise take the first candidate
                return candidates[0][1]
            
            # If we got here, try alternative: get features by admin_level
            gdf2 = ox.features_from_place(place, tags={"admin_level": ["7", "8"]})
            for idx, row in gdf2.iterrows():
                geom = row.geometry
                if geom is None or geom.is_empty:
                    continue
                
                osm_id = idx[1] if isinstance(idx, tuple) else idx
                name = row.get("name", "")
                
                matches = False
                if isinstance(osm_id, (int, np.integer)) and osm_id == rel_id:
                    matches = True
                elif name == jp_name:
                    matches = True
                
                if not matches:
                    continue
                
                if geom.geom_type in ("Polygon", "MultiPolygon"):
                    return geom
            
            print(f"    No polygon found for {jp_name}, trying fallback...")
            # Last resort: use the graph module geometry
            G = ox.graph_from_place(place, network_type="drive", simplify=True)
            area_poly = ox.graph_to_gdfs(G, nodes=False, edges=False)
            if area_poly is not None and not area_poly.empty:
                # The graph_from_place returns the boundary as the GeoDataFrame's geometry
                # Actually, graph_to_gdfs only returns edges and nodes
                # Let's try getting the boundary polygon
                pass
            
        except Exception as e:
            print(f"    Attempt {attempt+1} failed: {e}")
            time.sleep(3 * (attempt + 1))
    
    raise RuntimeError(f"Could not fetch boundary for {jp_name}")


# ---------------------------------------------------------------------------
# Step 1: Fetch ward boundaries
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching ward boundaries ===")

wards_list = []
for rel_id, jp_name, en_name in WARD_DATA:
    print(f"  Fetching {jp_name} ({en_name})...")
    try:
        geom = fetch_ward_boundary(jp_name, en_name, rel_id)
        if geom is not None:
            wards_list.append({
                "ward_id": rel_id,
                "ward_name": jp_name,
                "ward_name_en": en_name,
                "geometry": geom
            })
            print(f"    ✓ {jp_name}")
    except Exception as e:
        print(f"    ✗ {jp_name}: {e}")

print(f"\n  Got {len(wards_list)} / 23 ward boundaries")

if len(wards_list) < 23:
    # Try fetching missing wards using graph_from_address
    print("\n  Attempting fallback for missing wards...")
    fetched_ids = {w["ward_id"] for w in wards_list}
    for rel_id, jp_name, en_name in WARD_DATA:
        if rel_id in fetched_ids:
            continue
        print(f"  Trying fallback for {jp_name}...")
        try:
            G = ox.graph_from_place(f"{jp_name}, Tokyo, Japan", network_type="drive", simplify=True)
            # Get the boundary from the graph's nodes
            nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True)
            # The boundary isn't directly available, but we can use the union of buffered edges
            # Actually no, this approach won't give us ward boundaries
            
            # Instead, let's try fetching the OSM relation directly
            import requests
            url = "https://overpass-api.de/api/interpreter"
            q = f"[out:json];rel({rel_id});out body geom;"
            r = requests.get(url, params={"data": q}, timeout=120)
            if r.status_code == 200:
                data = r.json()
                for elem in data.get("elements", []):
                    if elem["type"] == "relation" and elem["id"] == rel_id:
                        members = elem.get("members", [])
                        outer_rings = []
                        inner_rings = []
                        for mem in members:
                            geom_data = mem.get("geometry")
                            if not geom_data or len(geom_data) < 3:
                                continue
                            coords = [(p["lon"], p["lat"]) for p in geom_data]
                            if coords[0] == coords[-1]:
                                if mem.get("role") == "outer":
                                    outer_rings.append(coords)
                                else:
                                    inner_rings.append(coords)
                        if outer_rings:
                            poly = Polygon(outer_rings[0])
                            for ic in inner_rings:
                                try:
                                    poly = poly.difference(Polygon(ic))
                                except Exception:
                                    pass
                            if poly.is_valid and not poly.is_empty:
                                wards_list.append({
                                    "ward_id": rel_id,
                                    "ward_name": jp_name,
                                    "ward_name_en": en_name,
                                    "geometry": poly
                                })
                                print(f"    ✓ {jp_name} (via Overpass detail)")
        except Exception as e:
            print(f"    ✗ {jp_name}: {e}")

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    exit(1)

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
print(f"\n  Final wards: {len(wards_gdf)}")

# ---------------------------------------------------------------------------
# Step 2: Fetch drivable roads (excluding service)
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads ===")

# Get the union of all ward boundaries
union_geom = unary_union(wards_gdf.geometry.values)

# Use OSMnx graph module to get the drivable road network
print("  Fetching drivable road network (this may take a while for Tokyo)...")
try:
    # Use a slightly buffered version of the union to ensure we capture boundary-crossing roads
    union_buffered = union_geom.buffer(0.01)  # ~1 km buffer in degrees
    
    G = ox.graph_from_polygon(union_buffered, network_type="drive", simplify=True, retain_all=True)
    print(f"  Graph: {len(G.nodes)} nodes, {len(G.edges)} edges")
    
    # Convert edges to GeoDataFrame
    edges_gdf = ox.graph_to_gdfs(G, nodes=False, edges=True)
    print(f"  Raw edges: {len(edges_gdf)}")
    
    # The edges already exclude service roads because we used network_type="drive"
    # But let's also explicitly exclude any service-type edges
    def filter_service(hw):
        if isinstance(hw, list):
            return not any("service" in str(h).lower() for h in hw)
        return "service" not in str(hw).lower()
    
    if "highway" in edges_gdf.columns:
        mask = edges_gdf["highway"].apply(filter_service)
        edges_gdf = edges_gdf[mask].copy()
    
    # Ensure LineString geometry
    edges_gdf = edges_gdf[edges_gdf.geometry.notna() & ~edges_gdf.geometry.is_empty].copy()
    
    def to_linestring(geom):
        if geom.geom_type == "MultiLineString":
            parts = list(geom.geoms)
            return max(parts, key=lambda x: x.length) if parts else geom
        return geom
    
    edges_gdf["geometry"] = edges_gdf["geometry"].apply(to_linestring)
    edges_gdf = edges_gdf[edges_gdf.geometry.geom_type == "LineString"].copy()
    edges_gdf = edges_gdf.set_geometry("geometry")
    
    print(f"  Drivable edges (no service): {len(edges_gdf)}")
    
except Exception as e:
    print(f"  Graph approach failed: {e}")
    print("  Trying features_from_polygon fallback...")
    
    try:
        roads = ox.features_from_polygon(union_geom, tags={"highway": True})
        
        drivable_highways = {
            "motorway", "motorway_link", "trunk", "trunk_link", "primary", "primary_link",
            "secondary", "secondary_link", "tertiary", "tertiary_link",
            "unclassified", "residential", "living_street", "road"
        }
        
        def is_drivable(hw):
            if isinstance(hw, list):
                return any(h in drivable_highways for h in hw)
            return hw in drivable_highways
        
        mask = roads["highway"].apply(is_drivable) if "highway" in roads.columns else pd.Series(True, index=roads.index)
        edges_gdf = roads[mask].copy()
        edges_gdf = edges_gdf[edges_gdf.geometry.notna() & ~edges_gdf.geometry.is_empty].copy()
        
        mask_geom = edges_gdf.geometry.geom_type.isin(["LineString", "MultiLineString"])
        edges_gdf = edges_gdf[mask_geom].copy()
        
        edges_gdf["geometry"] = edges_gdf["geometry"].apply(to_linestring)
        edges_gdf = edges_gdf[edges_gdf.geometry.geom_type == "LineString"].copy()
        edges_gdf = edges_gdf.set_geometry("geometry")
        
        print(f"  Drivable edges: {len(edges_gdf)}")
    except Exception as e2:
        print(f"  Fallback also failed: {e2}")
        exit(1)

# ---------------------------------------------------------------------------
# Step 3: Find crossing points
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing crossing points ===")

# Project to metric CRS
print("  Projecting to metric CRS...")
wards_metric = wards_gdf.to_crs(METRIC_CRS)
edges_metric = edges_gdf.to_crs(METRIC_CRS)

# Build boundary lines per ward
print("  Building ward boundary lines...")
boundary_parts = []  # list of (ward_id, LineString)
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_parts.append((ward_id, line))
    elif boundary.geom_type == "LineString":
        boundary_parts.append((ward_id, boundary))

print(f"  {len(boundary_parts)} boundary segments")

# Build spatial index for roads
edges_sindex = edges_metric.sindex

# Find intersections
crossings = []
for ward_id, boundary_line in boundary_parts:
    if boundary_line.length < 1.0:
        continue
    
    bounds = boundary_line.bounds
    candidates = list(edges_sindex.intersection(bounds))
    
    for ridx in candidates:
        road_geom = edges_metric.iloc[ridx].geometry
        
        if not boundary_line.intersects(road_geom):
            continue
        
        inter = boundary_line.intersection(road_geom)
        if inter.is_empty:
            continue
        
        if inter.geom_type == "Point":
            crossings.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": ward_id, "geometry": g})

print(f"  Raw crossings: {len(crossings)}")

# Deduplicate
if crossings:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["ward_id", "x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    # Add ward names
    name_map = {}
    for _, wrow in wards_metric.iterrows():
        name_map[wrow["ward_id"]] = (wrow["ward_name"], wrow["ward_name_en"])
    
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[1])
    
    print(f"  Unique crossings: {len(crossing_gdf)}")
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings found!")

# ---------------------------------------------------------------------------
# Step 4: 50m buffers
# ---------------------------------------------------------------------------
print("\n=== Step 4: 50m buffers ===")

if len(crossing_gdf) > 0:
    buffer_gdf = crossing_gdf.copy()
    buffer_gdf["geometry"] = buffer_gdf.geometry.buffer(50.0)
    buffer_gdf = buffer_gdf.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers")
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 5: Clip buffers to ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to ward ===")

if len(buffer_gdf) > 0:
    ward_geom_map = {}
    for _, wrow in wards_metric.iterrows():
        ward_geom_map[wrow["ward_id"]] = wrow.geometry
    
    clipped_geoms = []
    for _, brow in buffer_gdf.iterrows():
        ward_geom = ward_geom_map.get(brow["ward_id"])
        if ward_geom is not None:
            clipped = brow.geometry.intersection(ward_geom)
            clipped_geoms.append(clipped if not clipped.is_empty else None)
        else:
            clipped_geoms.append(None)
    
    clipped_gdf = buffer_gdf.copy()
    clipped_gdf["geometry"] = clipped_geoms
    clipped_gdf = clipped_gdf[clipped_gdf.geometry.notna() & ~clipped_gdf.geometry.is_empty].copy()
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    clipped_gdf["geometry"] = clipped_gdf.geometry.buffer(0)  # Fix invalid
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    
    print(f"  Created {len(clipped_gdf)} clipped buffers")
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 6: Density
# ---------------------------------------------------------------------------
print("\n=== Step 6: Crossing density ===")

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

ward_areas = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
ward_areas["ward_area_km2"] = ward_areas.geometry.area / 1_000_000

density_df = ward_areas.merge(counts, on="ward_id", how="left")
density_df["crossing_count"] = density_df["crossing_count"].fillna(0).astype(int)
density_df["crossings_per_km2"] = density_df["crossing_count"] / density_df["ward_area_km2"]

# WGS84 output
density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(density_df[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===")

# Layer 1: wards_jgd
wards_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wards_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wards_out)} features")

# Layer 2: crossing_points
if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id": [], "ward_name": [], "ward_name_en": []},
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_points", driver="GPKG")

# Layer 3: crossing_buffers_50m
if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id": [], "ward_name": [], "ward_name_en": []},
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

# Layer 4: buffer_ward_intersection
if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id": [], "ward_name": [], "ward_name_en": []},
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

# Layer 5: ward_crossing_density_wgs84
density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)} features")

print("\n=== DONE ===")

# Verify
for layer in ["wards_jgd", "crossing_points", "crossing_buffers_50m", "buffer_ward_intersection", "ward_crossing_density_wgs84"]:
    try:
        gdf = gpd.read_file(OUTPUT, layer=layer)
        crs_str = str(gdf.crs).split(":")[-1] if gdf.crs else "None"
        print(f"  {layer}: {len(gdf)} feat, EPSG:{crs_str}, cols={gdf.columns.tolist()}")
    except Exception as e:
        print(f"  {layer}: {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Uses OSMnx to fetch ward boundaries and drivable roads.
Output: tokyo_crossings.gpkg (5 layers).
"""

import time
import warnings

import geopandas as gpd
import numpy as np
import pandas as pd
import osmnx as ox
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiPoint
from shapely.ops import unary_union

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

ox.settings.timeout = 300
ox.settings.max_query_area_size = 100 * 1000 * 1000
ox.settings.memory = 8192
ox.settings.bulk_download = True

# Japanese names of the 23 special wards
WARD_NAMES_JP = [
    "千代田区", "中央区", "港区", "新宿区", "文京区", "台東区", "墨田区",
    "江東区", "品川区", "目黒区", "大田区", "世田谷区", "渋谷区", "中野区",
    "杉並区", "豊島区", "北区", "荒川区", "板橋区", "練馬区", "足立区",
    "葛飾区", "江戸川区"
]

# English names in same order
WARD_NAMES_EN = [
    "Chiyoda", "Chuo", "Minato", "Shinjuku", "Bunkyo", "Taito", "Sumida",
    "Koto", "Shinagawa", "Meguro", "Ota", "Setagaya", "Shibuya", "Nakano",
    "Suginami", "Toshima", "Kita", "Arakawa", "Itabashi", "Nerima", "Adachi",
    "Katsushika", "Edogawa"
]


def get_polygon_from_place(place_name):
    """Get the ward boundary polygon from OSMnx, returning the admin_level=7 polygon."""
    try:
        gdf = ox.features_from_place(place_name, tags={"boundary": "administrative"})
    except Exception as e:
        print(f"    ox.features_from_place failed: {e}")
        return None
    
    # Find the admin_level=7 polygon with matching name
    candidates = []
    for idx, row in gdf.iterrows():
        geom = row.geometry
        if geom is None or geom.is_empty:
            continue
        if geom.geom_type not in ("Polygon", "MultiPolygon"):
            continue
        
        admin = str(row.get("admin_level", "")).strip()
        name = row.get("name", "")
        
        if admin == "7":
            candidates.append((idx, geom, name))
    
    if not candidates:
        # Try without admin_level filter
        for idx, row in gdf.iterrows():
            geom = row.geometry
            if geom is None or geom.is_empty:
                continue
            if geom.geom_type not in ("Polygon", "MultiPolygon"):
                continue
            name = row.get("name", "")
            if name == place_name.split(",")[0]:
                candidates.append((idx, geom, name))
    
    if not candidates:
        return None
    
    # Prefer the one whose name matches best
    target = place_name.split(",")[0]
    for idx, geom, name in candidates:
        if name == target:
            osm_id = idx[1] if isinstance(idx, tuple) else idx
            return (osm_id, geom)
    
    # Return first candidate
    idx, geom, name = candidates[0]
    osm_id = idx[1] if isinstance(idx, tuple) else idx
    return (osm_id, geom)


# ---------------------------------------------------------------------------
# Step 1: Fetch all 23 ward boundaries
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching ward boundaries ===")

wards_list = []
for jp_name, en_name in zip(WARD_NAMES_JP, WARD_NAMES_EN):
    place = f"{jp_name}, Tokyo, Japan"
    print(f"  Fetching {jp_name} ({en_name})...", end=" ")
    
    result = get_polygon_from_place(place)
    if result is not None:
        osm_id, geom = result
        wards_list.append({
            "ward_id": int(osm_id),
            "ward_name": jp_name,
            "ward_name_en": en_name,
            "geometry": geom
        })
        print(f"✓ (osm_id={osm_id})")
    else:
        print(f"✗ NOT FOUND")
        # Fallback: try without comma
        try:
            gdf2 = ox.features_from_place(jp_name, tags={"boundary": "administrative"})
            for idx, row in gdf2.iterrows():
                geom = row.geometry
                if geom is None or geom.is_empty:
                    continue
                if geom.geom_type not in ("Polygon", "MultiPolygon"):
                    continue
                if str(row.get("admin_level", "")).strip() == "7":
                    osm_id = idx[1] if isinstance(idx, tuple) else idx
                    wards_list.append({
                        "ward_id": int(osm_id),
                        "ward_name": jp_name,
                        "ward_name_en": en_name,
                        "geometry": geom
                    })
                    print(f"    ✓ Found via fallback (osm_id={osm_id})")
                    break
        except Exception as e:
            print(f"    Fallback failed: {e}")

print(f"\n  Got {len(wards_list)} / 23 ward boundaries")

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    exit(1)

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
print(f"  Columns: {wards_gdf.columns.tolist()}")

# ---------------------------------------------------------------------------
# Step 2: Fetch drivable roads for all wards combined
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads ===")

union_geom = unary_union(wards_gdf.geometry.values)
# Buffer slightly to catch roads that cross at boundaries
union_geom_buffered = union_geom.buffer(0.005)  # ~500m

print("  Fetching drivable road network via OSMnx graph...")
try:
    G = ox.graph_from_polygon(union_geom_buffered, network_type="drive", 
                               simplify=True, retain_all=True)
    print(f"  Graph: {len(G.nodes)} nodes, {len(G.edges)} edges")
    
    edges_gdf = ox.graph_to_gdfs(G, nodes=False, edges=True)
    print(f"  Raw edges: {len(edges_gdf)}")
    
    # Filter out any service/highway types that slipped through
    # network_type="drive" already excludes service, but double-check
    def is_drivable(hw):
        if isinstance(hw, list):
            return not any("service" in str(h).lower() for h in hw)
        return "service" not in str(hw).lower()
    
    if "highway" in edges_gdf.columns:
        mask = edges_gdf["highway"].apply(is_drivable)
        edges_gdf = edges_gdf[mask].copy()
        print(f"  After service filter: {len(edges_gdf)}")
    
    # Ensure LineString
    edges_gdf = edges_gdf[edges_gdf.geometry.notna() & ~edges_gdf.geometry.is_empty].copy()
    
    def to_linestring(geom):
        if hasattr(geom, 'geom_type') and geom.geom_type == "MultiLineString":
            parts = list(geom.geoms)
            return max(parts, key=lambda x: x.length) if parts else geom
        return geom
    
    edges_gdf["geometry"] = edges_gdf["geometry"].apply(to_linestring)
    mask = edges_gdf.geometry.geom_type == "LineString"
    edges_gdf = edges_gdf[mask].copy()
    edges_gdf = edges_gdf.set_geometry("geometry")
    
    print(f"  Final drivable edges: {len(edges_gdf)}")
    
except Exception as e:
    print(f"  Graph approach failed: {e}")
    print("  Trying features_from_polygon fallback...")
    
    try:
        roads = ox.features_from_polygon(union_geom_buffered, tags={"highway": True})
        print(f"  Got {len(roads)} highway features")
        
        drivable_highways = {
            "motorway", "motorway_link", "trunk", "trunk_link", "primary", "primary_link",
            "secondary", "secondary_link", "tertiary", "tertiary_link",
            "unclassified", "residential", "living_street", "road"
        }
        
        def is_drivable(hw):
            if isinstance(hw, list):
                return any(h in drivable_highways for h in hw)
            return hw in drivable_highways
        
        if "highway" in roads.columns:
            mask = roads["highway"].apply(is_drivable)
            roads = roads[mask].copy()
        
        roads = roads[roads.geometry.notna() & ~roads.geometry.is_empty].copy()
        mask_geom = roads.geometry.geom_type.isin(["LineString", "MultiLineString"])
        roads = roads[mask_geom].copy()
        
        roads["geometry"] = roads["geometry"].apply(to_linestring)
        roads = roads[roads.geometry.geom_type == "LineString"].copy()
        edges_gdf = roads.set_geometry("geometry")
        
        print(f"  Drivable edges: {len(edges_gdf)}")
    except Exception as e2:
        print(f"  Fallback also failed: {e2}")
        exit(1)

# ---------------------------------------------------------------------------
# Step 3: Find road-boundary crossings
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing crossing points ===")

print("  Projecting to metric CRS...")
wards_metric = wards_gdf.to_crs(METRIC_CRS)
edges_metric = edges_gdf.to_crs(METRIC_CRS)

# Build boundary lines per ward
print("  Building boundary segments...")
boundary_parts = []
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_parts.append((ward_id, line))
    elif boundary.geom_type == "LineString":
        boundary_parts.append((ward_id, boundary))

print(f"  {len(boundary_parts)} boundary segments")

# Spatial index for roads
edges_sindex = edges_metric.sindex

# Find intersections
crossings = []
for ward_id, boundary_line in boundary_parts:
    if boundary_line.length < 1.0:
        continue
    
    bounds = boundary_line.bounds
    candidates = list(edges_sindex.intersection(bounds))
    
    for ridx in candidates:
        road_geom = edges_metric.iloc[ridx].geometry
        
        if not boundary_line.intersects(road_geom):
            continue
        
        inter = boundary_line.intersection(road_geom)
        if inter.is_empty:
            continue
        
        if inter.geom_type == "Point":
            crossings.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": ward_id, "geometry": g})

print(f"  Raw crossings: {len(crossings)}")

if crossings:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    # Deduplicate
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["ward_id", "x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    # Add ward names
    name_map = {}
    for _, wrow in wards_metric.iterrows():
        name_map[wrow["ward_id"]] = (wrow["ward_name"], wrow["ward_name_en"])
    
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[1])
    
    print(f"  Unique crossings: {len(crossing_gdf)}")
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings found!")

# ---------------------------------------------------------------------------
# Step 4: 50m buffers
# ---------------------------------------------------------------------------
print("\n=== Step 4: 50m buffers ===")

if len(crossing_gdf) > 0:
    buffer_gdf = crossing_gdf.copy()
    buffer_gdf["geometry"] = buffer_gdf.geometry.buffer(50.0)
    buffer_gdf = buffer_gdf.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers")
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 5: Clip buffers to ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to ward ===")

if len(buffer_gdf) > 0:
    ward_geom_map = {}
    for _, wrow in wards_metric.iterrows():
        ward_geom_map[wrow["ward_id"]] = wrow.geometry
    
    clipped_geoms = []
    for _, brow in buffer_gdf.iterrows():
        ward_geom = ward_geom_map.get(brow["ward_id"])
        if ward_geom is not None:
            clipped = brow.geometry.intersection(ward_geom)
            clipped_geoms.append(clipped if not clipped.is_empty else None)
        else:
            clipped_geoms.append(None)
    
    clipped_gdf = buffer_gdf.copy()
    clipped_gdf["geometry"] = clipped_geoms
    clipped_gdf = clipped_gdf[clipped_gdf.geometry.notna() & ~clipped_gdf.geometry.is_empty].copy()
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    clipped_gdf["geometry"] = clipped_gdf.geometry.buffer(0)
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    
    print(f"  Created {len(clipped_gdf)} clipped buffers")
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 6: Density
# ---------------------------------------------------------------------------
print("\n=== Step 6: Crossing density ===")

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

ward_areas = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
ward_areas["ward_area_km2"] = ward_areas.geometry.area / 1_000_000

density_df = ward_areas.merge(counts, on="ward_id", how="left")
density_df["crossing_count"] = density_df["crossing_count"].fillna(0).astype(int)
density_df["crossings_per_km2"] = density_df["crossing_count"] / density_df["ward_area_km2"]

# WGS84 output
density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(density_df[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===")

# Layer 1: wards_jgd (metric)
wards_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wards_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wards_out)} features (EPSG:{METRIC_CRS.split(':')[1]})")

# Layer 2: crossing_points (metric)
if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id": [], "ward_name": [], "ward_name_en": []},
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print("  ✓ crossing_points: 0 features (empty)")

# Layer 3: crossing_buffers_50m (metric)
if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id": [], "ward_name": [], "ward_name_en": []},
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print("  ✓ crossing_buffers_50m: 0 features (empty)")

# Layer 4: buffer_ward_intersection (metric)
if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id", "ward_name", "ward_name_en", "geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id": [], "ward_name": [], "ward_name_en": []},
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print("  ✓ buffer_ward_intersection: 0 features (empty)")

# Layer 5: ward_crossing_density_wgs84 (WGS84)
density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)} features (WGS84)")

print("\n=== DONE ===")

# Verify
print("\n=== Verification ===")
for layer in ["wards_jgd", "crossing_points", "crossing_buffers_50m", "buffer_ward_intersection", "ward_crossing_density_wgs84"]:
    try:
        gdf = gpd.read_file(OUTPUT, layer=layer)
        crs_epsg = str(gdf.crs).split(":")[-1] if gdf.crs else "None"
        print(f"  {layer}: {len(gdf)} feat, EPSG:{crs_epsg}, cols={gdf.columns.tolist()}")
    except Exception as e:
        print(f"  {layer}: ERROR - {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Per-ward OSMnx fetching (avoids large-area Overpass issues).
Output: tokyo_crossings.gpkg (5 layers).
"""

import time
import warnings

import geopandas as gpd
import numpy as np
import pandas as pd
import osmnx as ox
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiPoint
from shapely.ops import unary_union
from shapely import wkb

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

ox.settings.timeout = 300
ox.settings.max_query_area_size = 50 * 1000 * 1000
ox.settings.bulk_download = True
ox.settings.memory = 4096

WARD_NAMES_JP = [
    "千代田区","中央区","港区","新宿区","文京区","台東区","墨田区",
    "江東区","品川区","目黒区","大田区","世田谷区","渋谷区","中野区",
    "杉並区","豊島区","北区","荒川区","板橋区","練馬区","足立区",
    "葛飾区","江戸川区"
]
WARD_NAMES_EN = [
    "Chiyoda","Chuo","Minato","Shinjuku","Bunkyo","Taito","Sumida",
    "Koto","Shinagawa","Meguro","Ota","Setagaya","Shibuya","Nakano",
    "Suginami","Toshima","Kita","Arakawa","Itabashi","Nerima","Adachi",
    "Katsushika","Edogawa"
]


def fetch_ward(jp_name, en_name):
    """Fetch ward boundary polygon. Returns (osm_id, geometry) or None."""
    place = f"{jp_name}, Tokyo, Japan"
    for attempt in range(3):
        try:
            gdf = ox.features_from_place(place, tags={"boundary": "administrative"})
            # Find admin_level=7 polygon
            for idx, row in gdf.iterrows():
                geom = row.geometry
                if geom is None or geom.is_empty:
                    continue
                if geom.geom_type not in ("Polygon", "MultiPolygon"):
                    continue
                if str(row.get("admin_level", "")).strip() == "7":
                    osm_id = idx[1] if isinstance(idx, tuple) else idx
                    return (int(osm_id), geom)
            # Fallback: look for polygon whose name matches
            for idx, row in gdf.iterrows():
                geom = row.geometry
                if geom is None or geom.is_empty:
                    continue
                if geom.geom_type not in ("Polygon", "MultiPolygon"):
                    continue
                name = row.get("name", "")
                if name == jp_name:
                    osm_id = idx[1] if isinstance(idx, tuple) else idx
                    return (int(osm_id), geom)
            return None
        except Exception as e:
            print(f"    Attempt {attempt+1}: {e}")
            time.sleep(3 * (attempt + 1))
    return None


def fetch_roads(ward_geom):
    """Fetch drivable roads (no service) within a ward polygon."""
    try:
        G = ox.graph_from_polygon(ward_geom, network_type="drive", 
                                   simplify=True, retain_all=True)
        edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
        if len(edges) == 0:
            return []
        
        # Ensure geometry column
        if "geometry" not in edges.columns:
            return []
        
        edges = edges[edges.geometry.notna() & ~edges.geometry.is_empty].copy()
        
        # Filter out any service/highway that slipped through
        if "highway" in edges.columns:
            def ok(hw):
                if isinstance(hw, list):
                    return not any("service" in str(h).lower() for h in hw)
                return "service" not in str(hw).lower()
            mask = edges["highway"].apply(ok)
            edges = edges[mask].copy()
        
        # Convert MultiLineString to LineString
        def to_ls(g):
            if hasattr(g, 'geom_type') and g.geom_type == "MultiLineString":
                parts = list(g.geoms)
                return max(parts, key=lambda x: x.length) if parts else g
            return g
        
        edges["geometry"] = edges["geometry"].apply(to_ls)
        edges = edges[edges.geometry.geom_type == "LineString"].copy()
        
        return list(edges.geometry.values)
    except Exception as e:
        print(f"    Road fetch failed: {e}")
        return []


# ---------------------------------------------------------------------------
# Step 1: Fetch ward boundaries (per-ward)
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching ward boundaries ===")

wards_list = []
for jp_name, en_name in zip(WARD_NAMES_JP, WARD_NAMES_EN):
    print(f"  {jp_name} ({en_name})...", end=" ", flush=True)
    result = fetch_ward(jp_name, en_name)
    if result is not None:
        osm_id, geom = result
        wards_list.append({
            "ward_id": osm_id,
            "ward_name": jp_name,
            "ward_name_en": en_name,
            "geometry": geom
        })
        print(f"✓ osm_id={osm_id}")
    else:
        print("✗ FAILED")

print(f"\n  Got {len(wards_list)} / 23 ward boundaries")

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    exit(1)

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
print(f"  Columns: {wards_gdf.columns.tolist()}")

# ---------------------------------------------------------------------------
# Step 2: Fetch drivable roads (per-ward, merge)
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads ===")

# We'll collect unique road geometries from each ward and deduplicate
all_road_geoms = []
for idx, row in wards_gdf.iterrows():
    jp_name, en_name = row["ward_name"], row["ward_name_en"]
    ward_geom = row.geometry
    # Buffer slightly inward to avoid boundary overlaps
    ward_geom_inner = ward_geom.buffer(-0.0005)  # ~-50m in degrees
    
    print(f"  Fetching roads for {jp_name} ({en_name})...", end=" ", flush=True)
    road_geoms = fetch_roads(ward_geom_inner)
    print(f"{len(road_geoms)} roads")
    all_road_geoms.extend(road_geoms)

print(f"\n  Total raw road segments: {len(all_road_geoms)}")

if len(all_road_geoms) == 0:
    print("CRITICAL: No roads fetched.")
    exit(1)

# Deduplicate road geometries (within small tolerance)
# Use WKT representation for dedup
from shapely import wkt
seen = set()
unique_roads = []
for g in all_road_geoms:
    # Round coords for dedup
    simplified = wkt.dumps(g, rounding_precision=6)
    if simplified not in seen:
        seen.add(simplified)
        unique_roads.append(g)

print(f"  Unique road segments: {len(unique_roads)}")

roads_gdf = gpd.GeoDataFrame(geometry=unique_roads, crs=WGS84)
print(f"  Roads columns: {roads_gdf.columns.tolist()}")

# ---------------------------------------------------------------------------
# Step 3: Find crossing points
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing crossing points ===")

print("  Projecting to metric CRS...")
wards_metric = wards_gdf.to_crs(METRIC_CRS)
roads_metric = roads_gdf.to_crs(METRIC_CRS)

# Build boundary line segments per ward
print("  Building boundary segments...")
boundary_parts = []
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_parts.append((ward_id, line))
    elif boundary.geom_type == "LineString":
        boundary_parts.append((ward_id, boundary))

print(f"  {len(boundary_parts)} boundary segments")

# Spatial index for roads
roads_sindex = roads_metric.sindex

# Find intersections
crossings = []
for ward_id, boundary_line in boundary_parts:
    if boundary_line.length < 1.0:
        continue
    
    bounds = boundary_line.bounds
    candidates = list(roads_sindex.intersection(bounds))
    
    for ridx in candidates:
        road_geom = roads_metric.iloc[ridx].geometry
        
        if not boundary_line.intersects(road_geom):
            continue
        
        inter = boundary_line.intersection(road_geom)
        if inter.is_empty:
            continue
        
        if inter.geom_type == "Point":
            crossings.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": ward_id, "geometry": g})

print(f"  Raw crossings: {len(crossings)}")

if crossings:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["ward_id", "x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    name_map = {}
    for _, wrow in wards_metric.iterrows():
        name_map[wrow["ward_id"]] = (wrow["ward_name"], wrow["ward_name_en"])
    
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[1])
    
    print(f"  Unique crossings: {len(crossing_gdf)}")
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings found!")

# ---------------------------------------------------------------------------
# Step 4: 50m buffers
# ---------------------------------------------------------------------------
print("\n=== Step 4: 50m buffers ===")

if len(crossing_gdf) > 0:
    buffer_gdf = crossing_gdf.copy()
    buffer_gdf["geometry"] = buffer_gdf.geometry.buffer(50.0)
    buffer_gdf = buffer_gdf.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers")
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 5: Clip buffers to ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to ward ===")

if len(buffer_gdf) > 0:
    ward_geom_map = {}
    for _, wrow in wards_metric.iterrows():
        ward_geom_map[wrow["ward_id"]] = wrow.geometry
    
    clipped_geoms = []
    for _, brow in buffer_gdf.iterrows():
        wg = ward_geom_map.get(brow["ward_id"])
        if wg is not None:
            c = brow.geometry.intersection(wg)
            clipped_geoms.append(c if not c.is_empty else None)
        else:
            clipped_geoms.append(None)
    
    clipped_gdf = buffer_gdf.copy()
    clipped_gdf["geometry"] = clipped_geoms
    clipped_gdf = clipped_gdf[clipped_gdf.geometry.notna() & ~clipped_gdf.geometry.is_empty].copy()
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    clipped_gdf["geometry"] = clipped_gdf.geometry.buffer(0)
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    print(f"  Created {len(clipped_gdf)} clipped buffers")
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 6: Density
# ---------------------------------------------------------------------------
print("\n=== Step 6: Crossing density ===")

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

ward_areas = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
ward_areas["ward_area_km2"] = ward_areas.geometry.area / 1_000_000

density_df = ward_areas.merge(counts, on="ward_id", how="left")
density_df["crossing_count"] = density_df["crossing_count"].fillna(0).astype(int)
density_df["crossings_per_km2"] = density_df["crossing_count"] / density_df["ward_area_km2"]

# WGS84 output
density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(density_df[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===")

# Layer 1
wards_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wards_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wards_out)} features")

# Layer 2
if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_gdf)}")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_points", driver="GPKG")

# Layer 3
if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_gdf)}")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

# Layer 4
if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_gdf)}")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

# Layer 5
density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)}")

print("\n=== DONE ===")

# Verify
for layer in ["wards_jgd","crossing_points","crossing_buffers_50m","buffer_ward_intersection","ward_crossing_density_wgs84"]:
    try:
        gdf = gpd.read_file(OUTPUT, layer=layer)
        crs_s = str(gdf.crs).rsplit(":",1)[-1] if gdf.crs else "None"
        print(f"  {layer}: {len(gdf)} feat, EPSG:{crs_s}, cols={gdf.columns.tolist()}")
    except Exception as e:
        print(f"  {layer}: {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Uses OSMnx with kumi Overpass mirror for reliable data fetching.
Output: tokyo_crossings.gpkg (5 layers).
"""

import time
import warnings

import geopandas as gpd
import numpy as np
import pandas as pd
import osmnx as ox
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiPoint
from shapely.ops import unary_union

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

# Use kumi mirror for reliability
ox.settings.overpass_endpoint = "https://overpass.kumi.systems/api/interpreter"
ox.settings.timeout = 300
ox.settings.max_query_area_size = 50 * 1000 * 1000
ox.settings.bulk_download = True

WARD_NAMES_JP = [
    "千代田区","中央区","港区","新宿区","文京区","台東区","墨田区",
    "江東区","品川区","目黒区","大田区","世田谷区","渋谷区","中野区",
    "杉並区","豊島区","北区","荒川区","板橋区","練馬区","足立区",
    "葛飾区","江戸川区"
]
WARD_NAMES_EN = [
    "Chiyoda","Chuo","Minato","Shinjuku","Bunkyo","Taito","Sumida",
    "Koto","Shinagawa","Meguro","Ota","Setagaya","Shibuya","Nakano",
    "Suginami","Toshima","Kita","Arakawa","Itabashi","Nerima","Adachi",
    "Katsushika","Edogawa"
]


def fetch_ward(jp_name, en_name):
    """Fetch a single ward boundary polygon. Returns (osm_id, geometry) or None."""
    place = f"{jp_name}, Tokyo, Japan"
    for attempt in range(3):
        try:
            gdf = ox.features_from_place(place, tags={"boundary": "administrative"})
            for idx, row in gdf.iterrows():
                geom = row.geometry
                if geom is None or geom.is_empty:
                    continue
                if geom.geom_type not in ("Polygon", "MultiPolygon"):
                    continue
                if str(row.get("admin_level", "")).strip() == "7":
                    osm_id = idx[1] if isinstance(idx, tuple) else idx
                    return (int(osm_id), geom)
            # Fallback: best-effort name match
            for idx, row in gdf.iterrows():
                geom = row.geometry
                if geom is None or geom.is_empty:
                    continue
                if geom.geom_type not in ("Polygon", "MultiPolygon"):
                    continue
                if row.get("name", "") == jp_name:
                    osm_id = idx[1] if isinstance(idx, tuple) else idx
                    return (int(osm_id), geom)
            return None
        except Exception as e:
            print(f"    Attempt {attempt+1}: {e}")
            time.sleep(3)
    return None


# ---------------------------------------------------------------------------
# Step 1: Fetch ward boundaries
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching ward boundaries ===")

wards_list = []
for jp_name, en_name in zip(WARD_NAMES_JP, WARD_NAMES_EN):
    print(f"  {jp_name} ({en_name})...", end=" ", flush=True)
    result = fetch_ward(jp_name, en_name)
    if result is not None:
        osm_id, geom = result
        wards_list.append({
            "ward_id": osm_id,
            "ward_name": jp_name,
            "ward_name_en": en_name,
            "geometry": geom
        })
        print(f"✓ osm_id={osm_id}")
    else:
        print("✗ FAILED")

print(f"\n  Got {len(wards_list)} / 23 ward boundaries")

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    exit(1)

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)

# Print ward IDs
for _, row in wards_gdf.iterrows():
    print(f"    {row['ward_name_en']:>12s} (ID={row['ward_id']})")

# ---------------------------------------------------------------------------
# Step 2: Fetch ALL drivable roads for Tokyo in one pass
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads ===")

# Get bounding box of all wards, buffered
union_geom = unary_union(wards_gdf.geometry.values)
bbox = union_geom.bounds  # (minx, miny, maxx, maxy)
# Buffer ~1km in each direction
buf_deg = 0.01
query_poly = Polygon([
    (bbox[0]-buf_deg, bbox[1]-buf_deg),
    (bbox[2]+buf_deg, bbox[1]-buf_deg),
    (bbox[2]+buf_deg, bbox[3]+buf_deg),
    (bbox[0]-buf_deg, bbox[3]+buf_deg),
    (bbox[0]-buf_deg, bbox[1]-buf_deg)
])

print(f"  Bounding polygon area: {query_poly.area:.2f} deg²")
print("  Fetching drivable road network (single pass)...")

for attempt in range(3):
    try:
        # Use graph_from_polygon with the full bounding box
        G = ox.graph_from_polygon(query_poly, network_type="drive",
                                   simplify=True, retain_all=True)
        print(f"  Graph: {len(G.nodes)} nodes, {len(G.edges)} edges")
        break
    except Exception as e:
        print(f"  Attempt {attempt+1} failed: {e}")
        if attempt < 2:
            print("  Retrying...")
            time.sleep(5)
        else:
            print("  Falling back to per-ward road fetching...")
            G = None

if G is None:
    # Fallback: per-ward road fetching with dedup
    print("  Fetching roads per ward...")
    all_edges = []
    for idx, row in wards_gdf.iterrows():
        jp_name, en_name = row["ward_name"], row["ward_name_en"]
        ward_geom = row.geometry.buffer(-0.0003)  # slight inward buffer
        
        print(f"    {jp_name}...", end=" ", flush=True)
        for attempt in range(2):
            try:
                Gw = ox.graph_from_polygon(ward_geom, network_type="drive",
                                            simplify=True, retain_all=True)
                edges_w = ox.graph_to_gdfs(Gw, nodes=False, edges=True)
                all_edges.append(edges_w)
                print(f"{len(edges_w)} edges")
                break
            except Exception as e:
                print(f"attempt {attempt+1} fail: {e}", end=" ")
                time.sleep(3)
        else:
            print("FAILED")
    
    if all_edges:
        edges_gdf = gpd.pd.concat(all_edges, ignore_index=True)
        # Deduplicate by geometry
        edges_gdf["geom_wkt"] = edges_gdf.geometry.apply(
            lambda g: g.wkt if hasattr(g, 'wkt') else str(g)
        )
        edges_gdf = edges_gdf.drop_duplicates(subset=["geom_wkt"])
        edges_gdf = edges_gdf.drop(columns=["geom_wkt"])
        print(f"\n  Total unique edges: {len(edges_gdf)}")
    else:
        print("CRITICAL: No roads fetched.")
        exit(1)
else:
    edges_gdf = ox.graph_to_gdfs(G, nodes=False, edges=True)
    print(f"  Raw edges: {len(edges_gdf)}")

# Filter out service roads
if "highway" in edges_gdf.columns:
    def not_service(hw):
        if isinstance(hw, list):
            return not any("service" in str(h).lower() for h in hw)
        return "service" not in str(hw).lower()
    
    mask = edges_gdf["highway"].apply(not_service)
    edges_gdf = edges_gdf[mask].copy()
    print(f"  After service filter: {len(edges_gdf)}")

# Ensure LineString geometry
edges_gdf = edges_gdf[edges_gdf.geometry.notna() & ~edges_gdf.geometry.is_empty].copy()

def to_linestring(geom):
    if hasattr(geom, 'geom_type') and geom.geom_type == "MultiLineString":
        parts = list(geom.geoms)
        return max(parts, key=lambda x: x.length) if parts else geom
    return geom

edges_gdf["geometry"] = edges_gdf["geometry"].apply(to_linestring)
edges_gdf = edges_gdf[edges_gdf.geometry.geom_type == "LineString"].copy()
edges_gdf = edges_gdf.set_geometry("geometry")

print(f"  Final drivable edges: {len(edges_gdf)}")

# ---------------------------------------------------------------------------
# Step 3: Find road-boundary crossings
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing crossing points ===")

print("  Projecting to metric CRS...")
wards_metric = wards_gdf.to_crs(METRIC_CRS)
edges_metric = edges_gdf.to_crs(METRIC_CRS)

# Clip roads to the union of ward boundaries to avoid roads outside wards
union_metric = unary_union(wards_metric.geometry.values)
print("  Clipping roads to ward union...")
edges_metric["geometry"] = edges_metric.geometry.intersection(union_metric)
edges_metric = edges_metric[~edges_metric.geometry.is_empty & 
                             (edges_metric.geometry.geom_type == "LineString")].copy()
print(f"  Roads within wards: {len(edges_metric)}")

# Build boundary line segments per ward
print("  Building boundary segments...")
boundary_parts = []
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_parts.append((ward_id, line))
    elif boundary.geom_type == "LineString":
        boundary_parts.append((ward_id, boundary))

print(f"  {len(boundary_parts)} boundary segments")

# Spatial index for roads
edges_sindex = edges_metric.sindex

# Find intersections
crossings = []
for ward_id, boundary_line in boundary_parts:
    if boundary_line.length < 1.0:
        continue
    
    bounds = boundary_line.bounds
    candidates = list(edges_sindex.intersection(bounds))
    
    for ridx in candidates:
        road_geom = edges_metric.iloc[ridx].geometry
        
        if not boundary_line.intersects(road_geom):
            continue
        
        inter = boundary_line.intersection(road_geom)
        if inter.is_empty:
            continue
        
        if inter.geom_type == "Point":
            crossings.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": ward_id, "geometry": g})

print(f"  Raw crossings: {len(crossings)}")

if crossings:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["ward_id", "x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    name_map = {}
    for _, wrow in wards_metric.iterrows():
        name_map[wrow["ward_id"]] = (wrow["ward_name"], wrow["ward_name_en"])
    
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[1])
    
    print(f"  Unique crossings: {len(crossing_gdf)}")
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings found!")

# ---------------------------------------------------------------------------
# Step 4: 50m buffers
# ---------------------------------------------------------------------------
print("\n=== Step 4: 50m buffers ===")

if len(crossing_gdf) > 0:
    buffer_gdf = crossing_gdf.copy()
    buffer_gdf["geometry"] = buffer_gdf.geometry.buffer(50.0)
    buffer_gdf = buffer_gdf.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers")
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 5: Clip buffers to ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to ward ===")

if len(buffer_gdf) > 0:
    ward_geom_map = {}
    for _, wrow in wards_metric.iterrows():
        ward_geom_map[wrow["ward_id"]] = wrow.geometry
    
    clipped_geoms = []
    for _, brow in buffer_gdf.iterrows():
        wg = ward_geom_map.get(brow["ward_id"])
        if wg is not None:
            c = brow.geometry.intersection(wg)
            clipped_geoms.append(c if not c.is_empty else None)
        else:
            clipped_geoms.append(None)
    
    clipped_gdf = buffer_gdf.copy()
    clipped_gdf["geometry"] = clipped_geoms
    clipped_gdf = clipped_gdf[clipped_gdf.geometry.notna() & ~clipped_gdf.geometry.is_empty].copy()
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    clipped_gdf["geometry"] = clipped_gdf.geometry.buffer(0)
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    print(f"  Created {len(clipped_gdf)} clipped buffers")
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 6: Density
# ---------------------------------------------------------------------------
print("\n=== Step 6: Crossing density ===")

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

ward_areas = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
ward_areas["ward_area_km2"] = ward_areas.geometry.area / 1_000_000

density_df = ward_areas.merge(counts, on="ward_id", how="left")
density_df["crossing_count"] = density_df["crossing_count"].fillna(0).astype(int)
density_df["crossings_per_km2"] = density_df["crossing_count"] / density_df["ward_area_km2"]

# Convert to WGS84
density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(density_df[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===")

# Layer 1: wards_jgd (metric)
wards_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wards_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wards_out)} features")

# Layer 2: crossing_points (metric)
if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print("  ✓ crossing_points: 0 features (empty)")

# Layer 3: crossing_buffers_50m (metric)
if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

# Layer 4: buffer_ward_intersection (metric)
if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

# Layer 5: ward_crossing_density_wgs84 (WGS84)
density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)} features")

print("\n=== DONE ===")

# Verify
print("\n=== Verification ===")
for layer in ["wards_jgd","crossing_points","crossing_buffers_50m","buffer_ward_intersection","ward_crossing_density_wgs84"]:
    try:
        gdf = gpd.read_file(OUTPUT, layer=layer)
        crs_s = str(gdf.crs).rsplit(":",1)[-1] if gdf.crs else "None"
        print(f"  {layer}: {len(gdf)} feat, EPSG:{crs_s}, cols={gdf.columns.tolist()}")
        if len(gdf) > 0 and layer == "wards_jgd":
            bb = gdf.total_bounds
            print(f"    bounds: x=[{bb[0]:.1f}, {bb[2]:.1f}] y=[{bb[1]:.1f}, {bb[3]:.1f}]")
    except Exception as e:
        print(f"  {layer}: ERROR - {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Per-ward OSMnx fetching for boundaries and roads. Uses multiple Overpass
mirrors with fallbacks. Output: tokyo_crossings.gpkg (5 layers).
"""

import sys
import time
import warnings
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiPoint
from shapely.ops import unary_union

import geopandas as gpd
import numpy as np
import pandas as pd
import osmnx as ox

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

ox.settings.timeout = 300
ox.settings.max_query_area_size = 50 * 1000 * 1000
ox.settings.bulk_download = True
ox.settings.requests_kwargs = {"timeout": 300}

WARD_NAMES_JP = [
    "千代田区","中央区","港区","新宿区","文京区","台東区","墨田区",
    "江東区","品川区","目黒区","大田区","世田谷区","渋谷区","中野区",
    "杉並区","豊島区","北区","荒川区","板橋区","練馬区","足立区",
    "葛飾区","江戸川区"
]
WARD_NAMES_EN = [
    "Chiyoda","Chuo","Minato","Shinjuku","Bunkyo","Taito","Sumida",
    "Koto","Shinagawa","Meguro","Ota","Setagaya","Shibuya","Nakano",
    "Suginami","Toshima","Kita","Arakawa","Itabashi","Nerima","Adachi",
    "Katsushika","Edogawa"
]

# Try multiple Overpass URLs
OVERPASS_URLS = [
    "https://overpass-api.de/api",
    "https://overpass.kumi.systems/api",
    "https://lz4.overpass-api.de/api",
]


def try_fetch_ward(jp_name, en_name):
    """Fetch ward boundary via OSMnx, trying multiple Overpass endpoints."""
    place = f"{jp_name}, Tokyo, Japan"
    
    for overpass_url in OVERPASS_URLS:
        ox.settings.overpass_url = overpass_url
        for attempt in range(2):
            try:
                gdf = ox.features_from_place(place, tags={"boundary": "administrative"})
                for idx, row in gdf.iterrows():
                    geom = row.geometry
                    if geom is None or geom.is_empty:
                        continue
                    if geom.geom_type not in ("Polygon", "MultiPolygon"):
                        continue
                    if str(row.get("admin_level", "")).strip() == "7":
                        osm_id = idx[1] if isinstance(idx, tuple) else idx
                        return (int(osm_id), geom)
                # Fallback: any polygon with matching name
                for idx, row in gdf.iterrows():
                    geom = row.geometry
                    if geom is None or geom.is_empty:
                        continue
                    if geom.geom_type not in ("Polygon", "MultiPolygon"):
                        continue
                    if row.get("name", "") == jp_name:
                        osm_id = idx[1] if isinstance(idx, tuple) else idx
                        return (int(osm_id), geom)
            except Exception as e:
                err = str(e)
                if "No data" in err or "Connection refused" in err or "timed out" in err:
                    break  # Try next URL
                time.sleep(2)
    return None


def try_fetch_roads(polygon, label=""):
    """Fetch drivable roads within a polygon via OSMnx graph module."""
    for overpass_url in OVERPASS_URLS:
        ox.settings.overpass_url = overpass_url
        for attempt in range(2):
            try:
                G = ox.graph_from_polygon(polygon, network_type="drive",
                                           simplify=True, retain_all=True)
                edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
                if len(edges) > 0:
                    return edges
            except Exception as e:
                err = str(e)
                if "Connection refused" in err or "timed out" in err:
                    break
                time.sleep(2)
    return None


# ---------------------------------------------------------------------------
# Step 1: Ward boundaries
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching ward boundaries ===", flush=True)

wards_list = []
for jp_name, en_name in zip(WARD_NAMES_JP, WARD_NAMES_EN):
    print(f"  {jp_name} ({en_name})...", end=" ", flush=True)
    result = try_fetch_ward(jp_name, en_name)
    if result is not None:
        osm_id, geom = result
        wards_list.append({
            "ward_id": osm_id,
            "ward_name": jp_name,
            "ward_name_en": en_name,
            "geometry": geom
        })
        print(f"✓ osm_id={osm_id}")
    else:
        print("✗ FAILED")

print(f"\n  Got {len(wards_list)} / 23 ward boundaries", flush=True)

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    sys.exit(1)

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
for _, row in wards_gdf.iterrows():
    print(f"    {row['ward_name_en']:>12s} (ID={row['ward_id']})")

# ---------------------------------------------------------------------------
# Step 2: Fetch roads per-ward (with dedup)
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads (per-ward) ===", flush=True)

all_road_geometries = []
for idx, row in wards_gdf.iterrows():
    jp_name = row["ward_name"]
    en_name = row["ward_name_en"]
    print(f"  {jp_name} ({en_name})...", end=" ", flush=True)
    
    # Buffer ward slightly inward to avoid boundary overlap issues
    inner_poly = row.geometry.buffer(-0.0003)  # ~-30m in degrees
    
    edges = try_fetch_roads(inner_poly, f"{jp_name}")
    if edges is not None:
        # Filter out service roads
        if "highway" in edges.columns:
            def not_service(hw):
                if isinstance(hw, list):
                    return not any("service" in str(h).lower() for h in hw)
                return "service" not in str(hw).lower()
            mask = edges["highway"].apply(not_service)
            edges = edges[mask].copy()
        
        # Ensure LineString
        edges = edges[edges.geometry.notna() & ~edges.geometry.is_empty].copy()
        
        def to_ls(g):
            if hasattr(g, 'geom_type') and g.geom_type == "MultiLineString":
                parts = list(g.geoms)
                return max(parts, key=lambda x: x.length) if parts else g
            return g
        
        edges["geometry"] = edges["geometry"].apply(to_ls)
        edges = edges[edges.geometry.geom_type == "LineString"].copy()
        
        all_road_geometries.extend(edges.geometry.values)
        print(f"{len(edges)} roads")
    else:
        print("FAILED (no roads)")

print(f"\n  Total raw road segments: {len(all_road_geometries)}", flush=True)

if len(all_road_geometries) == 0:
    print("CRITICAL: No roads fetched.")
    sys.exit(1)

# Deduplicate road geometries
from shapely import wkt
seen = set()
unique_roads = []
for g in all_road_geometries:
    key = wkt.dumps(g, rounding_precision=6)
    if key not in seen:
        seen.add(key)
        unique_roads.append(g)

print(f"  Unique road segments: {len(unique_roads)}", flush=True)

roads_gdf = gpd.GeoDataFrame(geometry=unique_roads, crs=WGS84)

# ---------------------------------------------------------------------------
# Step 3: Find crossings
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing crossing points ===", flush=True)

print("  Projecting to metric CRS...", flush=True)
wards_metric = wards_gdf.to_crs(METRIC_CRS)
roads_metric = roads_gdf.to_crs(METRIC_CRS)

# Clip roads to union of ward boundaries (we buffered inward, so extend back)
union_metric = unary_union(wards_metric.geometry.values)
print("  Clipping roads to ward union...", flush=True)
roads_metric["geometry"] = roads_metric.geometry.intersection(union_metric)
roads_metric = roads_metric[~roads_metric.geometry.is_empty & 
                             (roads_metric.geometry.geom_type == "LineString")].copy()
print(f"  Roads within wards: {len(roads_metric)}", flush=True)

# Build boundary line segments per ward
print("  Building boundary segments...", flush=True)
boundary_parts = []
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_parts.append((ward_id, line))
    elif boundary.geom_type == "LineString":
        boundary_parts.append((ward_id, boundary))

print(f"  {len(boundary_parts)} boundary segments", flush=True)

# Spatial index
roads_sindex = roads_metric.sindex

# Find intersections
crossings = []
for ward_id, boundary_line in boundary_parts:
    if boundary_line.length < 1.0:
        continue
    
    bounds = boundary_line.bounds
    candidates = list(roads_sindex.intersection(bounds))
    
    for ridx in candidates:
        road_geom = roads_metric.iloc[ridx].geometry
        
        if not boundary_line.intersects(road_geom):
            continue
        
        inter = boundary_line.intersection(road_geom)
        if inter.is_empty:
            continue
        
        if inter.geom_type == "Point":
            crossings.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": ward_id, "geometry": g})

print(f"  Raw crossings: {len(crossings)}", flush=True)

if crossings:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["ward_id", "x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    name_map = {}
    for _, wrow in wards_metric.iterrows():
        name_map[wrow["ward_id"]] = (wrow["ward_name"], wrow["ward_name_en"])
    
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[1])
    
    print(f"  Unique crossings: {len(crossing_gdf)}", flush=True)
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings found!", flush=True)

# ---------------------------------------------------------------------------
# Step 4: 50m buffers
# ---------------------------------------------------------------------------
print("\n=== Step 4: 50m buffers ===", flush=True)

if len(crossing_gdf) > 0:
    buffer_gdf = crossing_gdf.copy()
    buffer_gdf["geometry"] = buffer_gdf.geometry.buffer(50.0)
    buffer_gdf = buffer_gdf.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers", flush=True)
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 5: Clip buffers to ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to ward ===", flush=True)

if len(buffer_gdf) > 0:
    ward_geom_map = {}
    for _, wrow in wards_metric.iterrows():
        ward_geom_map[wrow["ward_id"]] = wrow.geometry
    
    clipped_geoms = []
    for _, brow in buffer_gdf.iterrows():
        wg = ward_geom_map.get(brow["ward_id"])
        if wg is not None:
            c = brow.geometry.intersection(wg)
            clipped_geoms.append(c if not c.is_empty else None)
        else:
            clipped_geoms.append(None)
    
    clipped_gdf = buffer_gdf.copy()
    clipped_gdf["geometry"] = clipped_geoms
    clipped_gdf = clipped_gdf[clipped_gdf.geometry.notna() & ~clipped_gdf.geometry.is_empty].copy()
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    clipped_gdf["geometry"] = clipped_gdf.geometry.buffer(0)
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    print(f"  Created {len(clipped_gdf)} clipped buffers", flush=True)
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 6: Density
# ---------------------------------------------------------------------------
print("\n=== Step 6: Crossing density ===", flush=True)

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

ward_areas = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
ward_areas["ward_area_km2"] = ward_areas.geometry.area / 1_000_000

density_df = ward_areas.merge(counts, on="ward_id", how="left")
density_df["crossing_count"] = density_df["crossing_count"].fillna(0).astype(int)
density_df["crossings_per_km2"] = density_df["crossing_count"] / density_df["ward_area_km2"]

density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(density_df[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===", flush=True)

wards_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wards_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wards_out)} features")

if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print("  ✓ crossing_points: 0 features (empty)")

if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)} features")

print("\n=== DONE ===", flush=True)

# Verify
print("\n=== Verification ===")
for layer in ["wards_jgd","crossing_points","crossing_buffers_50m","buffer_ward_intersection","ward_crossing_density_wgs84"]:
    try:
        gdf = gpd.read_file(OUTPUT, layer=layer)
        crs_s = str(gdf.crs).rsplit(":",1)[-1] if gdf.crs else "None"
        print(f"  {layer}: {len(gdf)} feat, EPSG:{crs_s}, cols={gdf.columns.tolist()}")
    except Exception as e:
        print(f"  {layer}: ERROR - {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Per-ward OSMnx fetching. Output: tokyo_crossings.gpkg (5 layers).
"""

import sys
import time
import warnings
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiPoint
from shapely.ops import unary_union

import geopandas as gpd
import numpy as np
import pandas as pd
import osmnx as ox

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

ox.settings.timeout = 300
ox.settings.max_query_area_size = 50 * 1000 * 1000
ox.settings.bulk_download = True
ox.settings.requests_kwargs = {"timeout": 300}

# Exact Japanese names → English names mapping
WARD_MAP = {
    "千代田区": "Chiyoda",
    "中央区": "Chuo",
    "港区": "Minato",
    "新宿区": "Shinjuku",
    "文京区": "Bunkyo",
    "台東区": "Taito",
    "墨田区": "Sumida",
    "江東区": "Koto",
    "品川区": "Shinagawa",
    "目黒区": "Meguro",
    "大田区": "Ota",
    "世田谷区": "Setagaya",
    "渋谷区": "Shibuya",
    "中野区": "Nakano",
    "杉並区": "Suginami",
    "豊島区": "Toshima",
    "北区": "Kita",
    "荒川区": "Arakawa",
    "板橋区": "Itabashi",
    "練馬区": "Nerima",
    "足立区": "Adachi",
    "葛飾区": "Katsushika",
    "江戸川区": "Edogawa",
}

OVERPASS_URLS = [
    "https://overpass-api.de/api",
    "https://overpass.kumi.systems/api",
    "https://lz4.overpass-api.de/api",
]


def try_fetch_ward(jp_name, en_name):
    """Fetch ward boundary polygon, ensuring the name matches."""
    place = f"{jp_name}, Tokyo, Japan"
    
    for overpass_url in OVERPASS_URLS:
        ox.settings.overpass_url = overpass_url
        for attempt in range(2):
            try:
                gdf = ox.features_from_place(place, tags={"boundary": "administrative"})
                
                candidates = []
                for idx, row in gdf.iterrows():
                    geom = row.geometry
                    if geom is None or geom.is_empty:
                        continue
                    if geom.geom_type not in ("Polygon", "MultiPolygon"):
                        continue
                    
                    name = row.get("name", "")
                    name_ja = row.get("name:ja", "")
                    admin = str(row.get("admin_level", "")).strip()
                    
                    # Must match the Japanese name
                    name_match = (name == jp_name) or (name_ja == jp_name)
                    if not name_match:
                        continue
                    
                    osm_id = idx[1] if isinstance(idx, tuple) else idx
                    candidates.append((int(osm_id), geom, admin))
                
                if candidates:
                    # Prefer admin_level=7
                    for osm_id, geom, admin in candidates:
                        if admin == "7":
                            return (osm_id, geom)
                    return (candidates[0][0], candidates[0][1])
                    
            except Exception as e:
                err = str(e)
                if "Connection refused" in err or "timed out" in err:
                    break  # Try next URL
                time.sleep(2)
    return None


def try_fetch_roads(polygon):
    """Fetch drivable roads within a polygon. Returns list of LineString geometries."""
    for overpass_url in OVERPASS_URLS:
        ox.settings.overpass_url = overpass_url
        for attempt in range(2):
            try:
                G = ox.graph_from_polygon(polygon, network_type="drive",
                                           simplify=True, retain_all=True)
                edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
                if edges is None or len(edges) == 0:
                    return []
                
                if "geometry" not in edges.columns:
                    return []
                
                edges = edges[edges.geometry.notna() & ~edges.geometry.is_empty].copy()
                
                # Filter service roads
                if "highway" in edges.columns:
                    def not_service(hw):
                        if isinstance(hw, list):
                            return not any("service" in str(h).lower() for h in hw)
                        return "service" not in str(hw).lower()
                    mask = edges["highway"].apply(not_service)
                    edges = edges[mask].copy()
                
                # Normalize to LineString
                def to_ls(g):
                    if hasattr(g, 'geom_type') and g.geom_type == "MultiLineString":
                        parts = list(g.geoms)
                        return max(parts, key=lambda x: x.length) if parts else g
                    return g
                
                edges["geometry"] = edges["geometry"].apply(to_ls)
                edges = edges[edges.geometry.geom_type == "LineString"].copy()
                
                return list(edges.geometry.values)
                
            except Exception as e:
                err = str(e)
                if "Connection refused" in err or "timed out" in err:
                    break
                time.sleep(2)
    return []


# ---------------------------------------------------------------------------
# Step 1: Ward boundaries
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching ward boundaries ===", flush=True)

wards_list = []
for jp_name, en_name in WARD_MAP.items():
    print(f"  {jp_name} ({en_name})...", end=" ", flush=True)
    result = try_fetch_ward(jp_name, en_name)
    if result is not None:
        osm_id, geom = result
        wards_list.append({
            "ward_id": osm_id,
            "ward_name": jp_name,
            "ward_name_en": en_name,
            "geometry": geom
        })
        print(f"✓ osm_id={osm_id}")
    else:
        print("✗ FAILED")

print(f"\n  Got {len(wards_list)} / 23 ward boundaries", flush=True)

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    sys.exit(1)

# Check for duplicates
ward_ids = [w["ward_id"] for w in wards_list]
if len(set(ward_ids)) < len(ward_ids):
    print("  WARNING: Duplicate ward IDs detected!")
    from collections import Counter
    for wid, count in Counter(ward_ids).items():
        if count > 1:
            names = [w["ward_name"] for w in wards_list if w["ward_id"] == wid]
            print(f"    ID {wid} appears {count}x: {names}")

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
for _, row in wards_gdf.iterrows():
    print(f"    {row['ward_name_en']:>12s} (ID={row['ward_id']})")

# ---------------------------------------------------------------------------
# Step 2: Fetch roads per-ward
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads (per-ward) ===", flush=True)

all_road_geometries = []
for idx, row in wards_gdf.iterrows():
    jp_name = row["ward_name"]
    en_name = row["ward_name_en"]
    print(f"  {jp_name} ({en_name})...", end=" ", flush=True)
    
    # Use the full ward polygon (no inward buffer) since we need roads up to boundary
    ward_poly = row.geometry
    
    road_geoms = try_fetch_roads(ward_poly)
    print(f"{len(road_geoms)} roads")
    all_road_geometries.extend(road_geoms)

print(f"\n  Total raw road segments: {len(all_road_geometries)}", flush=True)

if len(all_road_geometries) == 0:
    print("CRITICAL: No roads fetched.")
    sys.exit(1)

# Deduplicate
from shapely import wkt
seen = set()
unique_roads = []
for g in all_road_geometries:
    key = wkt.dumps(g, rounding_precision=6)
    if key not in seen:
        seen.add(key)
        unique_roads.append(g)

print(f"  Unique road segments: {len(unique_roads)}", flush=True)

roads_gdf = gpd.GeoDataFrame(geometry=unique_roads, crs=WGS84)

# ---------------------------------------------------------------------------
# Step 3: Crossing points
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing crossing points ===", flush=True)

print("  Projecting to metric CRS...", flush=True)
wards_metric = wards_gdf.to_crs(METRIC_CRS)
roads_metric = roads_gdf.to_crs(METRIC_CRS)

# Since we fetched roads per-ward using the full polygon, roads should already
# be within ward boundaries. But clip to be safe.
union_metric = unary_union(wards_metric.geometry.values)
print("  Clipping roads to ward union...", flush=True)
roads_metric["geometry"] = roads_metric.geometry.intersection(union_metric)
roads_metric = roads_metric[~roads_metric.geometry.is_empty & 
                             (roads_metric.geometry.geom_type == "LineString")].copy()
print(f"  Roads within wards: {len(roads_metric)}", flush=True)

# Build boundary line segments per ward
print("  Building boundary segments...", flush=True)
boundary_parts = []
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_parts.append((ward_id, line))
    elif boundary.geom_type == "LineString":
        boundary_parts.append((ward_id, boundary))

print(f"  {len(boundary_parts)} boundary segments", flush=True)

# Spatial index
roads_sindex = roads_metric.sindex

# Find intersections
crossings = []
for ward_id, boundary_line in boundary_parts:
    if boundary_line.length < 1.0:
        continue
    
    bounds = boundary_line.bounds
    candidates = list(roads_sindex.intersection(bounds))
    
    for ridx in candidates:
        road_geom = roads_metric.iloc[ridx].geometry
        
        if not boundary_line.intersects(road_geom):
            continue
        
        inter = boundary_line.intersection(road_geom)
        if inter.is_empty:
            continue
        
        if inter.geom_type == "Point":
            crossings.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": ward_id, "geometry": g})

print(f"  Raw crossings: {len(crossings)}", flush=True)

if crossings:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["ward_id", "x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    name_map = {}
    for _, wrow in wards_metric.iterrows():
        name_map[wrow["ward_id"]] = (wrow["ward_name"], wrow["ward_name_en"])
    
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[1])
    
    print(f"  Unique crossings: {len(crossing_gdf)}", flush=True)
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings found!", flush=True)

# ---------------------------------------------------------------------------
# Step 4: 50m buffers
# ---------------------------------------------------------------------------
print("\n=== Step 4: 50m buffers ===", flush=True)

if len(crossing_gdf) > 0:
    buffer_gdf = crossing_gdf.copy()
    buffer_gdf["geometry"] = buffer_gdf.geometry.buffer(50.0)
    buffer_gdf = buffer_gdf.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers", flush=True)
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 5: Clip buffers to ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to ward ===", flush=True)

if len(buffer_gdf) > 0:
    ward_geom_map = {}
    for _, wrow in wards_metric.iterrows():
        ward_geom_map[wrow["ward_id"]] = wrow.geometry
    
    clipped_geoms = []
    for _, brow in buffer_gdf.iterrows():
        wg = ward_geom_map.get(brow["ward_id"])
        if wg is not None:
            c = brow.geometry.intersection(wg)
            clipped_geoms.append(c if not c.is_empty else None)
        else:
            clipped_geoms.append(None)
    
    clipped_gdf = buffer_gdf.copy()
    clipped_gdf["geometry"] = clipped_geoms
    clipped_gdf = clipped_gdf[clipped_gdf.geometry.notna() & ~clipped_gdf.geometry.is_empty].copy()
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    clipped_gdf["geometry"] = clipped_gdf.geometry.buffer(0)
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    print(f"  Created {len(clipped_gdf)} clipped buffers", flush=True)
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 6: Density
# ---------------------------------------------------------------------------
print("\n=== Step 6: Crossing density ===", flush=True)

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

ward_areas = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
ward_areas["ward_area_km2"] = ward_areas.geometry.area / 1_000_000

density_df = ward_areas.merge(counts, on="ward_id", how="left")
density_df["crossing_count"] = density_df["crossing_count"].fillna(0).astype(int)
density_df["crossings_per_km2"] = density_df["crossing_count"] / density_df["ward_area_km2"]

density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(density_df[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===", flush=True)

wards_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wards_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wards_out)} features")

if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print("  ✓ crossing_points: 0 features (empty)")

if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)} features")

print("\n=== DONE ===", flush=True)

# Verify
print("\n=== Verification ===")
for layer in ["wards_jgd","crossing_points","crossing_buffers_50m","buffer_ward_intersection","ward_crossing_density_wgs84"]:
    try:
        gdf = gpd.read_file(OUTPUT, layer=layer)
        crs_s = str(gdf.crs).rsplit(":",1)[-1] if gdf.crs else "None"
        print(f"  {layer}: {len(gdf)} feat, EPSG:{crs_s}, cols={gdf.columns.tolist()}")
    except Exception as e:
        print(f"  {layer}: ERROR - {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Direct Overpass queries for all data, geometric computation in-memory.
Output: tokyo_crossings.gpkg (5 layers).
"""

import json
import sys
import time
import warnings
from collections import defaultdict

import geopandas as gpd
import numpy as np
import pandas as pd
import requests
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiPoint, shape
from shapely.ops import unary_union

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

# The 23 special wards with their OSM relation IDs
WARD_RELATIONS = [
    (1759467, "千代田区", "Chiyoda"),
    (1760068, "中央区", "Chuo"),
    (1759475, "港区", "Minato"),
    (1759472, "新宿区", "Shinjuku"),
    (1759474, "文京区", "Bunkyo"),
    (1759471, "台東区", "Taito"),
    (1759470, "墨田区", "Sumida"),
    (1759469, "江東区", "Koto"),
    (1759477, "品川区", "Shinagawa"),
    (1759478, "目黒区", "Meguro"),
    (1759476, "大田区", "Ota"),
    (1759473, "世田谷区", "Setagaya"),
    (1759479, "渋谷区", "Shibuya"),
    (1759480, "中野区", "Nakano"),
    (1759481, "杉並区", "Suginami"),
    (1759484, "豊島区", "Toshima"),
    (1761848, "北区", "Kita"),
    (1759482, "荒川区", "Arakawa"),
    (1761847, "板橋区", "Itabashi"),
    (1761846, "練馬区", "Nerima"),
    (1761845, "足立区", "Adachi"),
    (1761844, "葛飾区", "Katsushika"),
    (1761843, "江戸川区", "Edogawa"),
]

OVERAPASS_URLS = [
    "https://overpass-api.de/api/interpreter",
    "https://overpass.kumi.systems/api/interpreter",
    "https://lz4.overpass-api.de/api/interpreter",
]


def overpass_query(query, timeout=300, max_retries=5):
    """Run an Overpass QL query with retries across multiple mirrors."""
    last_err = None
    for attempt in range(max_retries):
        for url in OVERAPASS_URLS:
            try:
                r = requests.get(url, params={"data": query}, timeout=timeout)
                r.raise_for_status()
                data = r.json()
                if "elements" in data:
                    return data
                if "remark" in data:
                    print(f"  Overpass remark: {data['remark']}", file=sys.stderr)
                return data
            except requests.exceptions.RequestException as e:
                last_err = e
                time.sleep(1)
        time.sleep(5 * (attempt + 1))
    raise RuntimeError(f"All Overpass attempts failed. Last error: {last_err}")


def build_polygon_from_outer_rings(outer_rings, inner_rings):
    """Build a Polygon from outer ring(s) with optional inner rings."""
    if not outer_rings:
        return None
    
    # Build main polygon from first outer ring
    main = Polygon(outer_rings[0])
    
    # Cut inner holes
    for ic in inner_rings:
        try:
            hole = Polygon(ic)
            if hole.is_valid:
                main = main.difference(hole)
        except Exception:
            pass
    
    # Add additional outer rings
    for oc in outer_rings[1:]:
        try:
            extra = Polygon(oc)
            if extra.is_valid:
                main = unary_union([main, extra])
        except Exception:
            pass
    
    if main.is_valid and not main.is_empty:
        return main
    return None


# ---------------------------------------------------------------------------
# Step 1: Fetch ward boundaries
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching ward boundaries ===", flush=True)

# Single query to get all 23 ward relations with full member geometries
rel_ids = ",".join(str(w[0]) for w in WARD_RELATIONS)
query_wards = f"""
[out:json][timeout:300];
(
  rel({rel_ids});
);
(._;>>;);
out body geom;
"""

print("  Fetching all 23 ward boundaries in one query...", flush=True)
try:
    ward_data = overpass_query(query_wards, timeout=300)
except RuntimeError as e:
    print(f"  Overpass failed: {e}")
    print("  Trying per-ward fallback...")
    ward_data = {"elements": []}
    for rid, jp_name, en_name in WARD_RELATIONS:
        q = f"[out:json][timeout:120];(rel({rid});(._;>>;);out body geom;);"
        for attempt in range(3):
            try:
                part = overpass_query(q, timeout=120)
                ward_data["elements"].extend(part["elements"])
                print(f"  Got {jp_name}")
                break
            except Exception as e:
                print(f"  {jp_name} attempt {attempt+1}: {e}")
                time.sleep(3)

print(f"  Got {len(ward_data.get('elements', []))} elements", flush=True)

# Collect all ways (by ID) that are part of relations
ways = {}
relations_found = []

for elem in ward_data.get("elements", []):
    if elem["type"] == "way":
        coords = [(p["lon"], p["lat"]) for p in elem.get("geometry", [])]
        ways[elem["id"]] = coords
    elif elem["type"] == "relation":
        relations_found.append(elem)

print(f"  Found {len(ways)} ways, {len(relations_found)} relations", flush=True)

# Build ward polygons from relations
wards_list = []
for rel in relations_found:
    rel_id = rel["id"]
    
    # Find the matching entry
    match = next((w for w in WARD_RELATIONS if w[0] == rel_id), None)
    if not match:
        continue
    
    _, jp_name, en_name = match
    
    # Collect outer and inner rings
    members = rel.get("members", [])
    outer_rings = []
    inner_rings = []
    
    for mem in members:
        way_id = mem.get("ref")
        role = mem.get("role", "")
        
        # Check if geometry is embedded on the member
        geom_coords = mem.get("geometry")
        
        coords = None
        if geom_coords and len(geom_coords) >= 3:
            coords = [(p["lon"], p["lat"]) for p in geom_coords]
        elif way_id in ways and len(ways[way_id]) >= 3:
            coords = ways[way_id]
        
        if coords is None:
            continue
        
        # Check if it's a closed ring
        if coords[0] != coords[-1]:
            continue
        
        if role == "outer":
            outer_rings.append(coords)
        elif role == "inner":
            inner_rings.append(coords)
    
    if outer_rings:
        poly = build_polygon_from_outer_rings(outer_rings, inner_rings)
        if poly is not None:
            wards_list.append({
                "ward_id": rel_id,
                "ward_name": jp_name,
                "ward_name_en": en_name,
                "geometry": poly
            })
            print(f"  ✓ {jp_name} ({en_name}) — osm_id={rel_id}")
        else:
            print(f"  ✗ {jp_name} — could not build polygon")
    else:
        print(f"  ✗ {jp_name} — no outer rings found")

print(f"\n  Built {len(wards_list)} / 23 ward polygons", flush=True)

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    sys.exit(1)

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)

# ---------------------------------------------------------------------------
# Step 2: Fetch drivable roads
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads ===", flush=True)

# Get bounding box
bbox = unary_union(wards_gdf.geometry.values).bounds
print(f"  Bounding box: {bbox}", flush=True)

# Query for all drivable road types (excluding service)
query_roads = f"""
[out:json][timeout:300];
(
  way["highway"~"^(motorway|motorway_link|trunk|trunk_link|primary|primary_link|secondary|secondary_link|tertiary|tertiary_link|unclassified|residential|living_street|road)$"]({bbox[1]},{bbox[0]},{bbox[3]},{bbox[2]});
);
out geom;
"""

print("  Fetching drivable roads...", flush=True)
try:
    road_data = overpass_query(query_roads, timeout=300)
    print(f"  Got {len(road_data.get('elements', []))} road elements", flush=True)
except RuntimeError as e:
    print(f"  Road fetch failed: {e}", flush=True)
    road_data = {"elements": []}

# Process road geometries
road_geoms = []
for elem in road_data.get("elements", []):
    if elem["type"] != "way":
        continue
    coords = [(p["lon"], p["lat"]) for p in elem.get("geometry", [])]
    if len(coords) < 2:
        continue
    road_geoms.append(LineString(coords))

print(f"  Found {len(road_geoms)} road segments", flush=True)

if len(road_geoms) == 0:
    print("CRITICAL: No roads fetched. Exiting.")
    sys.exit(1)

roads_gdf = gpd.GeoDataFrame(geometry=road_geoms, crs=WGS84)

# ---------------------------------------------------------------------------
# Step 3: Crossing points
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing crossing points ===", flush=True)

print("  Projecting to metric CRS...", flush=True)
wards_metric = wards_gdf.to_crs(METRIC_CRS)
roads_metric = roads_gdf.to_crs(METRIC_CRS)

# Clip roads to ward boundaries
union_metric = unary_union(wards_metric.geometry.values)
print("  Clipping roads to ward union...", flush=True)
roads_metric["geometry"] = roads_metric.geometry.intersection(union_metric)
roads_metric = roads_metric[~roads_metric.geometry.is_empty & 
                             (roads_metric.geometry.geom_type == "LineString")].copy()
print(f"  Roads within wards: {len(roads_metric)}", flush=True)

# Build boundary segments
print("  Building boundary segments...", flush=True)
boundary_parts = []
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_parts.append((ward_id, line))
    elif boundary.geom_type == "LineString":
        boundary_parts.append((ward_id, boundary))

print(f"  {len(boundary_parts)} boundary segments", flush=True)

# Spatial index
roads_sindex = roads_metric.sindex

# Find crossings
crossings = []
for ward_id, boundary_line in boundary_parts:
    if boundary_line.length < 1.0:
        continue
    
    bounds = boundary_line.bounds
    candidates = list(roads_sindex.intersection(bounds))
    
    for ridx in candidates:
        road_geom = roads_metric.iloc[ridx].geometry
        if not boundary_line.intersects(road_geom):
            continue
        
        inter = boundary_line.intersection(road_geom)
        if inter.is_empty:
            continue
        
        if inter.geom_type == "Point":
            crossings.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": ward_id, "geometry": g})

print(f"  Raw crossings: {len(crossings)}", flush=True)

if crossings:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["ward_id", "x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    name_map = {}
    for _, wrow in wards_metric.iterrows():
        name_map[wrow["ward_id"]] = (wrow["ward_name"], wrow["ward_name_en"])
    
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: name_map.get(x, ("", ""))[1])
    
    print(f"  Unique crossings: {len(crossing_gdf)}", flush=True)
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings found!", flush=True)

# ---------------------------------------------------------------------------
# Step 4: 50m buffers
# ---------------------------------------------------------------------------
print("\n=== Step 4: 50m buffers ===", flush=True)

if len(crossing_gdf) > 0:
    buffer_gdf = crossing_gdf.copy()
    buffer_gdf["geometry"] = buffer_gdf.geometry.buffer(50.0)
    buffer_gdf = buffer_gdf.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers", flush=True)
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 5: Clip to ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to ward ===", flush=True)

if len(buffer_gdf) > 0:
    ward_geom_map = {}
    for _, wrow in wards_metric.iterrows():
        ward_geom_map[wrow["ward_id"]] = wrow.geometry
    
    clipped_geoms = []
    for _, brow in buffer_gdf.iterrows():
        wg = ward_geom_map.get(brow["ward_id"])
        if wg is not None:
            c = brow.geometry.intersection(wg)
            clipped_geoms.append(c if not c.is_empty else None)
        else:
            clipped_geoms.append(None)
    
    clipped_gdf = buffer_gdf.copy()
    clipped_gdf["geometry"] = clipped_geoms
    clipped_gdf = clipped_gdf[clipped_gdf.geometry.notna() & ~clipped_gdf.geometry.is_empty].copy()
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    clipped_gdf["geometry"] = clipped_gdf.geometry.buffer(0)
    clipped_gdf = clipped_gdf.set_geometry("geometry")
    print(f"  Created {len(clipped_gdf)} clipped buffers", flush=True)
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 6: Density
# ---------------------------------------------------------------------------
print("\n=== Step 6: Crossing density ===", flush=True)

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

ward_areas = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
ward_areas["ward_area_km2"] = ward_areas.geometry.area / 1_000_000

density_df = ward_areas.merge(counts, on="ward_id", how="left")
density_df["crossing_count"] = density_df["crossing_count"].fillna(0).astype(int)
density_df["crossings_per_km2"] = density_df["crossing_count"] / density_df["ward_area_km2"]

density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(density_df[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===", flush=True)

wards_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wards_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wards_out)} features")

if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print("  ✓ crossing_points: 0 features (empty)")

if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_gdf)} features")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, 
                     geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    ).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)} features")

print("\n=== DONE ===", flush=True)

# Verify
print("\n=== Verification ===")
for layer in ["wards_jgd","crossing_points","crossing_buffers_50m","buffer_ward_intersection","ward_crossing_density_wgs84"]:
    try:
        gdf = gpd.read_file(OUTPUT, layer=layer)
        crs_s = str(gdf.crs).rsplit(":",1)[-1] if gdf.crs else "None"
        print(f"  {layer}: {len(gdf)} feat, EPSG:{crs_s}, cols={gdf.columns.tolist()}")
    except Exception as e:
        print(f"  {layer}: ERROR - {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Direct Overpass queries with proper headers.
Output: tokyo_crossings.gpkg (5 layers).
"""

import json
import sys
import time
import warnings
from collections import defaultdict

import geopandas as gpd
import numpy as np
import pandas as pd
import requests
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiPoint, shape
from shapely.ops import unary_union

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

# The 23 special wards with their OSM relation IDs
WARD_RELATIONS = [
    (1759467, "千代田区", "Chiyoda"),
    (1760068, "中央区", "Chuo"),
    (1759475, "港区", "Minato"),
    (1759472, "新宿区", "Shinjuku"),
    (1759474, "文京区", "Bunkyo"),
    (1759471, "台東区", "Taito"),
    (1759470, "墨田区", "Sumida"),
    (1759469, "江東区", "Koto"),
    (1759477, "品川区", "Shinagawa"),
    (1759478, "目黒区", "Meguro"),
    (1759476, "大田区", "Ota"),
    (1759473, "世田谷区", "Setagaya"),
    (1759479, "渋谷区", "Shibuya"),
    (1759480, "中野区", "Nakano"),
    (1759481, "杉並区", "Suginami"),
    (1759484, "豊島区", "Toshima"),
    (1761848, "北区", "Kita"),
    (1759482, "荒川区", "Arakawa"),
    (1761847, "板橋区", "Itabashi"),
    (1761846, "練馬区", "Nerima"),
    (1761845, "足立区", "Adachi"),
    (1761844, "葛飾区", "Katsushika"),
    (1761843, "江戸川区", "Edogawa"),
]

OVERPASS_URL = "https://overpass-api.de/api/interpreter"
HEADERS = {
    "User-Agent": "OSMnx/2.1.0 (https://github.com/gboeing/osmnx)",
    "Accept": "*/*",
}


def overpass_query(query_text, timeout=300, max_retries=5):
    """Run an Overpass QL query with retries and backoff."""
    last_err = None
    for attempt in range(max_retries):
        try:
            params = {"data": query_text}
            r = requests.get(OVERPASS_URL, params=params, headers=HEADERS, timeout=timeout)
            r.raise_for_status()
            data = r.json()
            if "elements" in data:
                return data
            if "remark" in data:
                print(f"  Remark: {data['remark']}", file=sys.stderr)
            return data
        except (requests.RequestException, json.JSONDecodeError) as e:
            last_err = e
            print(f"  Attempt {attempt + 1}/{max_retries} failed: {type(e).__name__}: {e}", file=sys.stderr)
            time.sleep(3 * (attempt + 1))
    raise RuntimeError(f"All Overpass attempts failed. Last: {last_err}")


def build_polygon_from_ring_sets(outer_rings, inner_rings):
    """Build a valid Polygon/MultiPolygon from ring coordinate lists."""
    if not outer_rings:
        return None
    
    main = Polygon(outer_rings[0])
    for ic in inner_rings:
        try:
            hole = Polygon(ic)
            if hole.is_valid:
                main = main.difference(hole)
        except Exception:
            pass
    for oc in outer_rings[1:]:
        try:
            extra = Polygon(oc)
            if extra.is_valid:
                main = unary_union([main, extra])
        except Exception:
            pass
    
    if main.is_valid and not main.is_empty:
        return main
    return None


# ---------------------------------------------------------------------------
# Step 1: Fetch ward boundaries
# ---------------------------------------------------------------------------
print("=== Step 1: Fetching ward boundaries ===", flush=True)

rel_ids = ",".join(str(w[0]) for w in WARD_RELATIONS)
q_wards = f"""
[out:json][timeout:300];
(
  rel({rel_ids});
);
(._;>>;);
out body geom;
"""

print("  Fetching ward boundaries...", flush=True)
try:
    ward_data = overpass_query(q_wards, timeout=300)
except RuntimeError as e:
    print(f"  Batch query failed: {e}", flush=True)
    print("  Trying per-ward...", flush=True)
    ward_data = {"elements": []}
    for rid, jp_name, en_name in WARD_RELATIONS:
        q = f"[out:json][timeout:180];(rel({rid});(._;>>;);out body geom;);"
        try:
            part = overpass_query(q, timeout=180)
            ward_data["elements"].extend(part["elements"])
            print(f"    ✓ {jp_name}", flush=True)
        except Exception as e:
            print(f"    ✗ {jp_name}: {e}", flush=True)
        time.sleep(1)

print(f"  Got {len(ward_data.get('elements', []))} elements", flush=True)

# Separate ways and relations
ways = {}
relations = []
for elem in ward_data.get("elements", []):
    if elem["type"] == "way":
        coords = [(p["lon"], p["lat"]) for p in elem.get("geometry", [])]
        ways[elem["id"]] = coords
    elif elem["type"] == "relation":
        relations.append(elem)

print(f"  {len(ways)} ways, {len(relations)} relations", flush=True)

# Build ward polygons
wards_list = []
for rel in relations:
    rel_id = rel["id"]
    match = next((w for w in WARD_RELATIONS if w[0] == rel_id), None)
    if not match:
        continue
    _, jp_name, en_name = match
    
    members = rel.get("members", [])
    outer_rings = []
    inner_rings = []
    
    for mem in members:
        way_id = mem.get("ref")
        role = mem.get("role", "")
        
        # Get coordinates from member geometry or from ways dict
        geom_coords = mem.get("geometry")
        coords = None
        if geom_coords and len(geom_coords) >= 3:
            coords = [(p["lon"], p["lat"]) for p in geom_coords]
        elif way_id in ways and len(ways[way_id]) >= 3:
            coords = ways[way_id]
        
        if not coords:
            continue
        if coords[0] != coords[-1]:
            continue  # Not a closed ring
        
        if role == "outer":
            outer_rings.append(coords)
        elif role == "inner":
            inner_rings.append(coords)
    
    if outer_rings:
        poly = build_polygon_from_ring_sets(outer_rings, inner_rings)
        if poly is not None:
            wards_list.append({
                "ward_id": rel_id,
                "ward_name": jp_name,
                "ward_name_en": en_name,
                "geometry": poly
            })
            print(f"  ✓ {jp_name} — osm_id={rel_id}", flush=True)
        else:
            print(f"  ✗ {jp_name} — invalid geometry", flush=True)
    else:
        print(f"  ✗ {jp_name} — no outer rings", flush=True)

print(f"\n  Built {len(wards_list)} / 23 ward polygons", flush=True)

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    sys.exit(1)

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
for _, row in wards_gdf.iterrows():
    print(f"    {row['ward_name_en']:>12s} (ID={row['ward_id']})")

# ---------------------------------------------------------------------------
# Step 2: Fetch drivable roads
# ---------------------------------------------------------------------------
print("\n=== Step 2: Fetching drivable roads ===", flush=True)

bbox = unary_union(wards_gdf.geometry.values).bounds
print(f"  BBox: {bbox}", flush=True)

q_roads = f"""
[out:json][timeout:300];
(
  way["highway"~"^(motorway|motorway_link|trunk|trunk_link|primary|primary_link|secondary|secondary_link|tertiary|tertiary_link|unclassified|residential|living_street|road)$"]({bbox[1]},{bbox[0]},{bbox[3]},{bbox[2]});
);
out geom;
"""

print("  Fetching drivable roads...", flush=True)
try:
    road_data = overpass_query(q_roads, timeout=300)
    print(f"  Got {len(road_data.get('elements', []))} road elements", flush=True)
except RuntimeError as e:
    print(f"  Road fetch failed: {e}", flush=True)
    road_data = {"elements": []}

road_geoms = []
for elem in road_data.get("elements", []):
    if elem["type"] != "way":
        continue
    coords = [(p["lon"], p["lat"]) for p in elem.get("geometry", [])]
    if len(coords) < 2:
        continue
    road_geoms.append(LineString(coords))

print(f"  Found {len(road_geoms)} road segments", flush=True)
if len(road_geoms) == 0:
    print("CRITICAL: No roads fetched. Exiting.")
    sys.exit(1)

roads_gdf = gpd.GeoDataFrame(geometry=road_geoms, crs=WGS84)

# ---------------------------------------------------------------------------
# Step 3: Crossing points
# ---------------------------------------------------------------------------
print("\n=== Step 3: Computing crossing points ===", flush=True)

print("  Projecting to metric CRS...", flush=True)
wards_metric = wards_gdf.to_crs(METRIC_CRS)
roads_metric = roads_gdf.to_crs(METRIC_CRS)

# Clip roads to ward union
union_metric = unary_union(wards_metric.geometry.values)
print("  Clipping roads to ward union...", flush=True)
roads_metric["geometry"] = roads_metric.geometry.intersection(union_metric)
roads_metric = roads_metric[~roads_metric.geometry.is_empty &
                             (roads_metric.geometry.geom_type == "LineString")].copy()
print(f"  Roads in wards: {len(roads_metric)}", flush=True)

# Boundary segments
print("  Building boundary segments...", flush=True)
boundary_parts = []
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_parts.append((ward_id, line))
    elif boundary.geom_type == "LineString":
        boundary_parts.append((ward_id, boundary))

print(f"  {len(boundary_parts)} segments", flush=True)

# Spatial index
roads_sindex = roads_metric.sindex

# Find crossings
crossings = []
for ward_id, boundary_line in boundary_parts:
    if boundary_line.length < 1.0:
        continue
    
    candidates = list(roads_sindex.intersection(boundary_line.bounds))
    for ridx in candidates:
        rg = roads_metric.iloc[ridx].geometry
        if not boundary_line.intersects(rg):
            continue
        inter = boundary_line.intersection(rg)
        if inter.is_empty:
            continue
        if inter.geom_type == "Point":
            crossings.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": ward_id, "geometry": g})

print(f"  Raw crossings: {len(crossings)}", flush=True)

if crossings:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["ward_id", "x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    nm = {}
    for _, wr in wards_metric.iterrows():
        nm[wr["ward_id"]] = (wr["ward_name"], wr["ward_name_en"])
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: nm.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: nm.get(x, ("", ""))[1])
    
    print(f"  Unique crossings: {len(crossing_gdf)}", flush=True)
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings!", flush=True)

# ---------------------------------------------------------------------------
# Step 4: 50m buffers
# ---------------------------------------------------------------------------
print("\n=== Step 4: 50m buffers ===", flush=True)

if len(crossing_gdf) > 0:
    bg = crossing_gdf.copy()
    bg["geometry"] = bg.geometry.buffer(50.0)
    buffer_gdf = bg.set_geometry("geometry")
    print(f"  Created {len(buffer_gdf)} buffers", flush=True)
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 5: Clip buffers to ward
# ---------------------------------------------------------------------------
print("\n=== Step 5: Clipping buffers to ward ===", flush=True)

if len(buffer_gdf) > 0:
    wgm = {wr["ward_id"]: wr.geometry for _, wr in wards_metric.iterrows()}
    clip_geoms = []
    for _, br in buffer_gdf.iterrows():
        wg = wgm.get(br["ward_id"])
        if wg is not None:
            c = br.geometry.intersection(wg)
            clip_geoms.append(c if not c.is_empty else None)
        else:
            clip_geoms.append(None)
    
    cg = buffer_gdf.copy()
    cg["geometry"] = clip_geoms
    cg = cg[cg.geometry.notna() & ~cg.geometry.is_empty].copy()
    cg = cg.set_geometry("geometry")
    cg["geometry"] = cg.geometry.buffer(0)
    clipped_gdf = cg.set_geometry("geometry")
    print(f"  Created {len(clipped_gdf)} clipped buffers", flush=True)
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ---------------------------------------------------------------------------
# Step 6: Density
# ---------------------------------------------------------------------------
print("\n=== Step 6: Crossing density ===", flush=True)

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

wa = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wa["ward_area_km2"] = wa.geometry.area / 1_000_000

dd = wa.merge(counts, on="ward_id", how="left")
dd["crossing_count"] = dd["crossing_count"].fillna(0).astype(int)
dd["crossings_per_km2"] = dd["crossing_count"] / dd["ward_area_km2"]

density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(dd[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ---------------------------------------------------------------------------
# Save
# ---------------------------------------------------------------------------
print("\n=== Saving to GeoPackage ===", flush=True)

wr_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wr_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wr_out)}")

if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_gdf)}")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print("  ✓ crossing_points: 0 (empty)")

if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_gdf)}")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_gdf)}")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)}")

print("\n=== DONE ===", flush=True)

# Verify
print("\n=== Verification ===")
for layer in ["wards_jgd","crossing_points","crossing_buffers_50m","buffer_ward_intersection","ward_crossing_density_wgs84"]:
    try:
        g = gpd.read_file(OUTPUT, layer=layer)
        cs = str(g.crs).rsplit(":",1)[-1] if g.crs else "None"
        print(f"  {layer}: {len(g)} feat, EPSG:{cs}, cols={g.columns.tolist()}")
    except Exception as e:
        print(f"  {layer}: {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Uses OSMnx for per-ward boundary and road fetching, with robust name matching.
Output: tokyo_crossings.gpkg (5 layers).
"""

import sys
import time
import warnings
from shapely.geometry import Point, LineString, MultiPoint
from shapely.ops import unary_union
from shapely import wkt

import geopandas as gpd
import numpy as np
import pandas as pd
import osmnx as ox

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

ox.settings.timeout = 300
ox.settings.max_query_area_size = 50 * 1000 * 1000
ox.settings.bulk_download = True
ox.settings.requests_kwargs = {"timeout": 300}

WARD_MAP = {
    "千代田区": "Chiyoda",
    "中央区": "Chuo",
    "港区": "Minato",
    "新宿区": "Shinjuku",
    "文京区": "Bunkyo",
    "台東区": "Taito",
    "墨田区": "Sumida",
    "江東区": "Koto",
    "品川区": "Shinagawa",
    "目黒区": "Meguro",
    "大田区": "Ota",
    "世田谷区": "Setagaya",
    "渋谷区": "Shibuya",
    "中野区": "Nakano",
    "杉並区": "Suginami",
    "豊島区": "Toshima",
    "北区": "Kita",
    "荒川区": "Arakawa",
    "板橋区": "Itabashi",
    "練馬区": "Nerima",
    "足立区": "Adachi",
    "葛飾区": "Katsushika",
    "江戸川区": "Edogawa",
}

OVERPASS_MIRRORS = [
    "https://overpass-api.de/api",
    "https://overpass.kumi.systems/api",
    "https://lz4.overpass-api.de/api",
]


def fetch_ward(jp_name, en_name):
    """Fetch one ward boundary. Tries multiple Overpass mirrors."""
    place = f"{jp_name}, Tokyo, Japan"
    
    for mirror in OVERPASS_MIRRORS:
        ox.settings.overpass_url = mirror
        for attempt in range(2):
            try:
                gdf = ox.features_from_place(place, tags={"boundary": "administrative"})
                
                # Find a polygon whose name matches this ward exactly
                for idx, row in gdf.iterrows():
                    geom = row.geometry
                    if geom is None or geom.is_empty:
                        continue
                    if geom.geom_type not in ("Polygon", "MultiPolygon"):
                        continue
                    
                    name = str(row.get("name", ""))
                    name_ja = str(row.get("name:ja", ""))
                    
                    # Name must match the Japanese ward name
                    if name != jp_name and name_ja != jp_name:
                        continue
                    
                    osm_id = idx[1] if isinstance(idx, tuple) else idx
                    return (int(osm_id), geom)
                
                return None  # No matching ward found
                
            except Exception as e:
                err_str = str(e)
                if "refused" in err_str or "timed out" in err_str or "timeout" in err_str:
                    break  # Try next mirror
                time.sleep(2)
    return None


def fetch_ward_roads(ward_polygon):
    """Fetch drivable (non-service) roads within a ward polygon."""
    for mirror in OVERPASS_MIRRORS:
        ox.settings.overpass_url = mirror
        for attempt in range(2):
            try:
                G = ox.graph_from_polygon(ward_polygon, network_type="drive",
                                           simplify=True, retain_all=True)
                edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
                if edges is None or len(edges) == 0:
                    return []
                
                if "geometry" not in edges.columns:
                    return []
                
                edges = edges[edges.geometry.notna() & ~edges.geometry.is_empty].copy()
                
                # Filter service roads
                if "highway" in edges.columns:
                    def ok(hw):
                        if isinstance(hw, list):
                            return not any("service" in str(h).lower() for h in hw)
                        return "service" not in str(hw).lower()
                    mask = edges["highway"].apply(ok)
                    edges = edges[mask].copy()
                
                # Normalize to LineString
                def to_ls(g):
                    if hasattr(g, 'geom_type') and g.geom_type == "MultiLineString":
                        parts = list(g.geoms)
                        return max(parts, key=lambda x: x.length) if parts else g
                    return g
                
                edges["geometry"] = edges["geometry"].apply(to_ls)
                edges = edges[edges.geometry.geom_type == "LineString"].copy()
                
                return list(edges.geometry.values)
                
            except Exception as e:
                err_str = str(e)
                if "refused" in err_str or "timed out" in err_str:
                    break
                time.sleep(2)
    return []


# ===========================================================================
# Step 1: Fetch ward boundaries
# ===========================================================================
print("=== Step 1: Fetching ward boundaries ===", flush=True)

wards_list = []
for jp_name, en_name in WARD_MAP.items():
    print(f"  {jp_name} ({en_name})...", end=" ", flush=True)
    result = fetch_ward(jp_name, en_name)
    if result is not None:
        osm_id, geom = result
        wards_list.append({
            "ward_id": osm_id,
            "ward_name": jp_name,
            "ward_name_en": en_name,
            "geometry": geom
        })
        print(f"✓ osm_id={osm_id}")
    else:
        print("✗ FAILED")

print(f"\n  Got {len(wards_list)} / 23", flush=True)

if len(wards_list) == 0:
    print("CRITICAL: Could not fetch any ward boundaries. Exiting.")
    sys.exit(1)

# Check for duplicates
ward_ids = [w["ward_id"] for w in wards_list]
if len(set(ward_ids)) < len(ward_ids):
    print("  WARNING: Duplicate ward IDs. Some wards may share boundaries.", flush=True)
    from collections import Counter
    for wid, cnt in Counter(ward_ids).items():
        if cnt > 1:
            names = [w["ward_name"] for w in wards_list if w["ward_id"] == wid]
            print(f"    ID {wid}: {names}")

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
for _, row in wards_gdf.iterrows():
    print(f"    {row['ward_name_en']:>12s} (ID={row['ward_id']})")

# ===========================================================================
# Step 2: Fetch drivable roads (per-ward, deduplicate)
# ===========================================================================
print("\n=== Step 2: Fetching drivable roads ===", flush=True)

# Use a slightly buffered union to get roads that cross ward boundaries
union_buffered = unary_union(wards_gdf.geometry.values).buffer(0.002)  # ~200m buffer

all_road_geoms = []
for idx, row in wards_gdf.iterrows():
    jp_name, en_name = row["ward_name"], row["ward_name_en"]
    print(f"  {jp_name} ({en_name})...", end=" ", flush=True)
    
    road_geoms = fetch_ward_roads(row.geometry)
    if road_geoms:
        print(f"{len(road_geoms)} roads", flush=True)
        all_road_geoms.extend(road_geoms)
    else:
        print("no roads", flush=True)

print(f"\n  Total raw: {len(all_road_geoms)}", flush=True)

if len(all_road_geoms) == 0:
    print("CRITICAL: No roads fetched. Exiting.")
    sys.exit(1)

# Deduplicate
seen = set()
unique_roads = []
for g in all_road_geoms:
    key = wkt.dumps(g, rounding_precision=6)
    if key not in seen:
        seen.add(key)
        unique_roads.append(g)

print(f"  Unique: {len(unique_roads)}", flush=True)

# Also fetch roads in the buffered border region to catch boundary-crossing roads
print("  Fetching roads in border region...", flush=True)
border_roads = fetch_ward_roads(union_buffered)
if border_roads:
    for g in border_roads:
        key = wkt.dumps(g, rounding_precision=6)
        if key not in seen:
            seen.add(key)
            unique_roads.append(g)
    print(f"  Added border roads, total unique: {len(unique_roads)}", flush=True)

roads_gdf = gpd.GeoDataFrame(geometry=unique_roads, crs=WGS84)

# ===========================================================================
# Step 3: Find crossing points
# ===========================================================================
print("\n=== Step 3: Computing crossing points ===", flush=True)

print("  Projecting to metric CRS...", flush=True)
wards_metric = wards_gdf.to_crs(METRIC_CRS)
roads_metric = roads_gdf.to_crs(METRIC_CRS)

# Build boundary segments (each ward's boundary line)
print("  Building boundary segments...", flush=True)
boundary_parts = []
for idx, row in wards_metric.iterrows():
    ward_id = row["ward_id"]
    boundary = row.geometry.boundary
    if boundary.geom_type == "MultiLineString":
        for line in boundary.geoms:
            boundary_parts.append((ward_id, line))
    elif boundary.geom_type == "LineString":
        boundary_parts.append((ward_id, boundary))

print(f"  {len(boundary_parts)} segments", flush=True)

# Spatial index on roads
roads_sindex = roads_metric.sindex

# Find intersections
crossings = []
for ward_id, boundary_line in boundary_parts:
    if boundary_line.length < 1.0:
        continue
    
    candidates = list(roads_sindex.intersection(boundary_line.bounds))
    for ridx in candidates:
        rg = roads_metric.iloc[ridx].geometry
        if not boundary_line.intersects(rg):
            continue
        inter = boundary_line.intersection(rg)
        if inter.is_empty:
            continue
        if inter.geom_type == "Point":
            crossings.append({"ward_id": ward_id, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": ward_id, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": ward_id, "geometry": g})

print(f"  Raw crossings: {len(crossings)}", flush=True)

if crossings:
    crossing_gdf = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    crossing_gdf["x"] = crossing_gdf.geometry.x.round(2)
    crossing_gdf["y"] = crossing_gdf.geometry.y.round(2)
    crossing_gdf = crossing_gdf.drop_duplicates(subset=["ward_id", "x", "y"])
    crossing_gdf = crossing_gdf.drop(columns=["x", "y"])
    
    nm = {}
    for _, wr in wards_metric.iterrows():
        nm[wr["ward_id"]] = (wr["ward_name"], wr["ward_name_en"])
    crossing_gdf["ward_name"] = crossing_gdf["ward_id"].map(lambda x: nm.get(x, ("", ""))[0])
    crossing_gdf["ward_name_en"] = crossing_gdf["ward_id"].map(lambda x: nm.get(x, ("", ""))[1])
    
    print(f"  Unique: {len(crossing_gdf)}", flush=True)
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings!", flush=True)

# ===========================================================================
# Step 4: 50m buffers
# ===========================================================================
print("\n=== Step 4: 50m buffers ===", flush=True)

if len(crossing_gdf) > 0:
    bg = crossing_gdf.copy()
    bg["geometry"] = bg.geometry.buffer(50.0)
    buffer_gdf = bg.set_geometry("geometry")
    print(f"  {len(buffer_gdf)} buffers", flush=True)
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ===========================================================================
# Step 5: Clip buffers to ward
# ===========================================================================
print("\n=== Step 5: Clipping to ward ===", flush=True)

if len(buffer_gdf) > 0:
    wgm = {wr["ward_id"]: wr.geometry for _, wr in wards_metric.iterrows()}
    clip_geoms = []
    for _, br in buffer_gdf.iterrows():
        wg = wgm.get(br["ward_id"])
        clip_geoms.append(br.geometry.intersection(wg) if wg is not None else None)
    
    cg = buffer_gdf.copy()
    cg["geometry"] = [None if (c is None or c.is_empty) else c for c in clip_geoms]
    cg = cg[cg.geometry.notna() & ~cg.geometry.is_empty].copy()
    cg = cg.set_geometry("geometry")
    cg["geometry"] = cg.geometry.buffer(0)
    clipped_gdf = cg.set_geometry("geometry")
    print(f"  {len(clipped_gdf)} clipped", flush=True)
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"), "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ===========================================================================
# Step 6: Density
# ===========================================================================
print("\n=== Step 6: Crossing density ===", flush=True)

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

wa = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wa["ward_area_km2"] = wa.geometry.area / 1_000_000

dd = wa.merge(counts, on="ward_id", how="left")
dd["crossing_count"] = dd["crossing_count"].fillna(0).astype(int)
dd["crossings_per_km2"] = dd["crossing_count"] / dd["ward_area_km2"]

density_gdf = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
density_gdf = density_gdf.merge(dd[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
density_gdf = density_gdf.to_crs(WGS84)

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ===========================================================================
# Save
# ===========================================================================
print("\n=== Saving ===", flush=True)

wr_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wr_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print(f"  ✓ wards_jgd: {len(wr_out)}")

if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points: {len(crossing_gdf)}")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_points", driver="GPKG")

if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m: {len(buffer_gdf)}")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection: {len(clipped_gdf)}")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84: {len(density_gdf)}")

print("\n=== DONE ===", flush=True)

# Verify
print("\n=== Verification ===")
for layer in ["wards_jgd","crossing_points","crossing_buffers_50m","buffer_ward_intersection","ward_crossing_density_wgs84"]:
    try:
        g = gpd.read_file(OUTPUT, layer=layer)
        cs = str(g.crs).rsplit(":",1)[-1] if g.crs else "None"
        print(f"  {layer}: {len(g)} feat, EPSG:{cs}, cols={g.columns.tolist()}")
    except Exception as e:
        print(f"  {layer}: {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Uses OSMnx per-ward. Skips wards that can't be fetched but tries all.
Output: tokyo_crossings.gpkg (5 layers).
"""

import sys
import time
import warnings
from collections import Counter
from shapely.geometry import Point, LineString, MultiPoint
from shapely.ops import unary_union
from shapely import wkt

import geopandas as gpd
import numpy as np
import pandas as pd
import osmnx as ox

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

WARD_MAP = {
    "千代田区": "Chiyoda", "中央区": "Chuo", "港区": "Minato",
    "新宿区": "Shinjuku", "文京区": "Bunkyo", "台東区": "Taito",
    "墨田区": "Sumida", "江東区": "Koto", "品川区": "Shinagawa",
    "目黒区": "Meguro", "大田区": "Ota", "世田谷区": "Setagaya",
    "渋谷区": "Shibuya", "中野区": "Nakano", "杉並区": "Suginami",
    "豊島区": "Toshima", "北区": "Kita", "荒川区": "Arakawa",
    "板橋区": "Itabashi", "練馬区": "Nerima", "足立区": "Adachi",
    "葛飾区": "Katsushika", "江戸川区": "Edogawa",
}

# OSM overpass mirrors
MIRRORS = [
    "https://overpass-api.de/api",
    "https://overpass.kumi.systems/api",
    "https://lz4.overpass-api.de/api",
]

ox.settings.timeout = 300
ox.settings.max_query_area_size = 50 * 1000 * 1000
ox.settings.bulk_download = True
ox.settings.requests_kwargs = {"timeout": 300}


def fetch_one_ward(jp_name, en_name):
    """Fetch boundary polygon for one ward. Returns (osm_id, geom) or None."""
    place = f"{jp_name}, Tokyo, Japan"
    for mirror in MIRRORS:
        ox.settings.overpass_url = mirror
        for attempt in range(2):
            try:
                gdf = ox.features_from_place(place, tags={"boundary": "administrative"})
                for idx, row in gdf.iterrows():
                    geom = row.geometry
                    if geom is None or geom.is_empty:
                        continue
                    if geom.geom_type not in ("Polygon", "MultiPolygon"):
                        continue
                    name = str(row.get("name", ""))
                    if name != jp_name:
                        continue
                    if str(row.get("admin_level", "")).strip() != "7":
                        continue
                    osm_id = idx[1] if isinstance(idx, tuple) else idx
                    return (int(osm_id), geom)
                return None  # No valid polygon found
            except Exception as e:
                err = str(e)
                if "refused" in err or "timed out" in err or "timeout" in err:
                    break
                time.sleep(2)
    return None


def fetch_ward_roads(polygon):
    """Fetch drivable (non-service) roads inside a polygon."""
    for mirror in MIRRORS:
        ox.settings.overpass_url = mirror
        for attempt in range(2):
            try:
                G = ox.graph_from_polygon(polygon, network_type="drive",
                                           simplify=True, retain_all=True)
                edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
                if edges is None or len(edges) == 0:
                    return []
                if "geometry" not in edges.columns:
                    return []
                edges = edges[edges.geometry.notna() & ~edges.geometry.is_empty].copy()
                
                # Filter service roads
                if "highway" in edges.columns:
                    def no_serv(hw):
                        if isinstance(hw, list):
                            return not any("service" in str(h).lower() for h in hw)
                        return "service" not in str(hw).lower()
                    edges = edges[edges["highway"].apply(no_serv)].copy()
                
                # Normalize to LineString
                def to_ls(g):
                    if hasattr(g, 'geom_type') and g.geom_type == "MultiLineString":
                        parts = list(g.geoms)
                        return max(parts, key=lambda x: x.length) if parts else g
                    return g
                edges["geometry"] = edges["geometry"].apply(to_ls)
                edges = edges[edges.geometry.geom_type == "LineString"].copy()
                return list(edges.geometry.values)
            except Exception as e:
                err = str(e)
                if "refused" in err or "timed out" in err:
                    break
                time.sleep(2)
    return []


# ===============================================================
# STEP 1: Fetch ward boundaries
# ===============================================================
print("=== Step 1: Fetching ward boundaries ===", flush=True)

wards_list = []
failed = []

for jp_name, en_name in WARD_MAP.items():
    print(f"  {jp_name} ({en_name})...", end=" ", flush=True)
    result = fetch_one_ward(jp_name, en_name)
    if result is not None:
        osm_id, geom = result
        wards_list.append({
            "ward_id": osm_id,
            "ward_name": jp_name,
            "ward_name_en": en_name,
            "geometry": geom
        })
        print(f"✓ osm_id={osm_id}")
    else:
        failed.append((jp_name, en_name))
        print("✗")

print(f"\n  Got {len(wards_list)}/23", flush=True)
if failed:
    print(f"  Failed: {[e for _, e in failed]}", flush=True)

if len(wards_list) == 0:
    print("CRITICAL: No wards. Exiting.")
    sys.exit(1)

# Check duplicates
ids = [w["ward_id"] for w in wards_list]
if len(set(ids)) < len(ids):
    print("  WARNING: Duplicate IDs detected", flush=True)
    for wid, cnt in Counter(ids).items():
        if cnt > 1:
            names = [w["ward_name"] for w in wards_list if w["ward_id"] == wid]
            print(f"    ID {wid}: {names}")

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
for _, r in wards_gdf.iterrows():
    print(f"    {r['ward_name_en']:>12s} (ID={r['ward_id']})")

# ===============================================================
# STEP 2: Fetch roads per ward, then merge
# ===============================================================
print("\n=== Step 2: Fetching drivable roads ===", flush=True)

all_geoms = []
for idx, row in wards_gdf.iterrows():
    jp, en = row["ward_name"], row["ward_name_en"]
    print(f"  {jp} ({en})...", end=" ", flush=True)
    geoms = fetch_ward_roads(row.geometry)
    if geoms:
        print(f"{len(geoms)} roads")
        all_geoms.extend(geoms)
    else:
        print("0 roads")

print(f"\n  Total: {len(all_geoms)}", flush=True)

if len(all_geoms) == 0:
    print("CRITICAL: No roads. Exiting.")
    sys.exit(1)

# Dedup
seen = set()
unique = []
for g in all_geoms:
    k = wkt.dumps(g, rounding_precision=6)
    if k not in seen:
        seen.add(k)
        unique.append(g)
print(f"  Unique: {len(unique)}", flush=True)

# Also get roads in the border region (to catch boundary-crossing roads)
union_geom = unary_union(wards_gdf.geometry.values)
border = union_geom.buffer(0.003)
print("  Fetching border roads...", end=" ", flush=True)
border_gs = fetch_ward_roads(border)
if border_gs:
    added = 0
    for g in border_gs:
        k = wkt.dumps(g, rounding_precision=6)
        if k not in seen:
            seen.add(k)
            unique.append(g)
            added += 1
    print(f"+{added} new")
else:
    print("0")

roads_gdf = gpd.GeoDataFrame(geometry=unique, crs=WGS84)
print(f"  Total road segments: {len(roads_gdf)}", flush=True)

# ===============================================================
# STEP 3: Find crossing points
# ===============================================================
print("\n=== Step 3: Computing crossing points ===", flush=True)

print("  Projecting to metric CRS...", flush=True)
wards_metric = wards_gdf.to_crs(METRIC_CRS)
roads_metric = roads_gdf.to_crs(METRIC_CRS)

# Clip roads to wards
union_m = unary_union(wards_metric.geometry.values)
roads_metric["geometry"] = roads_metric.geometry.intersection(union_m)
roads_metric = roads_metric[~roads_metric.geometry.is_empty &
                             (roads_metric.geometry.geom_type == "LineString")].copy()
print(f"  Roads in wards: {len(roads_metric)}", flush=True)

# Build boundary segments
parts = []
for idx, row in wards_metric.iterrows():
    wid = row["ward_id"]
    b = row.geometry.boundary
    if b.geom_type == "MultiLineString":
        for ln in b.geoms:
            parts.append((wid, ln))
    elif b.geom_type == "LineString":
        parts.append((wid, b))
print(f"  Segments: {len(parts)}", flush=True)

# Spatial index
sidx = roads_metric.sindex

# Find intersections
crossings = []
for wid, bline in parts:
    if bline.length < 1:
        continue
    cands = list(sidx.intersection(bline.bounds))
    for ridx in cands:
        rg = roads_metric.iloc[ridx].geometry
        if not bline.intersects(rg):
            continue
        inter = bline.intersection(rg)
        if inter.is_empty:
            continue
        if inter.geom_type == "Point":
            crossings.append({"ward_id": wid, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": wid, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": wid, "geometry": g})

print(f"  Raw: {len(crossings)}", flush=True)

if crossings:
    cg = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    cg["x"] = cg.geometry.x.round(2)
    cg["y"] = cg.geometry.y.round(2)
    cg = cg.drop_duplicates(subset=["ward_id", "x", "y"])
    cg = cg.drop(columns=["x", "y"])
    nm = {r["ward_id"]: (r["ward_name"], r["ward_name_en"]) for _, r in wards_metric.iterrows()}
    cg["ward_name"] = cg["ward_id"].map(lambda x: nm.get(x, ("", ""))[0])
    cg["ward_name_en"] = cg["ward_id"].map(lambda x: nm.get(x, ("", ""))[1])
    crossing_gdf = cg
    print(f"  Unique: {len(crossing_gdf)}", flush=True)
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"),
         "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings!", flush=True)

# ===============================================================
# STEP 4: Buffers
# ===============================================================
print("\n=== Step 4: 50m buffers ===", flush=True)

if len(crossing_gdf) > 0:
    bg = crossing_gdf.copy()
    bg["geometry"] = bg.geometry.buffer(50.0)
    buffer_gdf = bg.set_geometry("geometry")
    print(f"  {len(buffer_gdf)} buffers", flush=True)
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"),
         "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ===============================================================
# STEP 5: Clip buffers to ward
# ===============================================================
print("\n=== Step 5: Clipping to ward ===", flush=True)

if len(buffer_gdf) > 0:
    wgm = {r["ward_id"]: r.geometry for _, r in wards_metric.iterrows()}
    clip = []
    for _, br in buffer_gdf.iterrows():
        wg = wgm.get(br["ward_id"])
        if wg is not None:
            c = br.geometry.intersection(wg)
            clip.append(c if not c.is_empty else None)
        else:
            clip.append(None)
    cg = buffer_gdf.copy()
    cg["geometry"] = clip
    cg = cg[cg.geometry.notna() & ~cg.geometry.is_empty].copy()
    cg = cg.set_geometry("geometry")
    cg["geometry"] = cg.geometry.buffer(0)
    clipped_gdf = cg.set_geometry("geometry")
    print(f"  {len(clipped_gdf)} clipped", flush=True)
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"),
         "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ===============================================================
# STEP 6: Density
# ===============================================================
print("\n=== Step 6: Density ===", flush=True)

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

wa = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wa["ward_area_km2"] = wa.geometry.area / 1_000_000
dd = wa.merge(counts, on="ward_id", how="left")
dd["crossing_count"] = dd["crossing_count"].fillna(0).astype(int)
dd["crossings_per_km2"] = dd["crossing_count"] / dd["ward_area_km2"]

dg = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
dg = dg.merge(dd[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
dg = dg.to_crs(WGS84)
density_gdf = dg

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ===============================================================
# SAVE
# ===============================================================
print("\n=== Saving ===", flush=True)

wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print("  ✓ wards_jgd")

if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points ({len(crossing_gdf)})")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print("  ✓ crossing_points (empty)")

if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m ({len(buffer_gdf)})")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection ({len(clipped_gdf)})")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84 ({len(density_gdf)})")

print("\n=== DONE ===", flush=True)

# Verify
for lay in ["wards_jgd","crossing_points","crossing_buffers_50m","buffer_ward_intersection","ward_crossing_density_wgs84"]:
    try:
        g = gpd.read_file(OUTPUT, layer=lay)
        cs = str(g.crs).rsplit(":",1)[-1] if g.crs else "None"
        print(f"  {lay}: {len(g)} feat, EPSG:{cs}, cols={g.columns.tolist()}")
    except Exception as e:
        print(f"  {lay}: {e}")

#!/usr/bin/env python3
"""
Tokyo 23 Special Wards — Road crossing analysis pipeline.

Uses Overture Maps for ward boundaries (locality type) and roads (segment type).
Falls back to per-ward OSMnx for roads if Overture lacks detail.
Output: tokyo_crossings.gpkg (5 layers).
"""

import sys
import time
import warnings
from collections import Counter
from shapely.geometry import Point, LineString, MultiPoint, shape
from shapely.ops import unary_union
from shapely import wkt, wkb

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

warnings.filterwarnings("ignore")

METRIC_CRS = "EPSG:6675"
WGS84 = "EPSG:4326"
OUTPUT = "/work/tokyo_crossings.gpkg"

WARD_MAP = {
    "千代田区": "Chiyoda", "中央区": "Chuo", "港区": "Minato",
    "新宿区": "Shinjuku", "文京区": "Bunkyo", "台東区": "Taito",
    "墨田区": "Sumida", "江東区": "Koto", "品川区": "Shinagawa",
    "目黒区": "Meguro", "大田区": "Ota", "世田谷区": "Setagaya",
    "渋谷区": "Shibuya", "中野区": "Nakano", "杉並区": "Suginami",
    "豊島区": "Toshima", "北区": "Kita", "荒川区": "Arakawa",
    "板橋区": "Itabashi", "練馬区": "Nerima", "足立区": "Adachi",
    "葛飾区": "Katsushika", "江戸川区": "Edogawa",
}

# Tokyo bounding box (approx)
TOKYO_BBOX = (139.55, 35.55, 139.92, 35.82)


def overture_to_gdf(overture_type, bbox=TOKYO_BBOX):
    """Fetch Overture Maps data as a GeoDataFrame."""
    import overturemaps.core as oc
    reader = oc.record_batch_reader(overture_type, bbox=bbox)
    if reader is None:
        return gpd.GeoDataFrame()
    table = reader.read_all()
    df = table.to_pandas()
    # Decode WKB geometry
    df["geometry"] = df["geometry"].apply(
        lambda x: wkb.loads(bytes(x)) if x is not None else None
    )
    df = df[df["geometry"].notna()].copy()
    return gpd.GeoDataFrame(df, geometry="geometry", crs=WGS84)


# ===============================================================
# STEP 1: Fetch ward boundaries from Overture Maps
# ===============================================================
print("=== Step 1: Fetching ward boundaries (Overture) ===", flush=True)

# Overture has the wards as locality subtype with Japanese names
div = overture_to_gdf("division_area")
print(f"  Total division_area features: {len(div)}", flush=True)

# Filter for Japan
if "country" in div.columns:
    div = div[div["country"] == "JP"].copy()
    print(f"  Japan features: {len(div)}", flush=True)

# Filter to our 23 special wards by matching Japanese name
wards_list = []
for jp_name, en_name in WARD_MAP.items():
    # Check 'names' column
    for idx, row in div.iterrows():
        names = row.get("names", {})
        if isinstance(names, dict):
            primary = names.get("primary", "")
            common = names.get("common", {})
            if isinstance(common, dict):
                common_ja = common.get("ja", "")
            elif isinstance(common, list):
                common_ja = ""
            else:
                common_ja = ""
        else:
            primary = str(names) if names else ""
            common_ja = ""
        
        if primary == jp_name or common_ja == jp_name:
            geom = row.geometry
            if geom is not None and not geom.is_empty and geom.geom_type in ("Polygon", "MultiPolygon"):
                # Use the Overture ID as ward_id (it's a string, hash it)
                overture_id = row.get("id", f"ward_{jp_name}")
                wards_list.append({
                    "ward_id": hash(overture_id) % (10**10),  # numeric ID
                    "ward_name": jp_name,
                    "ward_name_en": en_name,
                    "geometry": geom,
                    "_overture_id": overture_id
                })
                print(f"  ✓ {jp_name} ({en_name})", flush=True)
                break
    else:
        print(f"  ✗ {jp_name} — not found in Overture", flush=True)

print(f"\n  Got {len(wards_list)} / 23 ward boundaries", flush=True)

if len(wards_list) == 0:
    print("CRITICAL: No ward boundaries found. Exiting.")
    sys.exit(1)

wards_gdf = gpd.GeoDataFrame(wards_list, crs=WGS84)
for _, r in wards_gdf.iterrows():
    print(f"    {r['ward_name_en']:>12s} (ID={r['ward_id']})")

# ===============================================================
# STEP 2: Fetch roads
# ===============================================================
print("\n=== Step 2: Fetching drivable roads ===", flush=True)

# First try Overture's transportation segments
print("  Fetching roads from Overture Maps...", flush=True)
try:
    seg = overture_to_gdf("segment")
    print(f"  Total segments: {len(seg)}", flush=True)
    
    # Filter to Japan - check country or region
    if "country" in seg.columns:
        seg = seg[seg["country"].isin(["JP", "JPN", ""])].copy() if any(seg["country"].notna()) else seg
    
    # Filter to drivable road types
    # Overture segment subtypes: road, highway, primary, secondary, etc.
    road_mask = pd.Series(True, index=seg.index)
    if "subtype" in seg.columns:
        drivable_types = {"road", "motorway", "trunk", "primary", "secondary", "tertiary", 
                          "motorway_link", "trunk_link", "primary_link", "secondary_link", "tertiary_link",
                          "unclassified", "residential", "living_street"}
        road_mask = seg["subtype"].apply(
            lambda x: str(x).lower() in drivable_types if pd.notna(x) else False
        )
    
    seg = seg[road_mask].copy()
    
    # Filter to LineString geometries
    seg = seg[seg.geometry.geom_type.isin(["LineString", "MultiLineString"])].copy()
    
    def to_ls(g):
        if g.geom_type == "MultiLineString":
            parts = list(g.geoms)
            return max(parts, key=lambda x: x.length) if parts else g
        return g
    
    seg["geometry"] = seg["geometry"].apply(to_ls)
    seg = seg[seg.geometry.geom_type == "LineString"].copy()
    
    print(f"  Drivable roads from Overture: {len(seg)}", flush=True)
    
    if len(seg) > 0:
        # Clip to ward union + small buffer
        union_buff = unary_union(wards_gdf.geometry.values).buffer(0.01)
        seg_clipped = seg[seg.geometry.intersects(union_buff)].copy()
        print(f"  Roads within Tokyo wards area: {len(seg_clipped)}", flush=True)
        
        if len(seg_clipped) > 0:
            all_road_geoms = list(seg_clipped.geometry.values)
        else:
            all_road_geoms = []
    else:
        all_road_geoms = []
        
except Exception as e:
    print(f"  Overture roads failed: {e}", flush=True)
    all_road_geoms = []

# If Overture didn't give us enough roads, use OSMnx per-ward with fallback
if len(all_road_geoms) < 100:
    print("  Getting roads from OSMnx (per-ward fallback)...", flush=True)
    import osmnx as ox
    
    MIRRORS = [
        "https://overpass-api.de/api",
        "https://overpass.kumi.systems/api",
        "https://lz4.overpass-api.de/api",
    ]
    ox.settings.timeout = 300
    ox.settings.max_query_area_size = 50 * 1000 * 1000
    ox.settings.bulk_download = True
    ox.settings.requests_kwargs = {"timeout": 300}
    
    for idx, row in wards_gdf.iterrows():
        jp, en = row["ward_name"], row["ward_name_en"]
        print(f"  {jp} ({en})...", end=" ", flush=True)
        
        got = False
        for mirror in MIRRORS:
            ox.settings.overpass_url = mirror
            for attempt in range(2):
                try:
                    G = ox.graph_from_polygon(row.geometry, network_type="drive",
                                               simplify=True, retain_all=True)
                    edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
                    if edges is not None and len(edges) > 0 and "geometry" in edges.columns:
                        edges = edges[edges.geometry.notna() & ~edges.geometry.is_empty].copy()
                        if "highway" in edges.columns:
                            def no_serv(hw):
                                if isinstance(hw, list):
                                    return not any("service" in str(h).lower() for h in hw)
                                return "service" not in str(hw).lower()
                            edges = edges[edges["highway"].apply(no_serv)].copy()
                        
                        def to_ls2(g):
                            if hasattr(g, 'geom_type') and g.geom_type == "MultiLineString":
                                parts = list(g.geoms)
                                return max(parts, key=lambda x: x.length) if parts else g
                            return g
                        edges["geometry"] = edges["geometry"].apply(to_ls2)
                        edges = edges[edges.geometry.geom_type == "LineString"].copy()
                        all_road_geoms.extend(edges.geometry.values)
                        print(f"{len(edges)} roads", flush=True)
                        got = True
                        break
                except Exception as e:
                    err = str(e)
                    if "refused" in err or "timed out" in err:
                        break
                    time.sleep(2)
            if got:
                break
        if not got:
            print("0 roads", flush=True)
    
    print(f"\n  Total from OSMnx: {len(all_road_geoms)}", flush=True)

# Deduplicate road geometries
seen = set()
unique_roads = []
for g in all_road_geoms:
    k = wkt.dumps(g, rounding_precision=6)
    if k not in seen:
        seen.add(k)
        unique_roads.append(g)

print(f"  Unique road segments: {len(unique_roads)}", flush=True)

if len(unique_roads) == 0:
    print("CRITICAL: No roads fetched. Exiting.")
    sys.exit(1)

roads_gdf = gpd.GeoDataFrame(geometry=unique_roads, crs=WGS84)

# ===============================================================
# STEP 3: Find crossing points
# ===============================================================
print("\n=== Step 3: Computing crossing points ===", flush=True)

print("  Projecting to metric CRS...", flush=True)
wards_metric = wards_gdf.to_crs(METRIC_CRS)
roads_metric = roads_gdf.to_crs(METRIC_CRS)

# Clip roads to ward union
union_m = unary_union(wards_metric.geometry.values)
roads_metric["geometry"] = roads_metric.geometry.intersection(union_m)
roads_metric = roads_metric[~roads_metric.geometry.is_empty &
                             (roads_metric.geometry.geom_type == "LineString")].copy()
print(f"  Roads in wards: {len(roads_metric)}", flush=True)

# Build boundary segments
parts = []
for idx, row in wards_metric.iterrows():
    wid = row["ward_id"]
    b = row.geometry.boundary
    if b.geom_type == "MultiLineString":
        for ln in b.geoms:
            parts.append((wid, ln))
    elif b.geom_type == "LineString":
        parts.append((wid, b))
print(f"  Segments: {len(parts)}", flush=True)

sidx = roads_metric.sindex
crossings = []
for wid, bline in parts:
    if bline.length < 1:
        continue
    cands = list(sidx.intersection(bline.bounds))
    for ridx in cands:
        rg = roads_metric.iloc[ridx].geometry
        if not bline.intersects(rg):
            continue
        inter = bline.intersection(rg)
        if inter.is_empty:
            continue
        if inter.geom_type == "Point":
            crossings.append({"ward_id": wid, "geometry": inter})
        elif inter.geom_type == "MultiPoint":
            for pt in inter.geoms:
                crossings.append({"ward_id": wid, "geometry": pt})
        elif inter.geom_type == "GeometryCollection":
            for g in inter.geoms:
                if g.geom_type == "Point":
                    crossings.append({"ward_id": wid, "geometry": g})

print(f"  Raw: {len(crossings)}", flush=True)

if crossings:
    cg = gpd.GeoDataFrame(crossings, crs=METRIC_CRS, geometry="geometry")
    cg["x"] = cg.geometry.x.round(2)
    cg["y"] = cg.geometry.y.round(2)
    cg = cg.drop_duplicates(subset=["ward_id", "x", "y"])
    cg = cg.drop(columns=["x", "y"])
    nm = {r["ward_id"]: (r["ward_name"], r["ward_name_en"]) for _, r in wards_metric.iterrows()}
    cg["ward_name"] = cg["ward_id"].map(lambda x: nm.get(x, ("", ""))[0])
    cg["ward_name_en"] = cg["ward_id"].map(lambda x: nm.get(x, ("", ""))[1])
    crossing_gdf = cg
    print(f"  Unique: {len(crossing_gdf)}", flush=True)
else:
    crossing_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"),
         "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )
    print("  No crossings!", flush=True)

# ===============================================================
# STEP 4: Buffers
# ===============================================================
print("\n=== Step 4: 50m buffers ===", flush=True)

if len(crossing_gdf) > 0:
    bg = crossing_gdf.copy()
    bg["geometry"] = bg.geometry.buffer(50.0)
    buffer_gdf = bg.set_geometry("geometry")
    print(f"  {len(buffer_gdf)} buffers", flush=True)
else:
    buffer_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"),
         "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ===============================================================
# STEP 5: Clip to ward
# ===============================================================
print("\n=== Step 5: Clipping to ward ===", flush=True)

if len(buffer_gdf) > 0:
    wgm = {r["ward_id"]: r.geometry for _, r in wards_metric.iterrows()}
    clip = []
    for _, br in buffer_gdf.iterrows():
        wg = wgm.get(br["ward_id"])
        clip.append(br.geometry.intersection(wg) if wg is not None else None)
    cg = buffer_gdf.copy()
    cg["geometry"] = [None if (c is None or c.is_empty) else c for c in clip]
    cg = cg[cg.geometry.notna() & ~cg.geometry.is_empty].copy()
    cg = cg.set_geometry("geometry")
    cg["geometry"] = cg.geometry.buffer(0)
    clipped_gdf = cg.set_geometry("geometry")
    print(f"  {len(clipped_gdf)} clipped", flush=True)
else:
    clipped_gdf = gpd.GeoDataFrame(
        {"ward_id": pd.Series(dtype="int64"), "ward_name": pd.Series(dtype="str"),
         "ward_name_en": pd.Series(dtype="str")},
        geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS
    )

# ===============================================================
# STEP 6: Density
# ===============================================================
print("\n=== Step 6: Density ===", flush=True)

if len(crossing_gdf) > 0:
    counts = crossing_gdf.groupby("ward_id").size().reset_index(name="crossing_count")
else:
    counts = pd.DataFrame({"ward_id": pd.Series(dtype="int64"), "crossing_count": pd.Series(dtype="int64")})

wa = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wa["ward_area_km2"] = wa.geometry.area / 1_000_000
dd = wa.merge(counts, on="ward_id", how="left")
dd["crossing_count"] = dd["crossing_count"].fillna(0).astype(int)
dd["crossings_per_km2"] = dd["crossing_count"] / dd["ward_area_km2"]

dg = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
dg = dg.merge(dd[["ward_id", "crossing_count", "ward_area_km2", "crossings_per_km2"]], on="ward_id")
dg = dg.to_crs(WGS84)
density_gdf = dg

print(density_gdf[["ward_id", "ward_name_en", "crossing_count", "ward_area_km2", "crossings_per_km2"]].to_string())

# ===============================================================
# SAVE
# ===============================================================
print("\n=== Saving ===", flush=True)

wr_out = wards_metric[["ward_id", "ward_name", "ward_name_en", "geometry"]].copy()
wr_out.to_file(OUTPUT, layer="wards_jgd", driver="GPKG")
print("  ✓ wards_jgd")

if len(crossing_gdf) > 0:
    crossing_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_points", driver="GPKG")
    print(f"  ✓ crossing_points ({len(crossing_gdf)})")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_points", driver="GPKG")

if len(buffer_gdf) > 0:
    buffer_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")
    print(f"  ✓ crossing_buffers_50m ({len(buffer_gdf)})")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="crossing_buffers_50m", driver="GPKG")

if len(clipped_gdf) > 0:
    clipped_gdf[["ward_id","ward_name","ward_name_en","geometry"]].to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")
    print(f"  ✓ buffer_ward_intersection ({len(clipped_gdf)})")
else:
    gpd.GeoDataFrame({"ward_id":[],"ward_name":[],"ward_name_en":[]}, geometry=gpd.GeoSeries(dtype="geometry"), crs=METRIC_CRS).to_file(OUTPUT, layer="buffer_ward_intersection", driver="GPKG")

density_gdf.to_file(OUTPUT, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f"  ✓ ward_crossing_density_wgs84 ({len(density_gdf)})")

print("\n=== DONE ===", flush=True)

# Verify
print("\n=== Verification ===")
for lay in ["wards_jgd","crossing_points","crossing_buffers_50m","buffer_ward_intersection","ward_crossing_density_wgs84"]:
    try:
        g = gpd.read_file(OUTPUT, layer=lay)
        cs = str(g.crs).rsplit(":",1)[-1] if g.crs else "None"
        print(f"  {lay}: {len(g)} feat, EPSG:{cs}, cols={g.columns.tolist()}")
    except Exception as e:
        print(f"  {lay}: {e}")