run-20260618-115720Z / crs-l3-tokyo-jgd-crossings

0:00

I'm rebuilding the road-safety dashboard for Tokyo's 23 special wards. For each ward I need to know how often the drivable road network crosses the ward boundary (please skip service roads), plus a 50 m buffer around every crossing for jurisdictional reporting. Pull the data from current OSM via Overpass. Run the engineering layers in the regional metric coordinate system so the distances and areas come out in real metres. Save everything to a single `tokyo_crossings.gpkg` with five layers: - `wards_jgd` (Polygon): the 23 ward boundaries. - `crossing_points` (Point): one point at each place a road crosses a ward boundary. - `crossing_buffers_50m` (Polygon): a 50 m buffer around each crossing point. - `buffer_ward_intersection` (Polygon): each buffer clipped to the ward whose boundary it crossed. - `ward_crossing_density_wgs84` (Polygon): one row per ward with `ward_id`, `crossing_count`, `ward_area_km2`, and `crossings_per_km2`. This is what the public dashboard reads, so it goes in WGS84 while the other four stay in the metric grid. Use each ward's OSM relation id as `ward_id`, and carry the English and native ward names as `ward_name_en` and `ward_name` on every layer that's tied to a specific ward.

0:00

System info

{
  "adapter": "openrouter",
  "model": "deepseek/deepseek-v4-flash",
  "base_url": "https://openrouter.ai/api/v1"
}

3:41

Let me break down this task:

1. Get the 23 special wards of Tokyo from OSM via Overpass
2. Get the drivable road network (excluding service roads) for Tokyo
3. Find where roads cross ward boundaries
4. Create 50m buffers around crossings
5. Clip buffers to the ward they crossed
6. Calculate crossing density per ward

Let me plan the approach:

1. Use Overpass API to fetch the 23 wards (they have relations with admin_level=9 or similar)
2. Use Overpass to fetch roads (highway=* excluding service roads)
3. Find intersections between roads and ward boundaries
4. Everything in a metric CRS - Tokyo uses JGD2011 / UTM zone 54N (EPSG:6676) or JGD2000 / UTM zone 54N. Let me check - Japan uses JGD2011/JGD2000. The appropriate UTM zone for Tokyo is 54N. JGD2011 / UTM zone 54N is EPSG:6676.

Actually, let me think about this more carefully. Tokyo is around 139.7°E longitude, which is UTM zone 54N (114°E to 120°E... wait, 54N covers 138°E to 144°E? Let me recalculate.)

UTM zones: 1 starts at 180°W going east, each zone is 6°. Zone 54: center at 54*6 - 183 = 324 - 183 = 141°E. So zone 54 covers 138°E to 144°E. Tokyo is at ~139.7°E, which is in zone 54. Yes.

So EPSG:6676 (JGD2011 / UTM zone 54N) or 3097 (JGD2000 / UTM zone 54N). Let me use JGD2011 - EPSG:6676.

Actually, let me use EPSG:6669 (JGD2011 / Japan Plane Rectangular CS IX) which is specific to Tokyo. CS IX covers Tokyo-to. Let me check - Japan Plane Rectangular Coordinate System:

Zone I: Kagoshima
Zone II: Fukuoka, Saga, etc.
Zone IX: Tokyo

EPSG:6674 is JGD2011 / Japan Plane Rectangular CS IX. Actually wait, let me verify:

Japan Plane Rectangular CS IX covers Tokyo-to. EPSG:6674 is JGD2011 / Japan Plane Rectangular CS IX. Yes, that's the one for Tokyo.

Actually, let me just use UTM 54N (EPSG:6676 for JGD2011) - it's simpler and covers the area well.

Let me start coding.

For Overpass queries:

1. Get the 23 special wards of Tokyo: these are relations with admin_level=9 or boundary=administrative, name:en contains "Special ward" or they have `is_in:city=Tokyo` or something. Actually, the 23 special wards of Tokyo are:
   - Adachi, Arakawa, Bunkyo, Chiyoda, Chuo, Edogawa, Itabashi, Katsushika, Kita, Koto, Meguro, Minato, Nakano, Nerima, Ota, Setagaya, Shibuya, Shinagawa, Shinjuku, Suginami, Sumida, Taito, Toshima

Let me query by relation IDs. Actually, I know the OSM relation IDs for Tokyo's wards:

Let me just use a query that gets all administrative boundaries with admin_level=9 or boundary=administrative within Tokyo. Or better, I can query specifically for the 23 special wards.

Actually, let's use a bounding box approach. Tokyo's 23 wards are roughly between 35.55°N to 35.8°N and 139.6°E to 139.9°E.

Let me write the Python script. I'll use osmnx to fetch the data since it handles Overpass well.

Actually, osmnx can be used to get administrative boundaries and road networks. Let me think about the approach.

For OSMnx:
- `osmnx.geometries_from_place` or `osmnx.geocode_to_gdf` for the wards
- `osmnx.graph_from_place` for the road network

Or I can use direct Overpass queries via requests.

Let me go with a direct Overpass approach since we need precise control over what we're fetching.

Let me write the Python script step by step:

1. Query Overpass for the 23 special ward boundaries
2. Query Overpass for drivable roads (excluding service roads) in the area
3. Find where roads cross ward boundaries
4. Create buffers and calculate densities

Let me write this.

For Overpass, I'll need to:
- Get all relations that are admin_level=9 within the Tokyo area (the 23 special wards)
- Get all highways except service roads within the same area

Let me think about the Overpass query for the wards:

```
[out:json][timeout:60];
area["name:en"="Tokyo"]["admin_level=4"]->.tokyo;
rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
out geom;
```

Wait, the special wards of Tokyo have admin_level=9. Let me verify... In Japan, admin levels are:
- 2: country
- 4: prefecture (Tokyo-to)
- 6: cities, special wards? No, special wards are admin_level=9 in OSM... 

Actually, let me check. In Japan OSM tagging, the 23 special wards of Tokyo are often tagged as admin_level=9 or admin_level=8. Let me check more carefully.

In OSM, Tokyo's 23 special wards are typically tagged with `admin_level=9` and `boundary=administrative`. But some might be admin_level=8. Let me use a broad query.

Actually, a simpler approach: use `osmnx.geocode_to_gdf("Tokyo, Japan")` to get Tokyo boundaries, then query for admin_level=9 boundaries within Tokyo.

Let me think of a simpler approach. I'll use Overpass to get the 23 wards directly by their relation IDs. Let me look them up.

Actually, I don't have the relation IDs memorized. Let me query for them.

Let me write the script using a combination of approaches. I'll use the Overpass API directly.

Here's my plan for the script:

```python
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPoint, MultiLineString, MultiPolygon
from shapely.ops import linemerge, unary_union
import requests
import json
import time

# Step 1: Fetch ward boundaries from Overpass
# Step 2: Fetch road network from Overpass
# Step 3: Find road-ward-boundary crossings
# Step 4: Create buffers
# Step 5: Calculate densities
# Step 6: Save to GPKG
```

Let me write the full script now.

For the Overpass query for wards, I'll use:

```
[out:json][timeout:120];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
(
  rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"]["name:en"];
  rel(area.tokyo)["admin_level"="8"]["boundary"="administrative"]["name:en"~"^(Adachi|Arakawa|Bunkyo|Chiyoda|Chuo|Edogawa|Itabashi|Katsushika|Kita|Koto|Meguro|Minato|Nakano|Nerima|Ota|Setagaya|Shibuya|Shinagawa|Shinjuku|Suginami|Sumida|Taito|Toshima)$"];
);
out geom;
```

Hmm, let me just use osmnx which is simpler. 

```python
import osmnx as ox
ox.config(use_cache=True, log_console=True)

# Get Tokyo special wards
# We can query by place
```

Let me try a different approach - just use direct Overpass queries with the `requests` library.

Actually, let me write a clean script. I'll use osmnx for the wards and roads.

Let me think about this more carefully.

For the 23 special wards, I can use osmnx's `ox.geometries_from_place` or `ox.geocode_to_gdf`.

For roads, I can use `ox.graph_from_place` with network_type='drive' which excludes service roads.

Let me write the script:

```python
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point
from shapely.ops import unary_union, linemerge
from shapely import wkt
import time
import requests

# Configure osmnx
ox.config(use_cache=True, log_console=True)

# Step 1: Get the 23 special wards
# Tokyo's 23 special wards can be fetched as administrative boundaries
# Let's use the place query for Tokyo and filter

# Actually, let's get wards by querying for the specific relations
# Use the Overpass API directly for more control
```

Let me just write a complete script using Overpass API requests directly. That gives me the most control.

Actually, I think the cleanest approach is to use osmnx to get the geometries. Let me try:

```python
import osmnx as ox
# Get Tokyo boundary first
tokyo = ox.geocode_to_gdf("Tokyo, Japan")
# Then get all admin boundaries within Tokyo
wards = ox.geometries_from_polygon(tokyo.geometry.iloc[0], tags={'boundary': 'administrative'})
```

But this might get more than we need. Let me be more specific.

OK let me just write a comprehensive script. I'll use Overpass API directly with queries.

For the metric CRS, I'll use EPSG:6676 (JGD2011 / UTM zone 54N) which covers the Tokyo area.

Let me write the script now:

```python
#!/usr/bin/env python3
"""
Tokyo 23 special wards road crossing analysis.
Fetches data from OSM via Overpass, computes crossings, buffers, and densities.
"""

import geopandas as gpd
import pandas as pd
import numpy as np
import requests
import json
import time
import sys
from shapely.geometry import shape, Point, LineString, MultiLineString, MultiPoint, box
from shapely.ops import unary_union, linemerge
from shapely import wkt

METRIC_CRS = "EPSG:6676"  # JGD2011 / UTM zone 54N
OUTPUT_GPKG = "/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",
]

def overpass_query(query, max_retries=3):
    """Run an Overpass query with retries."""
    for attempt in range(max_retries):
        for url in OVERPASS_URLS:
            try:
                resp = requests.post(url, data={'data': query}, timeout=180)
                if resp.status_code == 200:
                    return resp.json()
                print(f"  {url} returned {resp.status_code}", file=sys.stderr)
            except Exception as e:
                print(f"  {url} error: {e}", file=sys.stderr)
        if attempt < max_retries - 1:
            time.sleep(5)
    raise RuntimeError("All Overpass queries failed")

# Step 1: Fetch ward boundaries
print("Fetching ward boundaries from Overpass...")
wards_query = """
[out:json][timeout:120];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
(
  rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
);
out geom;
"""
wards_data = overpass_query(wards_query)

# Parse ward geometries
ward_features = []
for element in wards_data['elements']:
    tags = element.get('tags', {})
    name_en = tags.get('name:en', '')
    name = tags.get('name', '')
    ward_id = str(element['id'])
    
    # Get the geometry
    if 'members' in element:
        # For relations, construct the geometry from members
        ...
```

Wait, parsing relation geometries from Overpass JSON is complex. Let me use osmnx for this part instead.

Let me use a simpler approach with osmnx for fetching, and then do the analysis with geopandas.

```python
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, LineString, MultiLineString, box
from shapely.ops import unary_union
import sys

METRIC_CRS = "EPSG:6676"
OUTPUT_GPKG = "/work/tokyo_crossings.gpkg"
WGS84_CRS = "EPSG:4326"

ox.config(use_cache=True, log_console=True)

# Step 1: Fetch the 23 special wards
# Use osmnx to get administrative boundaries
print("Step 1: Fetching ward boundaries...")

# Get Tokyo boundaries to use as a clipping region
tokyo_gdf = ox.geocode_to_gdf("Tokyo, Japan")
tokyo_geom = tokyo_gdf.geometry.iloc[0]
tokyo_bounds = tokyo_geom.bounds  # (minx, miny, maxx, maxy)

# Get all admin boundaries within Tokyo
# The 23 special wards are admin_level=9 in OSM
wards = ox.geometries_from_bbox(
    tokyo_bounds[3], tokyo_bounds[1], tokyo_bounds[2], tokyo_bounds[0],
    tags={'boundary': 'administrative', 'admin_level': '9'}
)

# Filter to only keep polygons (the ward boundaries)
# ...
```

Hmm, actually `ox.geometries_from_bbox` with those tags might not directly get administrative boundary relations correctly. Let me think about this differently.

Actually, osmnx has a function `ox.geometries_from_place` that can get place boundaries.

Let me try:

```python
wards = ox.geometries_from_place("Tokyo, Japan", tags={'boundary': 'administrative'})
```

But this gives us all admin boundaries, not just the 23 wards.

Let me just use the Overpass API directly but with a simpler approach - fetch the ward boundaries as ways/relations with `out geom` and parse them.

Actually, the simplest approach might be:

1. Get ward boundaries from Overpass (which returns fully resolved geometries with `out geom`)
2. For roads, use osmnx to get the drivable network

Let me try the Overpass route for both.

For Overpass, when you use `out geom` on relations with `boundary=administrative` and `admin_level=9`, it returns the geometry as a polygon in the `geometry` field of the element if it's a multipolygon. Let me verify by checking what comes back.

Actually, for `rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"]; out geom;`, Overpass returns the relation elements. For multipolygon relations (which administrative boundaries are), the geometry is not directly in the element JSON but reconstructed from members when using `out geom`.

Let me use a different approach. Let me use `(rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];)->.wards; .wards out geom;` - this should give us the geometries.

Actually, in Overpass QL, when you do `out geom` on a relation that's a multipolygon/boundary, it reconstructs the geometry. The result JSON has `geometry` field on the element.

Let me just test this. Actually, let me write a script that first tries to get the data, and I'll handle the parsing carefully.

Let me write a comprehensive script:

For wards, I'll use this Overpass query:
```
[out:json][timeout:120];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
out geom;
```

For roads, I'll get all highways excluding service roads in the Tokyo area:
```
[out:json][timeout:180];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
(
  way(area.tokyo)["highway"]["highway"!="service"]["highway"!="proposed"]["highway"!="construction"]["highway"!="footway"]["highway"!="cycleway"]["highway"!="path"]["highway"!="track"]["highway"!="pedestrian"];
);
out geom;
```

Wait, I need to think about which highway values are "drivable". In OSM, drivable roads include:
- motorway, motorway_link
- trunk, trunk_link
- primary, primary_link
- secondary, secondary_link
- tertiary, tertiary_link
- unclassified
- residential
- living_street

And we want to exclude:
- service
- track
- footway
- cycleway
- path
- pedestrian
- steps
- proposed
- construction

So I should filter for these specific values. Or use a regex.

Let me use:
```
way(area.tokyo)["highway"~"^(motorway|motorway_link|trunk|trunk_link|primary|primary_link|secondary|secondary_link|tertiary|tertiary_link|unclassified|residential|living_street)$"];
```

This is more precise.

Now for parsing the Overpass JSON response:

For ways, each element has `geometry` which is a list of `{"lat": ..., "lon": ...}` points. I need to convert these to Shapely LineStrings.

For relations (wards), the geometry should be in the `geometry` field as well when using `out geom`. Let me check...

Actually, in Overpass JSON output, for relations with `out geom`, the geometry might be represented differently. Let me handle this carefully.

For a multipolygon relation, `out geom` should produce the geometry in the element. Let me check the structure...

Actually, I think the safest approach is to use osmnx for the wards since it handles the complex geometry parsing. Let me use osmnx.

```python
import osmnx as ox

# Get Tokyo boundary
tokyo = ox.geocode_to_gdf("Tokyo, Japan")
poly = tokyo.geometry.iloc[0]

# Get wards - admin_level=9 boundaries within Tokyo
# We can use osmnx's geometries_from_polygon
wards = ox.geometries_from_polygon(poly, tags={'boundary': 'administrative', 'admin_level': '9'})
```

But this might not get all 23 wards correctly since OSM tagging varies.

OK let me just go with a combined approach. I'll use osmnx where it works well and direct Overpass queries where needed.

Actually, let me just write a working script that uses both approaches and handles edge cases.

Let me look at this from a different angle. Let me use `osmnx.geocode_to_gdf` with specific queries for each ward.

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

ward_gdfs = []
for name in ward_names_en:
    gdf = ox.geocode_to_gdf(f"{name}, Tokyo, Japan")
    ward_gdfs.append(gdf)
   
wards = pd.concat(ward_gdfs)
```

But this might be slow with 23 separate geocode calls.

OK, I think the best approach is:

1. Use an Overpass query specifically for the 23 wards using their relation IDs or by querying for admin_level=9 in Tokyo
2. Parse the geometries from Overpass response
3. Use osmnx to get the road network (it handles graph construction, edge geometry, etc.)

Let me check if I can get the road network edges with osmnx:

```python
G = ox.graph_from_place("Tokyo, Japan", network_type='drive')
# Get edges as GeoDataFrame
edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
```

This excludes service roads by default for network_type='drive'. And we get a nice GeoDataFrame with edge geometries.

Then for the ward boundaries, I'll use Overpass.

Let me write the final script now.

Actually, let me also consider: "skip service roads" - in osmnx, `network_type='drive'` includes motorways, trunks, primary, secondary, tertiary, unclassified, residential, and living_street. It excludes service roads, tracks, pedestrian paths, etc. So that's perfect.

But wait, osmnx's network_type='drive' also includes motorway/trunk links, primary/secondary/tertiary links. That's fine for "drivable road network".

Let me now write the complete script.

Here's my final approach:

```python
#!/usr/bin/env python3
"""
Tokyo 23 special wards road crossing analysis.
"""

import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
import requests
import json
import time
import sys
from shapely.geometry import shape, Point, LineString, MultiLineString, box
from shapely.ops import unary_union, linemerge, split, snap
from shapely import wkt

METRIC_CRS = "EPSG:6676"
WGS84 = "EPSG:4326"
OUTPUT_GPKG = "/work/tokyo_crossings.gpkg"

ox.config(use_cache=True, log_console=True)

# Step 1: Get ward boundaries from Overpass
# Step 2: Get road network from osmnx
# Step 3: Find crossings
# ...

# ... etc
```

Let me just write the whole thing now and test it. Let me be very careful about the Overpass query syntax.

For ward boundaries, let me use this Overpass query:

```
[out:json][timeout:120];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
(
  rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
);
out geom;
```

This gets all relations with admin_level=9 and boundary=administrative within the Tokyo prefecture area. The 23 special wards of Tokyo should all have this tagging.

Actually wait, I should double check. In OSM, the 23 special wards might also be tagged as admin_level=8 in some cases. Let me check:
- Most Tokyo special wards in OSM are tagged admin_level=9
- Some might be tagged differently

Let me use a broader query:

```
[out:json][timeout:120];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
(
  rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
  rel(area.tokyo)["admin_level"="8"]["boundary"="administrative"]["name:en"~"^(Adachi|Arakawa|Bunkyo|Chiyoda|Chuo|Edogawa|Itabashi|Katsushika|Kita|Koto|Meguro|Minato|Nakano|Nerima|Ota|Setagaya|Shibuya|Shinagawa|Shinjuku|Suginami|Sumida|Taito|Toshima)$"];
);
out geom;
```

Hmm, but that regex is complex. Let me just fetch admin_level=9 and also possibly admin_level=8, then filter by name after.

Actually, let me just try fetching admin_level=9 first and see what we get.

OK let me write the complete script now. I'll handle the Overpass responses carefully.

For parsing the ward boundaries from Overpass:
- Relations with `out geom` return the geometry as part of the element
- For multipolygon relations, the geometry is typically in the `geometry` field as a list of polygon rings
- I need to construct Shapely polygons from these

Let me write the full script now and run it.

One more thing - for the road network from osmnx, I need to be careful about the area. The 23 wards occupy a specific area of Tokyo. If I query for the entire Tokyo prefecture, I'll get roads outside the 23 wards. So I should clip to the ward boundaries.

Let me structure the script:

1. Fetch ward boundaries from Overpass
2. Create a unified boundary of the 23 wards
3. Fetch road network using osmnx for the Tokyo area
4. Clip roads to the unified ward boundary
5. Split roads at ward boundaries to find crossings
6. Create buffers and calculate densities

Let me write this now.

Actually, one important consideration: to find where roads cross ward boundaries, I need to:

1. For each pair of adjacent wards, find the shared boundary
2. Find where roads intersect that shared boundary
3. The intersection points are the crossings

OR

1. Extract all ward boundary lines (the borders between wards)
2. Find where roads intersect these boundary lines
3. The intersection points are the crossings

The second approach is simpler. Let me use that.

Actually, the exact approach:
1. Get all ward boundary polygons
2. For each ward, extract its boundary as a line
3. Find all roads that cross the ward boundary line
4. The intersection points are the crossing points

Wait, but each ward boundary line includes both the outer boundary (between the ward and non-ward areas) and inner boundaries (between wards). The crossings we care about are where a road crosses from one ward to another, or from a ward to the outside. So we need the full ward boundaries.

Let me think again. "how often the drivable road network crosses the ward boundary" - this means for each ward, we count how many times a road crosses into or out of the ward. So we need the full boundary of the ward (the polygon's exterior ring).

So for each ward:
1. Extract the boundary line of the ward polygon
2. Find all roads that intersect this boundary line
3. For each intersection, create a point at the crossing location

But we need to be careful - if two adjacent wards share a boundary, a road crossing that boundary will be counted for both wards (once for each side of the crossing). That seems correct per the problem statement.

Actually, let me re-read: "how often the drivable road network crosses the ward boundary". For each ward, count how many times a road crosses its boundary. Yes, this is what I described.

Let me code this up.

For the implementation:

```python
# For each ward, get its boundary line
for idx, ward in wards.iterrows():
    boundary = ward.geometry.boundary  # This is a LineString or MultiLineString
    
    # Find all roads that intersect the boundary
    intersecting_roads = roads[roads.intersects(boundary)]
    
    for road_idx, road in intersecting_roads.iterrows():
        # Find intersection points
        intersection = road.geometry.intersection(boundary)
        if intersection.is_empty:
            continue
        if isinstance(intersection, (MultiPoint,)):
            for pt in intersection.geoms:
                crossing_points.append({'ward_id': ward.ward_id, 'geometry': pt})
        elif isinstance(intersection, Point):
            crossing_points.append({'ward_id': ward.ward_id, 'geometry': intersection})
        # Also handle cases where intersection is a LineString (road follows the boundary)
        # Skip those - we only want actual crossings
```

Wait, but if a road runs along a ward boundary for some distance, the intersection could be a LineString. In that case, the road doesn't really "cross" the boundary - it follows it. I need to handle this case.

Also, if a road goes from a ward to the outside (the coast or the edge of the 23-ward area), that's still a crossing. So we want all boundary intersections.

Let me think about this more carefully. The boundary of a ward polygon has:
- The outer ring (exterior) - separates the ward from outside areas (could be other wards, the sea, or the rest of Tokyo)
- Inner rings (interiors) - if the ward has holes (unlikely for Tokyo wards)

A road crossing the boundary is a road that goes from inside the ward to outside the ward (or vice versa). So I need to find where roads intersect the boundary.

For each road segment (edge in osmnx), if it intersects the ward boundary, the intersection geometry could be:
- A Point (road cleanly crosses the boundary)
- A MultiPoint (road crosses the boundary multiple times, e.g., a winding road)
- A LineString (road follows the boundary for some distance)

For LineString intersections, I should probably compute the start and end points of the road and check which side of the boundary they're on. If they're on different sides, the road crosses the boundary (but the intersection is just a line along the boundary).

Actually, for simplicity, let me handle it as follows:
- For each road and each ward, compute the intersection of the road with the ward boundary
- If the intersection is a Point, add it as a crossing
- If the intersection is a MultiPoint, add each point as a crossing
- If the intersection is a LineString (road follows the boundary), check if the road endpoints are on different sides of the boundary. If so, add the midpoint or endpoints as crossings.
- If the intersection is a GeometryCollection, extract Points and handle LineStrings similarly.

Actually, a simpler approach: for each road segment, check if it crosses the ward boundary by checking if the road's centroid is within the ward and an endpoint is outside, or vice versa. Then find the exact intersection point.

Hmm, this is getting complex. Let me just use a straightforward approach:

1. For each ward, get the boundary as a line
2. For each road, get the intersection with the boundary
3. Filter intersection geometries to only keep Points and MultiPoints (skip LineStrings)
4. Add those as crossing points

For LineString intersections (road follows the boundary), we can check if the road's endpoints are on different sides. If so, the start and end of the LineString intersection are the crossing points. But this is an edge case that might not happen often with well-mapped roads.

Actually, a better approach: for each road, find all points where it intersects the ward boundary. Use `unary_union` or split approach.

Simplest approach that handles the common case:

```python
# Extract ward boundary as a line
ward_boundary = ward.geometry.boundary

# For each road, find intersection
intersection = road.geometry.intersection(ward_boundary)

# Extract points
if intersection.geom_type == 'Point':
    points = [intersection]
elif intersection.geom_type == 'MultiPoint':
    points = list(intersection.geoms)
elif intersection.geom_type == 'GeometryCollection':
    points = [g for g in intersection.geoms if g.geom_type == 'Point']
elif intersection.geom_type in ('LineString', 'MultiLineString'):
    # Road follows boundary - check endpoints
    # If the road crosses (endpoints on different sides), take the intersection endpoints
    coords = list(intersection.coords if intersection.geom_type == 'LineString' else [c for ls in intersection.geoms for c in ls.coords])
    if coords:
        points = [Point(c) for c in coords]  # just use boundary crossing points
    else:
        points = []
else:
    points = []
```

Actually, let me just go with a simpler approach. Let me use the Overpass `around` operator or network approach.

You know what, let me just write the code and handle edge cases as they come up. The core logic is straightforward.

Let me finalize the script now.

For the road network, osmnx gives us edge geometries as LineStrings. Each edge goes from one node to another. If an edge crosses a ward boundary, I can find the intersection point(s).

Here's my refined approach:

```python
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPoint, LineString
from shapely.ops import unary_union
import requests
import json
import time
import sys

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

# Overpass query for ward boundaries
OVERPASS_URL = "https://overpass-api.de/api/interpreter"

def fetch_wards():
    """Fetch the 23 special ward boundaries from Overpass."""
    query = """
    [out:json][timeout:120];
    area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
    (
      rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
    );
    out geom;
    """
    resp = requests.post(OVERPASS_URL, data={'data': query}, timeout=180)
    resp.raise_for_status()
    data = resp.json()
    
    # Parse geometries from Overpass response
    features = []
    for el in data['elements']:
        if el['type'] != 'relation':
            continue
        tags = el.get('tags', {})
        name_en = tags.get('name:en', '')
        name = tags.get('name', '')
        ward_id = str(el['id'])
        
        # Get geometry - for boundaries, the geometry is reconstructed
        geom = None
        if 'geometry' in el:
            # This might be a list of polygon rings for multipolygon
            geom_data = el['geometry']
            # ...
        
        features.append({
            'ward_id': ward_id,
            'ward_name_en': name_en,
            'ward_name': name,
            'geometry': geom
        })
    
    return gpd.GeoDataFrame(features, crs=WGS84)
```

Hmm, the geometry parsing from Overpass JSON is tricky. Let me use an alternative approach.

I can use osmnx to get the ward boundaries directly:

```python
import osmnx as ox

# Get Tokyo boundaries
tokyo = ox.geocode_to_gdf("Tokyo, Japan")

# Get all admin_level=9 boundaries within Tokyo
# ox.geometries_from_bbox or ox.geometries_from_polygon
```

Let me try this approach. `ox.geometries_from_polygon` can fetch OSM elements within a polygon.

```python
tokyo_poly = tokyo.geometry.iloc[0]
wards = ox.geometries_from_polygon(tokyo_poly, tags={'boundary': 'administrative', 'admin_level': '9'})
```

This should give us the ward boundaries as polygons with their tags.

Let me write the complete script now and test it.

Let me also check: does osmnx's `graph_from_place` handle the 23 wards area properly? I should pass a query that confines to the 23 wards, not all of Tokyo prefecture.

Actually, I'll:
1. First get the ward boundaries
2. Create a dissolved boundary of all 23 wards
3. Use that as the polygon for osmnx road fetching

Let me write the script:

```python
#!/usr/bin/env python3
"""Tokyo 23 wards road crossing analysis."""

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

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

ox.config(use_cache=True, log_console=True)

# === STEP 1: Fetch ward boundaries ===
print("Fetching ward boundaries...")
tokyo = ox.geocode_to_gdf("Tokyo, Japan")
tokyo_poly = tokyo.geometry.iloc[0]

wards = ox.geometries_from_polygon(tokyo_poly, tags={
    'boundary': 'administrative',
    'admin_level': '9'
})

print(f"Found {len(wards)} ward boundaries")

# Filter to only keep polygon geometries
wards = wards[wards.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()

# Rename/clean columns
wards['ward_id'] = wards['osm_id'].astype(str)
wards['ward_name_en'] = wards['name:en'] if 'name:en' in wards.columns else wards.get('name', '')
wards['ward_name'] = wards['name'] if 'name' in wards.columns else ''
```

Wait, I need to be careful with osmnx's output. The columns and their names depend on the version and the data. Let me handle this more carefully when writing the script.

Let me also think about whether osmnx will actually be able to fetch all 23 special wards. The special wards in OSM are typically tagged as:
- boundary=administrative
- admin_level=9
- They're relations (type=boundary)

osmnx should be able to get these via `ox.geometries_from_polygon`.

OK, let me write the complete script now and run it. I'll handle errors as they come.

Actually, you know what, let me just write the script and test it. Let me be practical.

Here's the plan for the script:

1. Fetch 23 ward boundaries using osmnx or Overpass
2. Fetch drivable road network using osmnx (network_type='drive')
3. Reproject to metric CRS (EPSG:6676)
4. For each ward, find road crossings of its boundary
5. Create 50m buffers around crossing points
6. Clip buffers to the ward
7. Calculate crossing densities
8. Save to GPKG

Let me write it now.

One concern: osmnx.graph_from_place with "Tokyo, Japan" will fetch roads for the entire Tokyo prefecture, which is much larger than the 23 wards. I should clip to the 23 wards area.

Better approach: use the 23 wards boundaries (dissolved) as the query polygon for osmnx.

Or use `ox.graph_from_polygon` if available, or use `ox.graph_from_bbox`.

Actually, I can use `ox.graph_from_place` with a more specific place name, or I can use `ox.graph_from_polygon` which takes a polygon.

Wait, osmnx has `ox.graph_from_polygon` which accepts a polygon to limit the query area. Let me use that with the dissolved ward boundaries.

Let me finalize the script:

```python
# Get ward boundaries
# ... (using osmnx or Overpass)

# Create a unified boundary of all 23 wards
unified_boundary = unary_union(wards.geometry.values)

# Get road network within the unified boundary
G = ox.graph_from_polygon(unified_boundary, network_type='drive')
edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
```

This should work well.

For crossing detection, let me use a spatial join approach:
1. For each ward, extract its boundary line
2. Find all road edges that intersect the boundary line
3. Compute the intersection point(s)

Let me code this efficiently:

```python
# For each ward
crossings = []
for idx, ward in wards.iterrows():
    boundary = ward.geometry.boundary  # ward boundary as line
    # Find roads intersecting this boundary
    mask = edges.intersects(boundary)
    if mask.any():
        intersecting_roads = edges[mask]
        for road_idx, road in intersecting_roads.iterrows():
            inter = road.geometry.intersection(boundary)
            if inter.is_empty:
                continue
            # Extract points
            if inter.geom_type == 'Point':
                crossings.append({'ward_id': ward['ward_id'], 'ward_name_en': ward['ward_name_en'], 'ward_name': ward['ward_name'], 'geometry': inter})
            elif inter.geom_type == 'MultiPoint':
                for pt in inter.geoms:
                    crossings.append({'ward_id': ward['ward_id'], 'ward_name_en': ward['ward_name_en'], 'ward_name': ward['ward_name'], 'geometry': pt})
            elif inter.geom_type == 'GeometryCollection':
                for geom in inter.geoms:
                    if geom.geom_type == 'Point':
                        crossings.append({'ward_id': ward['ward_id'], 'ward_name_en': ward['ward_name_en'], 'ward_name': ward['ward_name'], 'geometry': geom})
```

But this could miss cases where the road runs along the boundary. Let me handle that by also checking if the road's endpoints are on different sides of the boundary.

Actually, for simplicity and the common case, let me just handle Point and MultiPoint for now. The LineString case (road follows boundary) is rare and would mean the road IS the boundary, which is unusual for mapped roads.

Let me write the full script now.

Wait, I also need to consider performance. For 23 wards and potentially thousands of road edges, doing O(n*m) intersection checks could be slow. Let me use spatial indexing.

Actually, geopandas has sjoin which uses spatial indexing. Let me use that.

Approach:
1. Create a GeoSeries of ward boundary lines
2. Find which roads intersect which boundary lines using spatial join
3. Compute intersection points

Let me restructure:

```python
# Create ward boundary lines
ward_boundaries = wards.geometry.boundary  # Series of LineStrings/MultiLineStrings
ward_boundaries_gdf = gpd.GeoDataFrame(wards[['ward_id', 'ward_name_en', 'ward_name']], 
                                        geometry=ward_boundaries, crs=wards.crs)

# Reproject everything to metric
wards_metric = wards.to_crs(METRIC_CRS)
edges_metric = edges.to_crs(METRIC_CRS)
ward_boundaries_metric = ward_boundaries_gdf.to_crs(METRIC_CRS)

# Spatial join: which roads intersect which ward boundaries
# But this gives us pairs (road, ward) where the road crosses the boundary
potential = gpd.sjoin(edges_metric, ward_boundaries_metric, how='inner', predicate='intersects')
```

Hmm, but `sjoin` with the boundary line might not work well because the boundary is a LineString and we're checking intersection with road edges. Let me think...

Actually, `gpd.sjoin` can join based on intersection. If we have roads (LineStrings) and ward boundary lines (LineStrings), it will join them where they intersect. That should work.

But wait, sjoin works with the full geometry, so a road that's entirely within a ward wouldn't intersect the boundary. Good, that's what we want.

But a road might run along the boundary for some distance, in which case they'd intersect as a LineString. That's a special case.

Let me just write the script and handle the common case. I can refine later if needed.

Actually, let me take a step back and think about what osmnx returns for edges. Each edge has a geometry (LineString) between two nodes. The edges are typically short segments between intersections. So the chance of an edge running along a ward boundary for a significant distance is very low. The intersection will almost always be a Point or empty.

OK let me just write the full script now.

Let me also think about the buffer and clipping:

```python
# Create 50m buffers
buffers = crossing_points_metric.copy()
buffers['geometry'] = buffers.geometry.buffer(50)
buffers = buffers.set_crs(METRIC_CRS)

# Clip buffers to their respective wards
# For each buffer, clip to the ward polygon
ward_intersections = []
for idx, buf in buffers.iterrows():
    ward_id = buf['ward_id']
    ward_geom = wards_metric[wards_metric['ward_id'] == ward_id].geometry.iloc[0]
    clipped = buf.geometry.intersection(ward_geom)
    ward_intersections.append(clipped)

buffer_ward_intersection = buffers.copy()
buffer_ward_intersection['geometry'] = ward_intersections
# Remove empty geometries
buffer_ward_intersection = buffer_ward_intersection[~buffer_ward_intersection.is_empty]
```

For the ward_crossing_density_wgs84 layer:
```python
# Count crossings per ward
crossing_counts = crossing_points.groupby('ward_id').size().reset_index(name='crossing_count')

# Merge with ward areas
wards_metric = wards.to_crs(METRIC_CRS)
wards_metric['ward_area_km2'] = wards_metric.geometry.area / 1_000_000

density = wards.merge(crossing_counts, on='ward_id', how='left')
density = density.merge(wards_metric[['ward_id', 'ward_area_km2']], on='ward_id')
density['crossing_count'] = density['crossing_count'].fillna(0).astype(int)
density['crossings_per_km2'] = density['crossing_count'] / density['ward_area_km2']
density = density.to_crs(WGS84)
```

OK I think I have a clear plan. Let me write the full script now.

But wait - there's a subtlety with the Overpass/osmnx fetching. osmnx's `ox.geometries_from_polygon` fetches OSM elements (ways/relations) that intersect the polygon. For admin boundaries (relations), it should return them as polygons.

Let me also consider: what if osmnx doesn't return all 23 wards? I might need to handle that case.

Let me also think about the order of wards. I hard-coded the list of 23 ward names in English, which I can use for validation.

OK, let me just write the script and run it. I'll debug as needed.

Let me also think about the Overpass query for the wards one more time. 

Actually, I just realized there might be an issue with osmnx and fetching admin boundaries. osmnx is primarily designed for street networks. Its `geometries_from_polygon` function can fetch OSM elements by tags, but it might not handle boundary relations perfectly.

Let me use a direct Overpass query for the wards and osmnx for the roads.

For the Overpass ward query, I need to properly parse the relation geometries. Let me handle this by using `out geom` on both the relation and its member ways.

Actually, the simplest approach for parsing boundary relation geometries from Overpass:

Use `out geom;` which returns the geometry of each element. For a multipolygon relation representing a boundary, the element should have a `geometry` field that's a list of coordinate arrays representing the polygon rings.

Let me check... In Overpass JSON output:

For a relation with `out geom`, the JSON looks like:
```json
{
  "type": "relation",
  "id": 12345,
  "members": [...],
  "tags": {...},
  "geometry": [{"lat": 35.6, "lon": 139.7}, ...]
}
```

Wait, the `geometry` for a relation is a flat list of coordinates that represents the boundary as a closed ring (for simple boundaries) or multiple rings (for multipolygons). But actually, for multipolygon boundaries, the geometry is more complex.

Let me use a different Overpass approach. I can use `make` or `convert` statements, or I can use the `out body;` and then construct geometries from ways.

Actually, the cleanest Overpass approach is to use:

```
[out:json][timeout:120];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
(
  rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
);
out geom;
```

And then parse the result. For relations with `out geom`, the geometry field should contain the reconstructed boundary.

Let me test this approach. If it doesn't work, I'll fall back to using osmnx.

Let me write a comprehensive script that handles this.

Actually, I just realized I can use a different approach altogether. Let me use `osmnx.geocode_to_gdf` for each ward individually. That way I know I'm getting correct geometries.

```python
ward_names = [
    "Adachi, Tokyo, Japan",
    "Arakawa, Tokyo, Japan",
    ...
]

ward_gdfs = []
for name in ward_names:
    gdf = ox.geocode_to_gdf(name)
    ward_gdfs.append(gdf)

wards = pd.concat(ward_gdfs, ignore_index=True)
```

This is slow (23 API calls) but reliable. However, each `geocode_to_gdf` call makes a Nominatim request which might be rate-limited.

OK, let me go with the Overpass query approach and parse it carefully.

Let me write the script now. I'll handle the Overpass response robustly.

Actually, I just thought of another approach: use `osmnx.geometries_from_place` with Tokyo and filter for admin level 9.

```python
wards = ox.geometries_from_place("Tokyo, Japan", tags={
    'boundary': 'administrative', 
    'admin_level': '9'
})
```

This should work reliably. Let me use this approach.

Then for the road network:
```python
G = ox.graph_from_polygon(unified_ward_boundary, network_type='drive')
edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
```

Let me write the complete script now.

One more thing: the `ox.graph_from_polygon` function accepts a Shapely polygon and fetches the street network within it. The `network_type='drive'` parameter ensures we only get drivable roads (excluding service roads).

Let me also think about the CRS. The task says:
- "Run the engineering layers in the regional metric coordinate system so the distances and areas come out in real metres."
- "crossing_buffers_50m (Polygon): a 50 m buffer around each crossing point" - so 50m in real metres
- "ward_crossing_density_wgs84 (Polygon) ... goes in WGS84 while the other four stay in the metric grid."

So the metric CRS should be:
- wards_jgd - metric
- crossing_points - metric
- crossing_buffers_50m - metric
- buffer_ward_intersection - metric
- ward_crossing_density_wgs84 - WGS84

The task says "wards_jgd" - the "jgd" suffix suggests Japan Geodetic Datum. I'll use EPSG:6676 (JGD2011 / UTM zone 54N) or EPSG:6674 (JGD2011 / Japan Plane Rectangular CS IX).

Let me check: Japan Plane Rectangular CS IX covers Tokyo-to. Zone IX's origin is at 36°N, 139°50'E. The scale factor at the origin is 0.9999. This is more accurate for Tokyo than UTM.

Let me use EPSG:6674 (JGD2011 / Japan Plane Rectangular CS IX). Or wait, let me check the exact EPSG codes.

JGD2011 / Japan Plane Rectangular CS IX - EPSG:6674
JGD2011 / UTM zone 54N - EPSG:6676

Both should work. Let me use EPSG:6676 (UTM zone 54N) as it's more commonly used.

Actually, let me use EPSG:6674 (Japan Plane Rectangular CS IX) since it's specifically designed for the Tokyo area and will give more accurate distances.

Wait, actually, let me just use a well-known metric CRS. The task mentions "the regional metric coordinate system" - I'll use EPSG:6676 (JGD2011 / UTM zone 54N).

Actually the user said "jgd" in the layer name `wards_jgd`. This suggests they want JGD (Japan Geodetic Datum). Let me use JGD2011.

EPSG:6676 - JGD2011 / UTM zone 54N

This is a metric CRS. Let me use this.

Let me finalize and write the script now.

Let me write the script step by step.

For the wards, I'll use Overpass directly to have full control. Let me handle the parsing of relation geometries carefully.

In Overpass, when you query `rel(...)["admin_level"="9"]["boundary"="administrative"]; out geom;`, the result elements have a `geometry` field (for ways in the relation that have been assembled). But actually, for relations, the `geometry` field in `out geom` is not always what we expect.

Looking at the Overpass documentation more carefully:

With `out geom`, for each element, the geometry is included as a list of coordinate pairs. For a way, this is straightforward. For a relation, the geometry is... complex.

Actually, let me look at this differently. The Overpass `out geom` on a relation that's tagged as a boundary (multipolygon) should reconstruct the multipolygon geometry. But the JSON format represents this as a single array of coordinates, which might be the outer ring.

Let me use a different approach. I'll fetch the member ways of the boundary relations and construct the polygon geometry myself.

Or, even simpler: use the Overpass `convert` statement to get the geometry as a GeoJSON feature.

Actually, the simplest reliable approach: use the Overpass `out geom` on ways within the boundary relation, or just use osmnx which handles all this complexity.

Let me just use osmnx to get the wards. It's more reliable.

Let me write the final script:

```python
#!/usr/bin/env python3
"""
Tokyo 23 special wards road crossing analysis.
Uses OSM data via osmnx for roads and Overpass API for ward boundaries.
"""

import os
import sys
import time
import requests
import json

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

# ============================================================
# CONFIGURATION
# ============================================================
WGS84 = "EPSG:4326"
METRIC_CRS = "EPSG:6676"  # JGD2011 / UTM zone 54N
OUTPUT_GPKG = "/work/tokyo_crossings.gpkg"

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

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

ox.config(use_cache=True, log_console=True)


def overpass_query(query, max_retries=3):
    """Run Overpass query with retries."""
    urls = [OVERAPASS_URL] + OVERAPASS_BACKUP_URLS
    for attempt in range(max_retries):
        for url in urls:
            try:
                resp = requests.post(url, data={"data": query}, timeout=300)
                if resp.status_code == 200:
                    return resp.json()
                print(f"  {url} HTTP {resp.status_code}", file=sys.stderr, flush=True)
            except Exception as e:
                print(f"  {url}: {e}", file=sys.stderr, flush=True)
        if attempt < max_retries - 1:
            wait = 5 * (attempt + 1)
            print(f"  Retrying in {wait}s...", file=sys.stderr, flush=True)
            time.sleep(wait)
    raise RuntimeError("All Overpass query attempts failed")


def fetch_wards_overpass():
    """Fetch 23 special ward boundaries via Overpass API."""
    print("Fetching ward boundaries from Overpass...", flush=True)
    
    # Query for all admin_level=9 boundaries in Tokyo
    query = """
    [out:json][timeout:120];
    area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
    (
      rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
    );
    out geom;
    """
    
    data = overpass_query(query)
    
    rows = []
    for el in data.get("elements", []):
        if el["type"] != "relation":
            continue
        
        tags = el.get("tags", {})
        ward_id = str(el["id"])
        ward_name_en = tags.get("name:en", "")
        ward_name = tags.get("name", "")
        
        # Skip non-Tokyo wards
        if ward_name_en not in TOKYO_WARDS_EN:
            continue
        
        # Build geometry from the relation geometry
        geom = None
        if "geometry" in el and len(el["geometry"]) > 0:
            coords = [(p["lon"], p["lat"]) for p in el["geometry"]]
            if len(coords) >= 3:
                geom = Polygon(coords)
                if not geom.is_valid:
                    geom = geom.buffer(0)
        
        if geom is not None:
            rows.append({
                "ward_id": ward_id,
                "ward_name_en": ward_name_en,
                "ward_name": ward_name,
                "geometry": geom
            })
    
    print(f"  Got {len(rows)} ward boundaries", flush=True)
    return gpd.GeoDataFrame(rows, crs=WGS84)


def main():
    # Step 1: Get ward boundaries
    wards = fetch_wards_overpass()
    
    # Check we got all 23 wards
    missing = set(TOKYO_WARDS_EN) - set(wards["ward_name_en"])
    if missing:
        print(f"WARNING: Missing wards: {missing}", flush=True)
    
    # Step 2: Get road network
    print("Fetching road network from osmnx...", flush=True)
    unified = unary_union(wards.geometry.values)
    
    G = ox.graph_from_polygon(unified, network_type='drive', simplify=False)
    edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
    print(f"  Got {len(edges)} road edges", flush=True)
    
    # Step 3: Reproject to metric CRS
    wards_metric = wards.to_crs(METRIC_CRS)
    edges_metric = edges.to_crs(METRIC_CRS)
    
    # Step 4: Find crossing points
    print("Finding boundary crossings...", flush=True)
    
    crossing_rows = []
    
    for idx, ward in wards_metric.iterrows():
        ward_id = ward["ward_id"]
        ward_name_en = ward["ward_name_en"]
        ward_name = ward["ward_name"]
        boundary = ward.geometry.boundary
        
        # Get all roads that intersect this ward boundary
        mask = edges_metric.intersects(boundary)
        intersecting_roads = edges_metric[mask]
        
        for road_idx, road in intersecting_roads.iterrows():
            inter = road.geometry.intersection(boundary)
            if inter.is_empty:
                continue
            
            if inter.geom_type == "Point":
                crossing_rows.append({
                    "ward_id": ward_id,
                    "ward_name_en": ward_name_en,
                    "ward_name": ward_name,
                    "geometry": inter
                })
            elif inter.geom_type == "MultiPoint":
                for pt in inter.geoms:
                    crossing_rows.append({
                        "ward_id": ward_id,
                        "ward_name_en": ward_name_en,
                        "ward_name": ward_name,
                        "geometry": pt
                    })
            elif inter.geom_type == "GeometryCollection":
                for g in inter.geoms:
                    if g.geom_type == "Point":
                        crossing_rows.append({
                            "ward_id": ward_id,
                            "ward_name_en": ward_name_en,
                            "ward_name": ward_name,
                            "geometry": g
                        })
                    elif g.geom_type in ("LineString", "MultiLineString"):
                        # Road follows boundary, take endpoints
                        if g.geom_type == "LineString":
                            for c in g.coords:
                                crossing_rows.append({
                                    "ward_id": ward_id,
                                    "ward_name_en": ward_name_en,
                                    "ward_name": ward_name,
                                    "geometry": Point(c)
                                })
                        else:
                            for ls in g.geoms:
                                for c in ls.coords:
                                    crossing_rows.append({
                                        "ward_id": ward_id,
                                        "ward_name_en": ward_name_en,
                                        "ward_name": ward_name,
                                        "geometry": Point(c)
                                    })
            elif inter.geom_type in ("LineString", "MultiLineString"):
                if inter.geom_type == "LineString":
                    for c in inter.coords:
                        crossing_rows.append({
                            "ward_id": ward_id,
                            "ward_name_en": ward_name_en,
                            "ward_name": ward_name,
                            "geometry": Point(c)
                        })
                else:
                    for ls in inter.geoms:
                        for c in ls.coords:
                            crossing_rows.append({
                                "ward_id": ward_id,
                                "ward_name_en": ward_name_en,
                                "ward_name": ward_name,
                                "geometry": Point(c)
                            })
    
    crossing_points = gpd.GeoDataFrame(crossing_rows, crs=METRIC_CRS)
    print(f"  Found {len(crossing_points)} crossing points", flush=True)
    
    # Step 5: Create 50m buffers
    print("Creating buffers...", flush=True)
    crossing_buffers = crossing_points.copy()
    crossing_buffers["geometry"] = crossing_buffers.geometry.buffer(50)
    
    # Step 6: Clip buffers to respective wards
    print("Clipping buffers to wards...", flush=True)
    ward_geom_lookup = dict(zip(wards_metric["ward_id"], wards_metric.geometry))
    
    clipped_geoms = []
    for idx, buf in crossing_buffers.iterrows():
        ward_geom = ward_geom_lookup.get(buf["ward_id"])
        if ward_geom is not None:
            clipped = buf.geometry.intersection(ward_geom)
            clipped_geoms.append(clipped)
        else:
            clipped_geoms.append(None)
    
    buffer_ward_intersection = crossing_buffers.copy()
    buffer_ward_intersection["geometry"] = clipped_geoms
    buffer_ward_intersection = buffer_ward_intersection[~buffer_ward_intersection.geometry.is_empty]
    buffer_ward_intersection = buffer_ward_intersection[buffer_ward_intersection.geometry.notna()]
    print(f"  {len(buffer_ward_intersection)} clipped buffer geometries", flush=True)
    
    # Step 7: Calculate crossing density
    print("Calculating ward crossing density...", flush=True)
    crossing_counts = crossing_points.groupby("ward_id").size().reset_index(name="crossing_count")
    
    wards_metric["ward_area_km2"] = wards_metric.geometry.area / 1_000_000
    
    density = wards.copy()
    density = density.merge(crossing_counts, on="ward_id", how="left")
    density = density.merge(wards_metric[["ward_id", "ward_area_km2"]], on="ward_id", how="left")
    density["crossing_count"] = density["crossing_count"].fillna(0).astype(int)
    density["crossings_per_km2"] = density["crossing_count"] / density["ward_area_km2"]
    density = density.to_crs(WGS84)
    
    # Step 8: Save to GPKG
    print("Saving to GPKG...", flush=True)
    
    wards_jgd = wards_metric.copy()
    
    # Ensure all required columns exist
    for gdf in [wards_jgd, crossing_points, crossing_buffers, buffer_ward_intersection]:
        for col in ["ward_id", "ward_name_en", "ward_name"]:
            if col not in gdf.columns:
                gdf[col] = None
    
    for col in ["ward_id", "ward_name_en", "ward_name"]:
        if col not in density.columns:
            density[col] = None
    
    wards_jgd.to_file(OUTPUT_GPKG, layer="wards_jgd", driver="GPKG")
    print("  Saved wards_jgd")
    
    crossing_points.to_file(OUTPUT_GPKG, layer="crossing_points", driver="GPKG")
    print("  Saved crossing_points")
    
    crossing_buffers.to_file(OUTPUT_GPKG, layer="crossing_buffers_50m", driver="GPKG")
    print("  Saved crossing_buffers_50m")
    
    buffer_ward_intersection.to_file(OUTPUT_GPKG, layer="buffer_ward_intersection", driver="GPKG")
    print("  Saved buffer_ward_intersection")
    
    density.to_file(OUTPUT_GPKG, layer="ward_crossing_density_wgs84", driver="GPKG")
    print("  Saved ward_crossing_density_wgs84")
    
    print("Done!", flush=True)


if __name__ == "__main__":
    main()
```

Wait, I just realized there might be issues with saving multiple layers to a single GPKG using GeoPandas. GeoPandas' `to_file` with the `driver="GPKG"` and `layer` parameter should work for appending to an existing GPKG. But the first call creates the file and subsequent calls append.

Actually, GeoPandas' `to_file` with `layer` parameter should work for creating/appending to a GPKG. Let me verify this is correct.

`gdf.to_file("file.gpkg", layer="layername", driver="GPKG")` should work. The first call creates the file, subsequent calls add layers.

OK, I think this script is mostly correct. Let me also think about edge cases:

1. What if a road crosses the boundary at a node (intersection)? In osmnx, road edges are between nodes, so the edge geometry starts and ends at nodes. If a node is exactly on the boundary, the edge might intersect the boundary at the node. This is actually a valid crossing.

2. What about roads that start at the boundary? If a road enters a ward from the boundary, the crossing point is at the boundary. This is correct.

3. What about roads that are entirely within the ward? They don't intersect the boundary, so they're correctly excluded.

4. What about the buffer? 50m buffer in the metric CRS should be 50 meters in real distance.

Let me also think about whether `simplify=False` is needed for osmnx. With simplify=False, we get all the original OSM nodes, which gives us more detailed geometry. This is important for accurate crossing detection.

Actually wait, `simplify=True` (default) simplifies the graph by removing interstitial nodes. This means a road edge might be a straight line between intersections, even if the original road curves. For crossing detection, this could miss some crossings or create inaccurate ones.

With `simplify=False`, we get all the original OSM nodes, so the road geometry follows the actual road more closely. This is better for crossing detection. But it means more edges and slower processing.

Let me use `simplify=False` for accuracy.

OK let me finalize and write the script.

Actually, wait. Let me reconsider the Overpass query for wards. The issue is that `rel(...) out geom;` on a boundary relation might not give us the polygon geometry directly. In Overpass JSON output, for a relation with `out geom`:

- If the relation is a multipolygon (type=multipolygon), the geometry field contains the polygon rings
- If the relation is a boundary (type=boundary), same thing hopefully

But the format of the geometry field for a relation is a flat list of lat/lon pairs, not separated into rings. For a simple boundary, this is the outer ring. For a boundary with multiple parts, it might be problematic.

Actually, looking at the Overpass JSON documentation more carefully, the `geometry` field for a relation with `out geom` is documented as follows:

> For a relation the geometry is defined as the convex hull of all its member coordinates, or as the multipolygon geometry if the relation is of type multipolygon or boundary.

Wait, I need to check this. Let me look at the actual Overpass behavior.

For a boundary relation with `out geom`, the geometry should be reconstructed as a multipolygon. But in the JSON output, it's represented as a flat array of coordinates.

Actually, I think the safest approach is to use Overpass to get the ways that make up the boundary, or use a different query approach.

Let me use this approach for wards:

```
[out:json][timeout:120];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
(
  rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
);
out body;
>;
out skel qt;
```

Then construct geometries from the ways. But this is complex.

Alternatively, I can use osmnx for the wards:

```python
wards = ox.geometries_from_place("Tokyo, Japan", tags={
    'boundary': 'administrative',
    'admin_level': '9'
})
```

Let me go with this approach. osmnx handles the geometry construction from OSM data.

Actually, I just realized there might be a simpler approach. Let me use `ox.geocode_to_gdf` with a specific query for each ward, or use `ox.geometries_from_bbox`.

Let me try this approach:

```python
# Get the bounding box of Tokyo's 23 wards
bbox = 35.55, 35.80, 139.55, 139.95  # rough bbox

wards = ox.geometries_from_bbox(35.80, 35.55, 139.95, 139.55, 
                                  tags={'boundary': 'administrative', 'admin_level': '9'})
```

This should get all admin_level=9 boundaries in the Tokyo area. Then I can filter to only keep the 23 special wards.

Let me go with this combined approach and write the final script.

OK, I'll write the script now. Let me be practical and just code it up.

A few more considerations:
1. The Overpass query might be large and slow due to the amount of data
2. The road network for 23 wards might have thousands of edges
3. The crossing detection might be slow without spatial indexing

I'll use spatial indexing in the form of geopandas operations.

Let me write the complete script now.

3:41

Write /work/solve.py

content · 319 lines · py

#!/usr/bin/env python3
"""
Tokyo 23 special wards road-crossing analysis.

Fetches ward boundaries and drivable roads (excluding service roads) from OSM,
finds every place a road crosses a ward boundary, creates 50 m buffers, clips
them to the ward, and computes crossing density per ward.

Output: /work/tokyo_crossings.gpkg with five layers
(wards_jgd, crossing_points, crossing_buffers_50m, buffer_ward_intersection,
 ward_crossing_density_wgs84).
"""

import sys
import time
import requests

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

# ---------------------------------------------------------------------------
#  Configuration
# ---------------------------------------------------------------------------
WGS84 = "EPSG:4326"
METRIC_CRS = "EPSG:6676"       # JGD2011 / UTM zone 54N – Tokyo metric grid
OUTPUT_GPKG = "/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",
]

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

ox.config(use_cache=True, log_console=True)


# ---------------------------------------------------------------------------
#  Overpass helpers
# ---------------------------------------------------------------------------
def overpass_query(query: str, max_retries: int = 3) -> dict:
    """Run an Overpass QL query, retrying on error with backoff."""
    for attempt in range(max_retries):
        for url in OVERPASS_URLS:
            try:
                resp = requests.post(url, data={"data": query}, timeout=300)
                if resp.status_code == 200:
                    return resp.json()
                print(f"  {url} → HTTP {resp.status_code}", file=sys.stderr, flush=True)
            except Exception as e:
                print(f"  {url}: {e}", file=sys.stderr, flush=True)
        if attempt < max_retries - 1:
            wait = 5 * (attempt + 1)
            print(f"  Retrying in {wait}s …", file=sys.stderr, flush=True)
            time.sleep(wait)
    raise RuntimeError("All Overpass query attempts failed")


# ---------------------------------------------------------------------------
#  Step 1 – Ward boundaries from Overpass
# ---------------------------------------------------------------------------
def fetch_wards() -> gpd.GeoDataFrame:
    """Fetch the 23 special-ward boundaries from Overpass.

    Returns a GeoDataFrame in WGS84 with columns
    ``ward_id``, ``ward_name_en``, ``ward_name``, ``geometry``.
    """
    print("Step 1 – Fetching ward boundaries from Overpass …", flush=True)

    # Get every admin_level=9 boundary relation inside Tokyo prefecture.
    # The special wards are tagged with admin_level=9 and boundary=administrative.
    query = """
    [out:json][timeout:120];
    area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
    (
      rel(area.tokyo)["admin_level"="9"]["boundary"="administrative"];
    );
    out geom;
    """
    data = overpass_query(query)

    rows = []
    for el in data.get("elements", []):
        if el.get("type") != "relation":
            continue

        tags = el.get("tags", {})
        ward_name_en = tags.get("name:en", "")
        ward_name = tags.get("name", "")

        # Only keep the 23 special wards
        if ward_name_en not in TOKYO_WARDS_EN:
            continue

        ward_id = str(el["id"])

        # The 'geometry' field from out geom on a boundary relation gives
        # the outer ring as a flat list of {lat,lon} dicts.
        raw_geom = el.get("geometry", [])
        if len(raw_geom) < 3:
            print(f"  Skip ward {ward_name_en} – insufficient geometry points",
                  file=sys.stderr, flush=True)
            continue

        coords = [(p["lon"], p["lat"]) for p in raw_geom]
        from shapely.geometry import Polygon as ShpPolygon
        poly = ShpPolygon(coords)
        if not poly.is_valid:
            poly = poly.buffer(0)

        rows.append({
            "ward_id": ward_id,
            "ward_name_en": ward_name_en,
            "ward_name": ward_name,
            "geometry": poly,
        })

    gdf = gpd.GeoDataFrame(rows, crs=WGS84)
    print(f"  → {len(gdf)} wards retrieved", flush=True)
    missing = set(TOKYO_WARDS_EN) - set(gdf["ward_name_en"])
    if missing:
        print(f"  ⚠ Missing: {missing}", flush=True)
    return gdf


# ---------------------------------------------------------------------------
#  Step 2 – Drivable road network (osmnx, excludes service roads)
# ---------------------------------------------------------------------------
def fetch_roads(polygon) -> gpd.GeoDataFrame:
    """Download drivable road edges inside *polygon* via osmnx.

    ``network_type='drive'`` excludes service roads, tracks, footways, etc.
    Returns a GeoDataFrame in WGS84.
    """
    print("Step 2 – Fetching drivable road network via osmnx …", flush=True)
    G = ox.graph_from_polygon(polygon, network_type="drive", simplify=False)
    edges = ox.graph_to_gdfs(G, nodes=False, edges=True)
    print(f"  → {len(edges)} road edges", flush=True)
    return edges


# ---------------------------------------------------------------------------
#  Step 3 – Find road-boundary crossing points
# ---------------------------------------------------------------------------
def find_crossings(wards_metric: gpd.GeoDataFrame,
                   edges_metric: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """For each ward find every place a road edge crosses its boundary.

    Returns a Point GeoDataFrame in METRIC_CRS.
    """
    print("Step 3 – Finding boundary crossings …", flush=True)

    rows = []
    for _, ward in wards_metric.iterrows():
        wid = ward["ward_id"]
        wen = ward["ward_name_en"]
        wn  = ward["ward_name"]
        boundary = ward.geometry.boundary

        # Edges that touch the boundary
        mask = edges_metric.intersects(boundary)
        if not mask.any():
            continue

        for _, road in edges_metric[mask].iterrows():
            inter = road.geometry.intersection(boundary)
            if inter.is_empty:
                continue

            it = inter.geom_type
            if it == "Point":
                rows.append({"ward_id": wid, "ward_name_en": wen,
                             "ward_name": wn, "geometry": inter})
            elif it == "MultiPoint":
                for pt in inter.geoms:
                    rows.append({"ward_id": wid, "ward_name_en": wen,
                                 "ward_name": wn, "geometry": pt})
            elif it == "GeometryCollection":
                for g in inter.geoms:
                    if g.geom_type == "Point":
                        rows.append({"ward_id": wid, "ward_name_en": wen,
                                     "ward_name": wn, "geometry": g})
                    elif g.geom_type in ("LineString",):
                        # road hugs boundary – count endpoints
                        for c in g.coords:
                            rows.append({"ward_id": wid, "ward_name_en": wen,
                                         "ward_name": wn,
                                         "geometry": Point(c)})
            elif it in ("LineString",):
                for c in inter.coords:
                    rows.append({"ward_id": wid, "ward_name_en": wen,
                                 "ward_name": wn, "geometry": Point(c)})

    result = gpd.GeoDataFrame(rows, crs=METRIC_CRS)
    print(f"  → {len(result)} crossing points", flush=True)
    return result


# ---------------------------------------------------------------------------
#  Step 4 – 50 m buffers
# ---------------------------------------------------------------------------
def create_buffers(crossings: gpd.GeoDataFrame,
                   radius: float = 50.0) -> gpd.GeoDataFrame:
    print(f"Step 4 – Creating {radius:.0f} m buffers …", flush=True)
    buf = crossings.copy()
    buf["geometry"] = buf.geometry.buffer(radius)
    return buf


# ---------------------------------------------------------------------------
#  Step 5 – Clip buffers to their ward
# ---------------------------------------------------------------------------
def clip_buffers_to_ward(buffers: gpd.GeoDataFrame,
                         wards_metric: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    print("Step 5 – Clipping buffers to ward polygon …", flush=True)
    ward_lookup = dict(zip(wards_metric["ward_id"], wards_metric.geometry))

    clipped = []
    for _, row in buffers.iterrows():
        ward_geom = ward_lookup.get(row["ward_id"])
        if ward_geom is not None:
            clipped.append(row.geometry.intersection(ward_geom))
        else:
            clipped.append(None)

    result = buffers.copy()
    result["geometry"] = clipped
    result = result[~result.geometry.is_empty & result.geometry.notna()].copy()
    print(f"  → {len(result)} clipped buffer geometries", flush=True)
    return result


# ---------------------------------------------------------------------------
#  Step 6 – Crossing-density table (WGS84 output)
# ---------------------------------------------------------------------------
def compute_density(wards_wgs84: gpd.GeoDataFrame,
                    wards_metric: gpd.GeoDataFrame,
                    crossings: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    print("Step 6 – Computing crossing densities …", flush=True)

    counts = crossings.groupby("ward_id").size().reset_index(name="crossing_count")

    # area from metric CRS
    wards_metric = wards_metric.copy()
    wards_metric["ward_area_km2"] = wards_metric.geometry.area / 1_000_000

    density = wards_wgs84[["ward_id", "ward_name_en", "ward_name", "geometry"]].copy()
    density = density.merge(counts, on="ward_id", how="left")
    density = density.merge(
        wards_metric[["ward_id", "ward_area_km2"]], on="ward_id", how="left"
    )
    density["crossing_count"] = density["crossing_count"].fillna(0).astype(int)
    density["crossings_per_km2"] = density["crossing_count"] / density["ward_area_km2"]
    return density


# ---------------------------------------------------------------------------
#  Main
# ---------------------------------------------------------------------------
def main():
    # ---- Fetch data ----
    wards_wgs84 = fetch_wards()
    if wards_wgs84.empty:
        print("ERROR: no ward boundaries fetched – cannot continue.", flush=True)
        sys.exit(1)

    # Build a unified polygon of the 23 wards to clip the road query.
    unified_poly = unary_union(wards_wgs84.geometry.values)
    # Buffer inward slightly to avoid fetching data just outside the wards
    unified_poly = unified_poly.buffer(-0.001)  # ~100 m inward at this latitude
    if unified_poly.is_empty:
        unified_poly = unary_union(wards_wgs84.geometry.values)

    roads_wgs84 = fetch_roads(unified_poly)

    # ---- Reproject to metric ----
    wards_metric = wards_wgs84.to_crs(METRIC_CRS)
    roads_metric = roads_wgs84.to_crs(METRIC_CRS)

    # ---- Analysis ----
    crossings = find_crossings(wards_metric, roads_metric)
    buffers = create_buffers(crossings, 50.0)
    clipped_buffers = clip_buffers_to_ward(buffers, wards_metric)
    density = compute_density(wards_wgs84, wards_metric, crossings)

    # ---- Write GPKG ----
    print("Writing /work/tokyo_crossings.gpkg …", flush=True)

    # Ensure all ward-specific layers carry the three id/name columns.
    for gdf_name, gdf in [
        ("wards_jgd",             wards_metric),
        ("crossing_points",       crossings),
        ("crossing_buffers_50m",  buffers),
        ("buffer_ward_intersection", clipped_buffers),
    ]:
        for col in ("ward_id", "ward_name_en", "ward_name"):
            if col not in gdf.columns:
                gdf[col] = None
        gdf.to_file(OUTPUT_GPKG, layer=gdf_name, driver="GPKG")
        print(f"  ✓ {gdf_name} ({len(gdf)} rows)")

    density.to_file(OUTPUT_GPKG, layer="ward_crossing_density_wgs84", driver="GPKG")
    print(f"  ✓ ward_crossing_density_wgs84 ({len(density)} rows)")

    print("All done!", flush=True)


if __name__ == "__main__":
    main()

Wrote /work/solve.py (12478 chars).