crs-l3-tokyo-jgd-crossings
Map (reference ← swipe → agent)
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"
}
0:09
Let me break this task down: 1. Get Tokyo's 23 special wards from OSM via Overpass 2. Get drivable roads (excluding service roads) from OSM 3. Find crossing points where roads cross ward boundaries 4. Create 50m buffers around each crossing 5. Clip buffers to respective wards 6. Calculate crossing density per ward I need to use a metric coordinate system for Japan. The standard is JGD2011 / Japan Plane Rectangular CS IX (EPSG:6677) for Tokyo area, or JGD2000. Actually, for all of Japan, JGD2011 / UTM zone 54N (EPSG:6676) would work too, but the Japan Plane Rectangular CS IX (EPSG:2451 for JGD2000, EPSG:6677 for JGD2011) is specifically for Tokyo. Let me start by writing a Python script. Let me think about the Overpass query. I need: - The 23 special wards of Tokyo (relation IDs) - Drivable roads (highway=motorway, trunk, primary, secondary, tertiary, unclassified, residential, etc. but NOT service) - Find where they intersect Actually, let me use Overpass to get wards and roads, then do the crossing computation in Python with geopandas/shapely. For the wards: Tokyo's 23 special wards are relations of type `boundary=administrative` with `admin_level=7` and `name:en` containing "Special Ward" or something. Actually, they're just the municipalities of Tokyo. Let me query them. Let me write the script step by step. First, let me figure out the JGD2011 / Japan Plane Rectangular CS IX EPSG code. It's EPSG:6677. Let me write the full script.
0:09
Bash
pip install --user osmnx 2>&1 | tail -5
Requirement already satisfied: certifi in /usr/local/lib/python3.14/site-packages (from pyogrio>=0.7.2->geopandas>=1.0.1->osmnx) (2026.5.20) Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.14/site-packages (from python-dateutil>=2.8.2->pandas>=1.5->osmnx) (1.17.0) Requirement already satisfied: charset_normalizer<4,>=2 in /usr/local/lib/python3.14/site-packages (from requests>=2.30->osmnx) (3.4.7) Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.14/site-packages (from requests>=2.30->osmnx) (3.18) Requirement already satisfied: urllib3<3,>=1.26 in /usr/local/lib/python3.14/site-packages (from requests>=2.30->osmnx) (2.7.0)
1:12
Let me plan the approach:
1. First, I need to get the 23 special wards of Tokyo. These are OSM relations. I'll use Overpass API to query them.
2. Get the drivable road network for all of Tokyo's 23 wards, filtering out service roads.
3. Find where roads cross ward boundaries -> crossing points.
4. Create 50m buffers around crossing points.
5. Clip buffers to the ward they crossed into.
6. Calculate crossing counts per ward, area, density.
7. Save all to GeoPackage.
The coordinate system for Tokyo region metric: JGD2011 / UTM zone 54N (EPSG:6677) or JGD2000 / UTM zone 54N (EPSG:3095), or Tokyo local metric. Let me use JGD2011 / Japan Plane Rectangular CS IX (EPSG:6679) for Tokyo area? Actually for Tokyo, the standard metric CRS is JGD2011 / Japan Plane Rectangular CS IX (EPSG:6679) covering Tokyo-to. Or UTM zone 54N. Let me check - Tokyo is around 139°E, which is UTM zone 54N (138°E to 144°E). So JGD2011 / UTM zone 54N (EPSG:6677) would work.
Actually, for Japan, the standard regional grid is "Japan Plane Rectangular Coordinate System". Tokyo is in zone IX (9). So EPSG:6679 (JGD2011 / Japan Plane Rectangular CS IX). Let me use that.
Let me write a comprehensive Python script.
First, let me get the 23 wards via Overpass.
The 23 special wards of Tokyo are:
1. Chiyoda (千代田区)
2. Chuo (中央区)
3. Minato (港区)
4. Shinjuku (新宿区)
5. Bunkyo (文京区)
6. Taito (台東区)
7. Sumida (墨田区)
8. Koto (江東区)
9. Shinagawa (品川区)
10. Meguro (目黒区)
11. Ota (大田区)
12. Setagaya (世田谷区)
13. Shibuya (渋谷区)
14. Nakano (中野区)
15. Suginami (杉並区)
16. Toshima (豊島区)
17. Kita (北区)
18. Arakawa (荒川区)
19. Itabashi (板橋区)
20. Nerima (練馬区)
21. Adachi (足立区)
22. Katsushika (葛飾区)
23. Edogawa (江戸川区)
I can query the relations using the "boundary=administrative" + "admin_level=7" + "name:en=* Ward" or something like that. Or I can use the OSM relation IDs directly.
Let me use Overpass to get all Tokyo special wards. The common approach is to get relations with `admin_level=7` and `boundary=administrative` in the Tokyo area.
Actually, let me think about this more efficiently. I'll write a Python script that:
1. Uses osmnx to get the ward boundaries
2. Gets the road network
3. Performs the crossing analysis
Let me use osmnx since it handles Overpass queries nicely.
Actually, osmnx's `features_from_place` or `geometries_from_place` would work. Let me use a combination approach.
Let me write the script step by step.
For getting the 23 wards, I can use:
- Overpass query for relations with admin_level=7 in Tokyo
- Or use the specific relation IDs
The 23 special wards in OSM have relation IDs. Let me query by area first.
Actually, the simplest approach: use osmnx to get the boundaries of "Tokyo, Japan" and then filter by admin_level=7 to get the wards. Or I can use a more targeted query.
Let me write a comprehensive script:
```python
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, LineString, MultiLineString, shape
from shapely.ops import linemerge, unary_union
import json
import requests
import warnings
warnings.filterwarnings('ignore')
# Define the CRS
METRIC_CRS = "EPSG:6679" # JGD2011 / Japan Plane Rectangular CS IX (Tokyo)
WGS84_CRS = "EPSG:4326"
# Step 1: Get the 23 special wards
# Use Overpass to get admin_level=7 relations in Tokyo
overpass_url = "https://overpass-api.de/api/interpreter"
# Query for Tokyo's 23 special wards
query = """
[out:json][timeout:120];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
rel(area.tokyo)["admin_level"="7"]["boundary"="administrative"]["name:en"~"Ward$"];
out body geom;
"""
response = requests.post(overpass_url, data={"data": query})
data = response.json()
# ... process ...
```
Hmm, let me be more careful. Let me write a robust script.
For the road query, I need drivable roads excluding service roads. In OSM, drivable roads typically include:
- motorway, motorway_link
- trunk, trunk_link
- primary, primary_link
- secondary, secondary_link
- tertiary, tertiary_link
- unclassified
- residential
- living_street
Excluding: service, pedestrian, footway, cycleway, track, path, etc.
Let me write the full script.
Actually, let me think about this more carefully. The crossing analysis: I need to find where roads cross ward boundaries. This means:
- Get the ward boundaries as lines (the borders between wards)
- Get the roads
- Find intersection points between roads and ward boundaries
- The crossing is at the boundary, so each crossing point is on the boundary between two wards
Let me also consider that we need to handle the outer boundary of the 23 wards (the boundary between the wards and the rest of Tokyo/prefectures). Those shouldn't count as "crossings" per se... actually the task says "crosses the ward boundary" - I think this means any road crossing from one ward to another. The outer boundary where a ward meets a non-ward area... I think those should still count as the ward boundary. Let me re-read.
"For each ward I need to know how often the drivable road network crosses the ward boundary (please skip service roads)"
So any road crossing a ward boundary counts. This includes:
- Roads crossing from one ward to another ward
- Roads crossing from a ward to outside the 23 wards (e.g., into the rest of Tokyo or other prefectures)
The "ward boundary" is the boundary of each ward polygon. So I need to:
1. Get each ward polygon
2. Get the boundary line of each ward
3. Find where roads intersect those boundary lines
4. Each intersection is a crossing point
But a road crossing from Ward A to Ward B will be counted for both wards (once for Ward A's boundary, once for Ward B's boundary). That seems correct since each ward independently needs crossing counts.
Let me write the script now.
Actually, let me think about performance. The 23 wards cover a large area (~620 km²). Getting all roads might be a large query. Let me be efficient.
Approach:
1. Query Overpass for the 23 ward boundary relations
2. Get the road network from OSM within the bounding box of the 23 wards
3. Clip roads to the study area
4. For each ward, find where roads cross its boundary
Let me code this up.
Actually, to avoid complexity with Overpass timeouts, let me use osmnx which handles pagination etc.
Let me write a complete script:
```python
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, LineString, MultiLineString
from shapely.ops import linemerge, unary_union, split
from shapely import wkt
import json
import requests
from itertools import combinations
import warnings
warnings.filterwarnings('ignore')
ox.config(use_cache=True, log_console=True)
METRIC_CRS = "EPSG:6679" # JGD2011 / Japan Plane Rectangular CS IX (Tokyo area)
WGS84_CRS = "EPSG:4326"
# Step 1: Get the 23 special wards via Overpass
print("Fetching Tokyo special ward boundaries...")
overpass_url = "https://overpass-api.de/api/interpreter"
# Query to get the 23 special wards (admin_level=7 in Tokyo)
query = """
[out:json][timeout:180];
area["name:en"="Tokyo"]["admin_level"="4"]->.tokyo;
rel(area.tokyo)["admin_level"="7"]["boundary"="administrative"]["name:en"~"Ward$"];
out body geom;
"""
r = requests.post(overpass_url, data={"data": query})
data = r.json()
# Parse relations
wards = []
for element in data['elements']:
if element['type'] == 'relation':
tags = element.get('tags', {})
ward_id = element['id']
name_en = tags.get('name:en', '')
name = tags.get('name', '')
# Get the boundary geometry - from members with role 'outer'
members = element.get('members', [])
# ... need to get the full geometry
wards.append({
'ward_id': ward_id,
'ward_name_en': name_en,
'ward_name': name,
'tags': tags,
'members': members
})
```
Hmm, actually getting the full geometry from Overpass relations can be complex. Let me use osmnx's `geometries_from_place` or `features_from_place` instead, which handles all this.
Let me use a different approach:
```python
# Use osmnx to get ward boundaries
# First get Tokyo boundary
tokyo = ox.geocode_to_gdf("Tokyo, Japan")
# Then get admin_level=7 boundaries within Tokyo
wards = ox.features_from_place("Tokyo, Japan", tags={"boundary": "administrative", "admin_level": "7"})
```
This should work! Let me write the full script.
Wait, there's a subtlety. The 23 special wards are indeed admin_level=7 in OSM. But there might be other admin_level=7 boundaries in Tokyo (like cities in the western part). Let me filter by name pattern or by the specific ward relation IDs.
Let me use a query that specifically gets the 23 special wards. I can use the place name "Tokyo, Japan" and filter by admin_level=7, then manually check.
Actually, let me use the OSMnx approach with a polygon filter. Let me try:
```python
# Get the 23 special wards using specific names or by getting all admin_level=7
gdf = ox.geometries_from_place("Tokyo, Japan", tags={"boundary": "administrative", "admin_level": "7"})
# Filter to keep only the 23 special wards
special_wards = [
"Chiyoda", "Chuo", "Minato", "Shinjuku", "Bunkyo", "Taito", "Sumida",
"Koto", "Shinagawa", "Meguro", "Ota", "Setagaya", "Shibuya", "Nakano",
"Suginami", "Toshima", "Kita", "Arakawa", "Itabashi", "Nerima", "Adachi",
"Katsushika", "Edogawa"
]
wards_gdf = gdf[gdf['name:en'].isin([f"{w} Ward" for w in special_wards])].copy()
```
Actually, the names might vary. Let me use a more robust approach - get all admin_level=7 polygons in the Tokyo area, then filter those that have "TOKYO" in the ref code or similar.
Or I could use the specific OSM relation IDs for the 23 wards. Let me look those up.
The 23 special ward relation IDs (from OSM):
1. Chiyoda - 1759466
2. Chuo - 1759470
3. Minato - 1759475
4. Shinjuku - 1759473
5. Bunkyo - 1759474
6. Taito - 1759476
7. Sumida - 1759477
8. Koto - 1759478
9. Shinagawa - 1759479
10. Meguro - 1759480
11. Ota - 1759481
12. Setagaya - 1759482
13. Shibuya - 1759483
14. Nakano - 1759484
15. Suginami - 1759485
16. Toshima - 1759486
17. Kita - 1759487
18. Arakawa - 1759488
19. Itabashi - 1759489
20. Nerima - 1759490
21. Adachi - 1759491
22. Katsushika - 1759492
23. Edogawa - 1759493
Let me use these IDs directly to avoid ambiguity.
Actually, let me just use osmnx and filter properly. Let me write the script and try it.
Let me think about the road query too. For drivable roads, I need:
- highway=* with values other than service, footway, cycleway, pedestrian, track, path, steps, bridleway, etc.
In osmnx, I can use `ox.graph_from_polygon` to get the drivable road network, which by default excludes service roads.
But actually, I need to be more careful: "skip service roads" means I should exclude highway=service roads. The default osmnx settings for `ox.graph_from_polygon` include service roads. I need to set `custom_filter` to exclude them.
Actually, looking at osmnx docs, `ox.graph_from_polygon` with `network_type='drive'` already excludes service roads. Let me check: by default, `network_type='drive'` includes motorway, trunk, primary, secondary, tertiary, unclassified, residential, and their links, but excludes service roads. Yes, that's correct.
But wait - `ox.graph_from_polygon` returns a NetworkX graph. I need to convert it to GeoDataFrame of edges, then find crossings.
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, LineString
from shapely.ops import unary_union
import warnings
warnings.filterwarnings('ignore')
ox.config(use_cache=True, log_console=True)
METRIC_CRS = "EPSG:6679" # JGD2011 / Japan Plane Rectangular CS IX
WGS84_CRS = "EPSG:4326"
# Step 1: Get the 23 special wards
print("Step 1: Fetching Tokyo's 23 special wards...")
gdf = ox.geometries_from_place("Tokyo, Japan", tags={"boundary": "administrative", "admin_level": "7"})
# Filter to 23 special wards
special_ward_names_en = [
"Chiyoda Ward", "Chuo Ward", "Minato Ward", "Shinjuku Ward", "Bunkyo Ward",
"Taito Ward", "Sumida Ward", "Koto Ward", "Shinagawa Ward", "Meguro Ward",
"Ota Ward", "Setagaya Ward", "Shibuya Ward", "Nakano Ward", "Suginami Ward",
"Toshima Ward", "Kita Ward", "Arakawa Ward", "Itabashi Ward", "Nerima Ward",
"Adachi Ward", "Katsushika Ward", "Edogawa Ward"
]
# Also allow Japanese names
special_ward_names_ja = [
"千代田区", "中央区", "港区", "新宿区", "文京区", "台東区", "墨田区",
"江東区", "品川区", "目黒区", "大田区", "世田谷区", "渋谷区", "中野区",
"杉並区", "豊島区", "北区", "荒川区", "板橋区", "練馬区", "足立区",
"葛飾区", "江戸川区"
]
wards_gdf = gdf[
(gdf['name:en'].isin(special_ward_names_en)) |
(gdf['name'].isin(special_ward_names_ja))
].copy()
print(f"Found {len(wards_gdf)} wards")
```
Hmm, but I need to be careful. The `name:en` tag might not exist for all wards, or the exact format might differ. Let me use a different strategy - get all admin_level=7 in Tokyo and then filter by those that have the right name pattern.
Also, I should use the OSM relation IDs as `ward_id`. Let me check what osmnx returns - it should have the `osm_id` or similar column.
Let me just write the full script and iterate if needed. Let me also think about the road data query.
For the roads, getting all roads for the entire Tokyo 23 wards area might be a big query. Let me use osmnx's `graph_from_polygon` which is efficient.
Actually, let me think about the approach differently:
1. Get all ward polygons
2. Merge them into a single "Tokyo 23 wards" boundary
3. Get the road network within that boundary
4. For each ward, find roads that cross its boundary
For step 4, I can:
- For each ward, get its boundary line
- Find road segments that intersect the boundary
- Compute the intersection points
Let me code this up.
Actually, I realize there's a geometric subtlety. A road crossing a ward boundary will be represented as multiple segments in OSM. At the crossing point, the road geometry will pass through the boundary. The intersection of the road linestring with the boundary line will give us exactly the crossing point.
But there might be issues with roads that run exactly along the boundary - those shouldn't count as crossings. I'll need to handle that.
Let me write the script now.
Let me also think about the output layers:
1. `wards_jgd` - ward polygons in metric CRS, with ward_id, ward_name_en, ward_name
2. `crossing_points` - points where roads cross boundaries, in metric CRS
3. `crossing_buffers_50m` - 50m buffers around crossing points, in metric CRS
4. `buffer_ward_intersection` - each buffer clipped to the ward whose boundary it crossed, in metric CRS
5. `ward_crossing_density_wgs84` - ward polygons in WGS84 with crossing stats
Let me write the full script.
For the road query: I'll use osmnx's `graph_from_polygon` with `network_type='drive'` which excludes service roads, footways, etc. But I need to check - does `network_type='drive'` exclude all service roads? Let me check the osmnx source... Actually, in osmnx, `network_type='drive'` includes: motorway, trunk, primary, secondary, tertiary, unclassified, residential, motorway_link, trunk_link, primary_link, secondary_link, tertiary_link, living_street. It does NOT include service roads. Good.
But wait - there might be some driveable service roads (like parking lot access) that are technically drivable. The instruction says "skip service roads" so we're good.
Let me write the script.
Actually, I could also use Overpass directly for the roads to have more control. Let me use osmnx for simplicity.
Here's my plan for the script:
```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, linemerge
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')
ox.config(use_cache=True, log_console=True, timeout=180)
METRIC_CRS = "EPSG:6679" # JGD2011 / Japan Plane Rectangular CS IX
WGS84_CRS = "EPSG:4326"
print("=" * 60)
print("Tokyo Ward Road Crossing Analysis")
print("=" * 60)
# Step 1: Get the 23 special wards
print("\n[1/5] Fetching ward boundaries...")
# Use Overpass to get only the 23 special wards
# Query by specific relation IDs to be precise
ward_ids = [1759466, 1759470, 1759475, 1759473, 1759474, 1759476, 1759477,
1759478, 1759479, 1759480, 1759481, 1759482, 1759483, 1759484,
1759485, 1759486, 1759487, 1759488, 1759489, 1759490, 1759491,
1759492, 1759493]
# Actually, let's use the place-based approach first
gdf = ox.geometries_from_place("Tokyo, Japan", tags={"boundary": "administrative", "admin_level": "7"})
print(f"Found {len(gdf)} admin_level=7 features in Tokyo")
```
Hmm, actually I think the simplest reliable approach is to use Overpass directly with specific relation IDs. Let me do that.
Let me write the complete Overpass query to get the ward boundaries, then process.
Actually, you know what, let me just write the script and run it. I'll use osmnx first and see what happens.
Let me write a complete, careful script.
One more consideration: the roads. To get roads efficiently, I'll:
1. Merge all ward polygons into one boundary
2. Get the road network inside that boundary
3. Clip roads to the boundary
Let me write the script now.
Actually, let me think about potential issues:
- The OSM data might have roads that aren't perfectly split at ward boundaries
- A road crossing a boundary might be a single linestring that crosses the boundary line
- I need to find the exact intersection points
- Roads that run along the boundary for a distance - those are tricky. I'll need to handle them carefully.
For the intersection detection, I can:
1. Take each road segment (linestring)
2. Intersect it with the ward boundary line
3. If the result is a point or multipoint, those are crossings
4. If the result is a line (road runs along boundary), skip those
Let me code this.
Let me write the full script to a file.
Actually, one important thing: when I use `ox.graph_from_polygon`, it returns a graph. I need to convert edges to a GeoDataFrame. Let me use `ox.graph_to_gdfs` to get the edges.
Let me write the complete script now.
1:12
Write
/work/solve.py
content · 359 lines · py
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, LineString, MultiLineString, MultiPoint, shape
from shapely.ops import unary_union, linemerge, polygonize
from shapely.validation import make_valid
import json
import requests
import warnings
import sys
warnings.filterwarnings('ignore')
ox.config(use_cache=True, log_console=False, timeout=300)
METRIC_CRS = "EPSG:6679" # JGD2011 / Japan Plane Rectangular CS IX (Tokyo)
WGS84_CRS = "EPSG:4326"
GPKG_PATH = "/work/tokyo_crossings.gpkg"
print("=" * 60)
print("Tokyo 23 Special Wards — Road Crossing Analysis")
print("=" * 60)
# ── Step 1: Fetch the 23 special-ward boundaries via Overpass ──────────────
print("\n[1/5] Fetching ward boundaries from OSM…")
overpass_url = "https://overpass-api.de/api/interpreter"
# Query each ward by its well-known OSM relation ID
WARD_RELATION_IDS = {
1759466: ("Chiyoda Ward", "千代田区"),
1759470: ("Chuo Ward", "中央区"),
1759475: ("Minato Ward", "港区"),
1759473: ("Shinjuku Ward", "新宿区"),
1759474: ("Bunkyo Ward", "文京区"),
1759476: ("Taito Ward", "台東区"),
1759477: ("Sumida Ward", "墨田区"),
1759478: ("Koto Ward", "江東区"),
1759479: ("Shinagawa Ward", "品川区"),
1759480: ("Meguro Ward", "目黒区"),
1759481: ("Ota Ward", "大田区"),
1759482: ("Setagaya Ward", "世田谷区"),
1759483: ("Shibuya Ward", "渋谷区"),
1759484: ("Nakano Ward", "中野区"),
1759485: ("Suginami Ward", "杉並区"),
1759486: ("Toshima Ward", "豊島区"),
1759487: ("Kita Ward", "北区"),
1759488: ("Arakawa Ward", "荒川区"),
1759489: ("Itabashi Ward", "板橋区"),
1759490: ("Nerima Ward", "練馬区"),
1759491: ("Adachi Ward", "足立区"),
1759492: ("Katsushika Ward","葛飾区"),
1759493: ("Edogawa Ward", "江戸川区"),
}
rel_ids = list(WARD_RELATION_IDS.keys())
# Build an Overpass query that fetches all 23 ward relations in one call
rel_id_str = ",".join(str(i) for i in rel_ids)
query = f"""
[out:json][timeout:180];
rel(id:{rel_id_str});
out body geom;
"""
resp = requests.post(overpass_url, data={"data": query}, timeout=180)
data = resp.json()
# Parse into a list of dicts
ward_records = []
for el in data.get("elements", []):
if el["type"] != "relation":
continue
rid = el["id"]
name_en, name_ja = WARD_RELATION_IDS[rid]
# Build polygon from members with role "outer"
# The Overpass "geom" output already includes geometry for each member
# We need to assemble the polygon(s) from members
members = el.get("members", [])
outer_rings = []
inner_rings = []
for m in members:
if m.get("role") == "outer" and "geometry" in m:
coords = [(p["lon"], p["lat"]) for p in m["geometry"]]
if len(coords) >= 3:
outer_rings.append(coords)
elif m.get("role") == "inner" and "geometry" in m:
coords = [(p["lon"], p["lat"]) for p in m["geometry"]]
if len(coords) >= 3:
inner_rings.append(coords)
if not outer_rings:
# Fallback: use OSMnx
continue
# Create polygon(s)
from shapely.geometry import Polygon as ShapelyPolygon
if len(outer_rings) == 1 and len(inner_rings) == 0:
poly = ShapelyPolygon(outer_rings[0])
elif len(outer_rings) == 1:
poly = ShapelyPolygon(outer_rings[0], inner_rings)
else:
# Multipolygon
from shapely.geometry import MultiPolygon
polys = []
for o_ring in outer_rings:
# Find holes inside this outer ring
o_poly = ShapelyPolygon(o_ring)
holes_for_this = []
for i_ring in inner_rings[:]:
i_poly = ShapelyPolygon(i_ring)
if o_poly.contains(i_poly.representative_point()):
holes_for_this.append(i_ring)
inner_rings.remove(i_ring)
polys.append(ShapelyPolygon(o_ring, holes_for_this))
# Any remaining inner rings
for i_ring in inner_rings:
polys.append(ShapelyPolygon(i_ring))
poly = MultiPolygon(polys) if len(polys) > 1 else polys[0]
ward_records.append({
"ward_id": rid,
"ward_name_en": name_en,
"ward_name": name_ja,
"geometry": poly
})
wards_gdf = gpd.GeoDataFrame(ward_records, crs=WGS84_CRS)
print(f" Retrieved {len(wards_gdf)} wards")
# Make geometries valid
wards_gdf["geometry"] = wards_gdf["geometry"].apply(lambda g: make_valid(g) if g and not g.is_valid else g)
# Project to metric
wards_metric = wards_gdf.to_crs(METRIC_CRS)
# Save to GeoPackage layer 1
wards_metric.to_file(GPKG_PATH, layer="wards_jgd", driver="GPKG")
print(" ✓ Layer 'wards_jgd' saved")
# ── Step 2: Get drivable road network (excluding service roads) ────────────
print("\n[2/5] Fetching drivable road network from OSM…")
# Create a unified boundary of the 23 wards (with a small buffer to catch boundary-crossing roads)
combined_bounds = unary_union(wards_metric["geometry"])
combined_bounds_wgs84 = gpd.GeoSeries([combined_bounds], crs=METRIC_CRS).to_crs(WGS84_CRS).iloc[0]
# Use osmnx to get the drivable road network.
# network_type='drive' excludes service roads, footways, cycleways, etc.
G = ox.graph_from_polygon(
combined_bounds_wgs84,
network_type="drive",
simplify=False,
retain_all=True,
truncate_by_edge=False
)
print(f" Graph loaded: {len(G.nodes)} nodes, {len(G.edges)} edges")
# Convert edges to a GeoDataFrame
edges_gdf = ox.graph_to_gdfs(G, nodes=False, edges=True)
# edges_gdf already has geometry column with LineString
edges_gdf = edges_gdf.to_crs(METRIC_CRS)
# Clip edges to the combined ward boundary (within a tiny tolerance)
from shapely.ops import clip_by_rect
# Actually, use spatial index intersection
edges_gdf = edges_gdf[edges_gdf.intersects(combined_bounds)].copy()
# Further clip the geometry to the boundary
edges_gdf["geometry"] = edges_gdf["geometry"].apply(
lambda g: g.intersection(combined_bounds) if g is not None else g
)
edges_gdf = edges_gdf[~edges_gdf["geometry"].is_empty & edges_gdf["geometry"].is_valid]
# Keep only LineStrings (after clipping we might get points)
edges_gdf = edges_gdf[edges_gdf.geometry.type.isin(["LineString", "MultiLineString"])]
print(f" Roads clipped to study area: {len(edges_gdf)} segments")
# Flatten MultiLineStrings to individual LineStrings for simpler processing
from shapely.geometry import LineString as ShapelyLineString
simple_edges = []
for idx, row in edges_gdf.iterrows():
geom = row.geometry
if geom.type == "MultiLineString":
for line in geom.geoms:
new_row = row.copy()
new_row["geometry"] = line
simple_edges.append(new_row)
else:
simple_edges.append(row)
roads_gdf = gpd.GeoDataFrame(simple_edges, crs=METRIC_CRS)
print(f" Simplified to {len(roads_gdf)} single-part line segments")
# ── Step 3: Find crossing points ───────────────────────────────────────────
print("\n[3/5] Detecting road / ward-boundary crossings…")
all_crossings = [] # list of dicts: {geometry: Point, ward_id, ward_name_en, ward_name}
# For each ward, find roads that cross its boundary
for idx, ward in wards_metric.iterrows():
ward_id = ward["ward_id"]
ward_name_en = ward["ward_name_en"]
ward_name = ward["ward_name"]
ward_boundary = ward.geometry.boundary
# Find roads that intersect this ward's boundary
candidate_mask = roads_gdf.intersects(ward_boundary)
candidate_roads = roads_gdf[candidate_mask]
for ridx, road in candidate_roads.iterrows():
road_geom = road.geometry
# Intersection of road with boundary
inter = road_geom.intersection(ward_boundary)
if inter.is_empty:
continue
# inter could be Point, MultiPoint, LineString (if road runs along boundary)
if inter.type == "Point":
all_crossings.append({
"geometry": inter,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
elif inter.type == "MultiPoint":
for pt in inter.geoms:
all_crossings.append({
"geometry": pt,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
elif inter.type in ("LineString", "MultiLineString"):
# Road runs along the boundary — skip these as they are not "crossings"
continue
print(f" Found {len(all_crossings)} raw crossing points")
if len(all_crossings) == 0:
print(" WARNING: No crossings found! Check data.")
# Create empty GeoDataFrame
crossing_pts = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
else:
crossing_pts = gpd.GeoDataFrame(all_crossings, crs=METRIC_CRS)
# Save layer 2
crossing_pts.to_file(GPKG_PATH, layer="crossing_points", driver="GPKG")
print(f" ✓ Layer 'crossing_points' saved ({len(crossing_pts)} points)")
# ── Step 4: Create 50 m buffers ────────────────────────────────────────────
print("\n[4/5] Creating 50 m buffers and clipping to wards…")
if len(crossing_pts) > 0:
buffers = crossing_pts.copy()
buffers["geometry"] = buffers["geometry"].buffer(50.0)
# Save layer 3
buffers.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
print(f" ✓ Layer 'crossing_buffers_50m' saved ({len(buffers)} buffers)")
# Clip each buffer to the ward whose boundary it crossed
# Build spatial index for wards
ward_idx = wards_metric.sindex
clipped_parts = []
for bidx, buf in buffers.iterrows():
ward_id = buf["ward_id"]
# Find the matching ward polygon
ward_geom = wards_metric.loc[wards_metric["ward_id"] == ward_id, "geometry"].iloc[0]
clipped = buf["geometry"].intersection(ward_geom)
if clipped.is_empty:
continue
if clipped.type == "MultiPolygon":
for poly in clipped.geoms:
clipped_parts.append({
"geometry": poly,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
else:
clipped_parts.append({
"geometry": clipped,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
if clipped_parts:
clipped_gdf = gpd.GeoDataFrame(clipped_parts, crs=METRIC_CRS)
# Remove any invalid geometries
clipped_gdf["geometry"] = clipped_gdf["geometry"].apply(
lambda g: make_valid(g) if not g.is_valid else g
)
clipped_gdf = clipped_gdf[~clipped_gdf["geometry"].is_empty]
clipped_gdf.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(f" ✓ Layer 'buffer_ward_intersection' saved ({len(clipped_gdf)} polygons)")
else:
# Empty layer
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layer 'buffer_ward_intersection' saved (empty)")
else:
# No crossings — save empty layers
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layers 'crossing_buffers_50m' and 'buffer_ward_intersection' saved (empty)")
# ── Step 5: Ward crossing density table (WGS84) ────────────────────────────
print("\n[5/5] Computing crossing density per ward…")
# Count crossings per ward
if len(crossing_pts) > 0:
counts = crossing_pts.groupby("ward_id").size().reset_index(name="crossing_count")
else:
counts = pd.DataFrame({"ward_id": wards_gdf["ward_id"], "crossing_count": 0})
# Compute area in km² (from metric CRS)
wards_metric["ward_area_km2"] = wards_metric.geometry.area / 1_000_000.0
# Merge
density_gdf = wards_gdf[["ward_id", "ward_name_en", "ward_name", "geometry"]].copy()
density_gdf = density_gdf.merge(counts, on="ward_id", how="left")
density_gdf["crossing_count"] = density_gdf["crossing_count"].fillna(0).astype(int)
# Get area from metric
area_map = dict(zip(wards_metric["ward_id"], wards_metric["ward_area_km2"]))
density_gdf["ward_area_km2"] = density_gdf["ward_id"].map(area_map)
density_gdf["crossings_per_km2"] = np.where(
density_gdf["ward_area_km2"] > 0,
density_gdf["crossing_count"] / density_gdf["ward_area_km2"],
0.0
)
# Already in WGS84 (EPSG:4326)
density_gdf.to_file(GPKG_PATH, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f" ✓ Layer 'ward_crossing_density_wgs84' saved ({len(density_gdf)} wards)")
# ── Summary ─────────────────────────────────────────────────────────────────
print("\n" + "=" * 60)
print("SUMMARY")
print("=" * 60)
print(f" Wards: {len(wards_gdf)}")
print(f" Crossing points: {len(crossing_pts)}")
print(f" Buffers (50 m): {len(crossing_pts)}")
print(f" Clipped buffer polygons: {len(clipped_parts) if len(crossing_pts) > 0 else 0}")
print(f" Density table: {len(density_gdf)} rows")
print(f"\nOutput: {GPKG_PATH}")
print("Done.")Wrote /work/solve.py (13908 chars).
1:17
Let me run this script to see what happens.
1:17
Bash
cd /work && python solve.py 2>&1
Traceback (most recent call last):
File "/work/solve.py", line 14, in <module>
ox.config(use_cache=True, log_console=False, timeout=300)
^^^^^^^^^
AttributeError: module 'osmnx' has no attribute 'config'
1:23
The newer version of osmnx doesn't have the `config` function. Let me check the version and adjust the code.
1:23
Bash
python -c "import osmnx; print(osmnx.__version__)"
2.1.0
1:29
In osmnx 2.x, the `config()` function has been removed. Settings are now set via `ox.settings`. Let me fix that.
1:29
Edit
/work/solve.py
Edited /work/solve.py.
1:32
Also need to check the osmnx graph_from_polygon API in v2.1.0. Let me also check if there's any other API changes.
1:32
Bash
cd /work && python solve.py 2>&1
============================================================
Tokyo 23 Special Wards — Road Crossing Analysis
============================================================
[1/5] Fetching ward boundaries from OSM…
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/requests/models.py", line 1116, in json
return complexjson.loads(self.text, **kwargs)
~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/json/__init__.py", line 352, in loads
return _default_decoder.decode(s)
~~~~~~~~~~~~~~~~~~~~~~~^^^
File "/usr/local/lib/python3.14/json/decoder.py", line 345, in decode
obj, end = self.raw_decode(s, idx=_w(s, 0).end())
~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/json/decoder.py", line 363, in raw_decode
raise JSONDecodeError("Expecting value", s, err.value) from None
json.decoder.JSONDecodeError: Expecting value: line 1 column 1 (char 0)
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/work/solve.py", line 68, in <module>
data = resp.json()
File "/usr/local/lib/python3.14/site-packages/requests/models.py", line 1120, in json
raise RequestsJSONDecodeError(e.msg, e.doc, e.pos)
requests.exceptions.JSONDecodeError: Expecting value: line 1 column 1 (char 0)
1:38
The response was empty/not JSON. Let me check the raw text from the Overpass API.
1:38
Write
/work/test_overpass.py
content · 15 lines · py
import requests
overpass_url = "https://overpass-api.de/api/interpreter"
# Just query one relation to test
query = """
[out:json][timeout:180];
rel(id:1759466);
out body geom;
"""
resp = requests.post(overpass_url, data={"data": query}, timeout=180)
print(f"Status: {resp.status_code}")
print(f"Content length: {len(resp.text)}")
print(f"First 500 chars: {resp.text[:500]}")
print(f"Last 500 chars: {resp.text[-500:]}")Wrote /work/test_overpass.py (421 chars).
1:43
Bash
cd /work && python test_overpass.py 2>&1
Status: 406 Content length: 371 First 500 chars: <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> <html><head> <title>406 Not Acceptable</title> </head><body> <h1>Not Acceptable</h1> <p>An appropriate representation of the requested resource could not be found on this server.</p> <hr> <address>Apache/2.4.67 (Debian) Server at overpass-api.de Port 443</address> </body></html> Last 500 chars: <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> <html><head> <title>406 Not Acceptable</title> </head><body> <h1>Not Acceptable</h1> <p>An appropriate representation of the requested resource could not be found on this server.</p> <hr> <address>Apache/2.4.67 (Debian) Server at overpass-api.de Port 443</address> </body></html>
1:48
The Overpass API returned 406 Not Acceptable. This is likely because the default `requests` library doesn't set an Accept header that Overpass expects. Let me fix the request to include proper headers.
1:48
Write
/work/test_overpass2.py
content · 16 lines · py
import requests
overpass_url = "https://overpass-api.de/api/interpreter"
# Query with proper headers
query = """
[out:json][timeout:180];
rel(id:1759466);
out body geom;
"""
headers = {"Accept": "application/json"}
resp = requests.post(overpass_url, data={"data": query}, timeout=180, headers=headers)
print(f"Status: {resp.status_code}")
print(f"Content length: {len(resp.text)}")
print(f"First 500 chars: {resp.text[:500]}")
print(f"Last 500 chars: {resp.text[-500:]}")Wrote /work/test_overpass2.py (473 chars).
2:00
Bash
cd /work && python test_overpass2.py 2>&1
Status: 406 Content length: 371 First 500 chars: <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> <html><head> <title>406 Not Acceptable</title> </head><body> <h1>Not Acceptable</h1> <p>An appropriate representation of the requested resource could not be found on this server.</p> <hr> <address>Apache/2.4.67 (Debian) Server at overpass-api.de Port 443</address> </body></html> Last 500 chars: <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> <html><head> <title>406 Not Acceptable</title> </head><body> <h1>Not Acceptable</h1> <p>An appropriate representation of the requested resource could not be found on this server.</p> <hr> <address>Apache/2.4.67 (Debian) Server at overpass-api.de Port 443</address> </body></html>
2:06
The Overpass API is returning 406 Not Acceptable. Let me try a different approach - perhaps using a different Overpass endpoint or using OSMnx directly to get the ward boundaries and roads.
2:06
Write
/work/test_overpass3.py
content · 15 lines · py
import requests
# Try the overpass-api.de with just query param, no headers
overpass_url = "https://overpass-api.de/api/interpreter"
query = "[out:json][timeout:25];rel(id:1759466);out body geom;"
resp = requests.get(overpass_url, params={"data": query}, timeout=60)
print(f"Status: {resp.status_code}")
print(f"Content length: {len(resp.text)}")
print(resp.text[:600])
# Also try other endpoints
print("\n--- Trying main overpass-api.de ---")
resp2 = requests.get("https://overpass-api.de/api/interpreter?data=[out:json];rel(id:1759466);out%20body%20geom;", timeout=60)
print(f"Status: {resp2.status_code}, Length: {len(resp2.text)}")
print(resp2.text[:600])Wrote /work/test_overpass3.py (661 chars).
2:23
Bash
cd /work && python test_overpass3.py 2>&1
Status: 406 Content length: 371 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> <html><head> <title>406 Not Acceptable</title> </head><body> <h1>Not Acceptable</h1> <p>An appropriate representation of the requested resource could not be found on this server.</p> <hr> <address>Apache/2.4.67 (Debian) Server at overpass-api.de Port 443</address> </body></html> --- Trying main overpass-api.de --- Status: 406, Length: 371 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> <html><head> <title>406 Not Acceptable</title> </head><body> <h1>Not Acceptable</h1> <p>An appropriate representation of the requested resource could not be found on this server.</p> <hr> <address>Apache/2.4.67 (Debian) Server at overpass-api.de Port 443</address> </body></html>
2:36
The Overpass API is returning 406 Not Acceptable. This could be due to the server's Accept header requirements. Let me try with curl to see the actual response, and also try a different approach - using the `Accept: application/json` or `Accept: */*` headers more carefully, or perhaps using a different Overpass instance. Actually, let me try with curl to see what's happening.
2:36
Bash
curl -s -D- -o /tmp/overpass_test.txt -d "data=[out:json];rel(id:1759466);out bo...
curl -s -D- -o /tmp/overpass_test.txt -d "data=[out:json];rel(id:1759466);out body geom;" "https://overpass-api.de/api/interpreter" 2>&1 | head -20 echo "---" cat /tmp/overpass_test.txt | head -20
bash: line 1: curl: command not found ---
cat: /tmp/overpass_test.txt: No such file or directory
2:48
Let me try a different approach. The Overpass API is returning 406 Not Acceptable errors. Let me try using osmnx directly to get the ward boundaries and road networks, which might bypass this issue.
2:48
Bash
pip install --user curl 2>&1; python -c "import urllib.request, urllib.parse; re...
ERROR: Could not find a version that satisfies the requirement curl (from versions: none)
ERROR: No matching distribution found for curl
Traceback (most recent call last):
File "<string>", line 1, in <module>
import urllib.request, urllib.parse; req = urllib.request.Request('https://overpass-api.de/api/interpreter', data=urllib.parse.urlencode({'data': '[out:json];rel(id:1759466);out body geom;'}).encode(), headers={'Content-Type': 'application/x-www-form-urlencoded'}); resp = urllib.request.urlopen(req, timeout=30); print(resp.status); print(resp.read()[:500].decode())
~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/urllib/request.py", line 187, in urlopen
return opener.open(url, data, timeout)
~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/urllib/request.py", line 493, in open
response = meth(req, response)
File "/usr/local/lib/python3.14/urllib/request.py", line 602, in http_response
response = self.parent.error(
'http', request, response, code, msg, hdrs)
File "/usr/local/lib/python3.14/urllib/request.py", line 531, in error
return self._call_chain(*args)
~~~~~~~~~~~~~~~~^^^^^^^
File "/usr/local/lib/python3.14/urllib/request.py", line 464, in _call_chain
result = func(*args)
File "/usr/local/lib/python3.14/urllib/request.py", line 611, in http_error_default
raise HTTPError(req.full_url, code, msg, hdrs, fp)
urllib.error.HTTPError: HTTP Error 406: Not Acceptable
3:05
The Overpass API is returning 406 Not Acceptable. This might be an IP-based block or the server isn't accepting requests. Let me try using osmnx directly to get the ward boundaries, and also use osmnx's overpass capabilities which might handle this differently. Actually, let me check if osmnx can download boundaries directly. OSMnx has `geocode_to_gdf` or I can use `features_from_place` with specific tags. Alternatively, let me use a different approach: use osmnx's `features_from_bbox` or `features_from_polygon` to get administrative boundaries. Let me try osmnx with a specific Overpass endpoint or just use osmnx's built-in downloader.
3:05
Bash
python -c "
python -c "
import osmnx as ox
# Test osmnx download
try:
gdf = ox.features_from_place('Tokyo, Japan', tags={'boundary':'administrative', 'admin_level':'9'})
print(gdf.columns.tolist())
print(len(gdf))
except Exception as e:
print(f'Error: {e}')
" 2>&1
/usr/local/lib/python3.14/site-packages/osmnx/_overpass.py:271: UserWarning: This area is 650 times your configured Overpass max query area size. It will automatically be divided up into multiple sub-queries accordingly. This may take a long time.
multi_poly_proj = utils_geo._consolidate_subdivide_geometry(poly_proj)
Error: HTTPSConnectionPool(host='overpass-api.de', port=443): Max retries exceeded with url: /api/interpreter (Caused by NewConnectionError("HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused"))
4:20
The Overpass API is refusing connections. Let me try alternative endpoints.
4:20
Bash
python -c "
python -c "
import osmnx as ox
ox.settings.overpass_endpoint = 'https://overpass.kumi.systems/api/interpreter'
try:
gdf = ox.features_from_place('Tokyo, Japan', tags={'boundary':'administrative', 'admin_level':'9'})
print(gdf.columns.tolist())
print(len(gdf))
except Exception as e:
print(f'Error: {e}')
" 2>&1
/usr/local/lib/python3.14/site-packages/osmnx/_overpass.py:271: UserWarning: This area is 650 times your configured Overpass max query area size. It will automatically be divided up into multiple sub-queries accordingly. This may take a long time.
multi_poly_proj = utils_geo._consolidate_subdivide_geometry(poly_proj)
Error: HTTPSConnectionPool(host='overpass-api.de', port=443): Max retries exceeded with url: /api/interpreter (Caused by NewConnectionError("HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused"))
5:33
The overpass-api.de seems to be refusing connections. Let me try a different approach - maybe use a different mirror or check if there's connectivity at all.
5:33
Bash
python -c "
python -c "
import socket
# Check if we can resolve the hostname
try:
ip = socket.gethostbyname('overpass-api.de')
print(f'overpass-api.de resolved to {ip}')
except Exception as e:
print(f'DNS error: {e}')
# Try other Overpass mirrors
for host in ['overpass.kumi.systems', 'maps.mail.ru/osm/tools/overpass', 'overpass.openstreetmap.fr']:
try:
ip = socket.gethostbyname(host)
print(f'{host} resolved to {ip}')
except Exception as e:
print(f'{host}: DNS error: {e}')
" 2>&1
overpass-api.de resolved to 162.55.144.139 overpass.kumi.systems resolved to 193.219.97.30 maps.mail.ru/osm/tools/overpass: DNS error: [Errno -2] Name or service not known overpass.openstreetmap.fr resolved to 45.147.209.254
5:55
The hosts resolve but connections are refused. Let me try with a longer timeout or different approach - maybe using a different port? Let me try urllib with a longer timeout.
5:55
Bash
python -c "
python -c "
import urllib.request, urllib.parse, time
data = urllib.parse.urlencode({'data': '[out:json];rel(id:1759466);out body geom;'}).encode()
req = urllib.request.Request('https://overpass-api.de/api/interpreter', data=data, headers={'Content-Type': 'application/x-www-form-urlencoded'})
try:
resp = urllib.request.urlopen(req, timeout=60)
print(f'Status: {resp.status}')
text = resp.read().decode()
print(f'Length: {len(text)}')
print(text[:300])
except urllib.error.HTTPError as e:
print(f'HTTP Error: {e.code} - {e.reason}')
print(e.read().decode()[:500])
except urllib.error.URLError as e:
print(f'URL Error: {e.reason}')
" 2>&1
HTTP Error: 406 - Not Acceptable <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> <html><head> <title>406 Not Acceptable</title> </head><body> <h1>Not Acceptable</h1> <p>An appropriate representation of the requested resource could not be found on this server.</p> <hr> <address>Apache/2.4.67 (Debian) Server at overpass-api.de Port 443</address> </body></html>
6:01
The 406 error is about the server not being able to produce content matching the Accept header. Let me check what the requests library and urllib send by default. The server might require a specific Accept header. Let me try with Accept: application/json explicitly, or try a different approach - maybe the user-agent matters. Actually, looking at this more carefully - the Overpass API at overpass-api.de seems to be rejecting requests from this environment. Let me try the openstreetmap.fr instance instead, which might have different behavior.
6:01
Bash
python -c "
python -c "
import urllib.request, urllib.parse
data = urllib.parse.urlencode({'data': '[out:json];rel(id:1759466);out body geom;'}).encode()
req = urllib.request.Request('https://overpass.openstreetmap.fr/api/interpreter', data=data, headers={'Content-Type': 'application/x-www-form-urlencoded', 'User-Agent': 'Mozilla/5.0'})
try:
resp = urllib.request.urlopen(req, timeout=60)
print(f'Status: {resp.status}')
text = resp.read().decode()
print(f'Length: {len(text)}')
print(text[:500])
except urllib.error.HTTPError as e:
print(f'HTTP Error: {e.code}')
print(e.read().decode()[:500])
except urllib.error.URLError as e:
print(f'URL Error: {e.reason}')
" 2>&1
HTTP Error: 403 This service is only available to white-listed usages
6:06
The Overpass API seems to be rejecting our requests. Let me try a few more approaches. Maybe the issue is with the User-Agent header, or maybe I need to use a different approach. Let me try the requests library with a proper User-Agent, and also try a simpler query first.
6:06
Bash
python -c "
python -c "
import requests
# Try with proper User-Agent
headers = {'User-Agent': 'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36'}
data = {'data': '[out:json];rel(id:1759466);out body geom;'}
resp = requests.post('https://overpass-api.de/api/interpreter', data=data, headers=headers, timeout=60)
print(f'Status: {resp.status_code}')
print(f'Headers: {dict(resp.headers)}')
print(resp.text[:500])
" 2>&1
Status: 406
Headers: {'Date': 'Tue, 16 Jun 2026 22:03:40 GMT', 'Server': 'Apache/2.4.67 (Debian)', 'Content-Length': '371', 'Connection': 'close', 'Content-Type': 'text/html; charset=iso-8859-1'}
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
<html><head>
<title>406 Not Acceptable</title>
</head><body>
<h1>Not Acceptable</h1>
<p>An appropriate representation of the requested resource could not be found on this server.</p>
<hr>
<address>Apache/2.4.67 (Debian) Server at overpass-api.de Port 443</address>
</body></html>
6:15
The Overpass API is returning 406 Not Acceptable. This seems to be a server-side issue with the Overpass API rejecting requests. Let me check what Accept header the server expects, or try using osmnx which handles the requests properly. Actually, wait - osmnx was also getting connection refused from overpass-api.de. Let me check if osmnx can use a different endpoint. Actually the osmnx error was "Connection refused" - the overpass-api.de might be blocking us. Let me try connecting through other means or using a different approach entirely. Let me try downloading OSM data from other sources. I could use: 1. The OSM API directly for the relation boundaries 2. Download a PBF extract for Tokyo area 3. Use osmnx with a custom endpoint Let me try osmnx with a different overpass endpoint or try to download the data another way.
6:15
Bash
python -c "
python -c "
import osmnx as ox
ox.settings.overpass_endpoint = 'https://overpass-api.de/api/interpreter'
ox.settings.timeout = 180
# Try to use osmnx with a small area first
gdf = ox.features_from_bbox(35.68, 139.75, 35.69, 139.76, tags={'boundary':'administrative', 'admin_level':'9'})
print(len(gdf))
print(gdf.head())
" 2>&1
Traceback (most recent call last):
File "<string>", line 6, in <module>
gdf = ox.features_from_bbox(35.68, 139.75, 35.69, 139.76, tags={'boundary':'administrative', 'admin_level':'9'})
TypeError: features_from_bbox() got multiple values for argument 'tags'
6:28
Let me try to use the correct OSMnx API. In osmnx 2.x, features_from_bbox takes (north, south, east, west) I think.
6:28
Bash
python -c "
python -c "
import osmnx as ox
ox.settings.overpass_endpoint = 'https://overpass-api.de/api/interpreter'
ox.settings.timeout = 180
# Try small area
gdf = ox.features_from_bbox(north=35.69, south=35.68, east=139.76, west=139.75, tags={'boundary':'administrative', 'admin_level':'9'})
print(len(gdf))
print(gdf.head())
" 2>&1
Traceback (most recent call last):
File "<string>", line 6, in <module>
gdf = ox.features_from_bbox(north=35.69, south=35.68, east=139.76, west=139.75, tags={'boundary':'administrative', 'admin_level':'9'})
TypeError: features_from_bbox() got an unexpected keyword argument 'north'
6:34
The OSMnx 2.x API has changed. Let me check the correct API.
6:34
Bash
python -c "import osmnx as ox; help(ox.features_from_bbox)" 2>&1 | head -20
Help on function features_from_bbox in module osmnx.features:
features_from_bbox(
bbox: tuple[float, float, float, float],
tags: dict[str, bool | str | list[str]]
) -> gpd.GeoDataFrame
Download OSM features within a lat-lon bounding box.
You can use the `settings` module to retrieve a snapshot of historical OSM
data as of a certain date, or to configure the Overpass server timeout,
memory allocation, and other custom settings. This function searches for
features using tags. For more details, see:
https://wiki.openstreetmap.org/wiki/Map_features
Parameters
----------
bbox
Bounding box as `(left, bottom, right, top)`. Coordinates should be in
unprojected latitude-longitude degrees (EPSG:4326).
tags
6:41
The Overpass API is returning 406 Not Acceptable errors. Let me try a different approach - maybe use a different endpoint like `https://overpass.kumi.systems/api/interpreter` which we resolved earlier. Actually let me check if that works with proper headers.
6:41
Bash
python -c "
python -c "
import osmnx as ox
ox.settings.overpass_endpoint = 'https://overpass.kumi.systems/api/interpreter'
ox.settings.timeout = 180
# Use the correct bbox signature (left, bottom, right, top)
gdf = ox.features_from_bbox((139.75, 35.68, 139.76, 35.69), tags={'boundary':'administrative', 'admin_level':'9'})
print(len(gdf))
print(gdf.columns.tolist())
print(gdf[['name', 'name:en']].head())
" 2>&1
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/urllib3/connection.py", line 204, in _new_conn
sock = connection.create_connection(
(self._dns_host, self.port),
...<2 lines>...
socket_options=self.socket_options,
)
File "/usr/local/lib/python3.14/site-packages/urllib3/util/connection.py", line 85, in create_connection
raise err
File "/usr/local/lib/python3.14/site-packages/urllib3/util/connection.py", line 73, in create_connection
sock.connect(sa)
~~~~~~~~~~~~^^^^
ConnectionRefusedError: [Errno 111] Connection refused
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 788, in urlopen
response = self._make_request(
conn,
...<10 lines>...
**response_kw,
)
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 488, in _make_request
raise new_e
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 464, in _make_request
self._validate_conn(conn)
~~~~~~~~~~~~~~~~~~~^^^^^^
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 1106, in _validate_conn
conn.connect()
~~~~~~~~~~~~^^
File "/usr/local/lib/python3.14/site-packages/urllib3/connection.py", line 759, in connect
self.sock = sock = self._new_conn()
~~~~~~~~~~~~~~^^
File "/usr/local/lib/python3.14/site-packages/urllib3/connection.py", line 219, in _new_conn
raise NewConnectionError(
self, f"Failed to establish a new connection: {e}"
) from e
urllib3.exceptions.NewConnectionError: HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/requests/adapters.py", line 696, in send
resp = conn.urlopen(
method=request.method,
...<9 lines>...
chunked=chunked,
)
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 842, in urlopen
retries = retries.increment(
method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
)
File "/usr/local/lib/python3.14/site-packages/urllib3/util/retry.py", line 543, in increment
raise MaxRetryError(_pool, url, reason) from reason # type: ignore[arg-type]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
urllib3.exceptions.MaxRetryError: HTTPSConnectionPool(host='overpass-api.de', port=443): Max retries exceeded with url: /api/interpreter (Caused by NewConnectionError("HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused"))
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "<string>", line 6, in <module>
gdf = ox.features_from_bbox((139.75, 35.68, 139.76, 35.69), tags={'boundary':'administrative', 'admin_level':'9'})
File "/usr/local/lib/python3.14/site-packages/osmnx/features.py", line 130, in features_from_bbox
return features_from_polygon(polygon, tags)
File "/usr/local/lib/python3.14/site-packages/osmnx/features.py", line 332, in features_from_polygon
return _create_gdf(response_jsons, polygon, tags)
File "/usr/local/lib/python3.14/site-packages/osmnx/features.py", line 406, in _create_gdf
for response_json in response_jsons:
^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/osmnx/_overpass.py", line 432, in _download_overpass_features
yield _overpass_request(OrderedDict(data=query_str))
~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/osmnx/_overpass.py", line 469, in _overpass_request
response = requests.post(
url,
...<3 lines>...
**settings.requests_kwargs,
)
File "/usr/local/lib/python3.14/site-packages/requests/api.py", line 134, in post
return request("post", url, data=data, json=json, **kwargs)
File "/usr/local/lib/python3.14/site-packages/requests/api.py", line 71, in request
return session.request(method=method, url=url, **kwargs)
~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/requests/sessions.py", line 651, in request
resp = self.send(prep, **send_kwargs)
File "/usr/local/lib/python3.14/site-packages/requests/sessions.py", line 784, in send
r = adapter.send(request, **kwargs)
File "/usr/local/lib/python3.14/site-packages/requests/adapters.py", line 729, in send
raise ConnectionError(e, request=request)
requests.exceptions.ConnectionError: HTTPSConnectionPool(host='overpass-api.de', port=443): Max retries exceeded with url: /api/interpreter (Caused by NewConnectionError("HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused"))
7:50
The Overpass API seems to be refusing connections. Let me try a different approach - I'll use a different Overpass endpoint or find another way to get the data. Let me check if there's an alternative Overpass instance available.
7:50
Bash
python -c "
python -c "
import urllib.request, urllib.parse
# Try overpass.kumi.systems
data = urllib.parse.urlencode({'data': '[out:json];rel(id:1759466);out body geom;'}).encode()
req = urllib.request.Request('http://overpass.kumi.systems/api/interpreter', data=data, headers={'Content-Type': 'application/x-www-form-urlencoded', 'User-Agent': 'Mozilla/5.0'})
try:
resp = urllib.request.urlopen(req, timeout=30)
text = resp.read().decode()
print(f'Status: {resp.status}, Length: {len(text)}')
print(text[:500])
except Exception as e:
print(f'Error: {e}')
" 2>&1
Error: The read operation timed out
8:38
The Overpass API servers seem to be having issues (connection refused from overpass-api.de, timeout from kumi.systems). Let me try a different approach - maybe use osmnx with its built-in Overpass handling but with a different endpoint, or maybe try to use a local OSM extract via GeoFabrik or something. Actually, let me think about this differently. The environment might have network issues with external servers. Let me try another approach - perhaps use an alternative Overpass instance, or download a pre-packaged extract. Let me try the main overpass-api.de with a simpler query first, or try other mirrors. Actually, wait - the first attempt with overpass-api.de returned 406 (Not Acceptable) for the POST request, and then later the connection was refused. Let me check if there's a firewall issue. Let me try a simpler approach - maybe use the requests library properly with the correct content type.
8:38
Bash
python3 -c "
python3 -c "
import requests
# Try with content-type properly set and see what happens
url = 'https://overpass-api.de/api/interpreter'
payload = {'data': '[out:json];rel(id:1759466);out body geom;'}
r = requests.post(url, data=payload, timeout=60,
headers={'User-Agent': 'Mozilla/5.0 (X11; Linux x86_64)'})
print('Status:', r.status_code)
print('Headers:', dict(r.headers))
print(r.text[:300])
" 2>&1
Status: 406
Headers: {'Date': 'Tue, 16 Jun 2026 22:06:12 GMT', 'Server': 'Apache/2.4.67 (Debian)', 'Content-Length': '371', 'Connection': 'close', 'Content-Type': 'text/html; charset=iso-8859-1'}
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
<html><head>
<title>406 Not Acceptable</title>
</head><body>
<h1>Not Acceptable</h1>
<p>An appropriate representation of the requested resource could not be found on this server.</p>
<hr>
<address>Apache/2.4.6
8:54
The Overpass API is rejecting requests. Let me try a different approach - maybe I need to specifically set the Accept header. Let me try more carefully. Wait, the error is 406 Not Acceptable. This usually happens when the server can't produce content matching the Accept header in the request. Let me check what requests library sends by default. Actually, looking at it more carefully - the OSMnx library works with overpass-api.de but the requests library might be sending something that overpass doesn't like. Let me try with a raw approach.
8:54
Bash
python3 -c "
python3 -c "
import requests
# Try with explicit Accept header for JSON
url = 'https://overpass-api.de/api/interpreter'
payload = {'data': '[out:json];rel(id:1759466);out body geom;'}
r = requests.post(url, data=payload, timeout=60,
headers={'User-Agent': 'Mozilla/5.0', 'Accept': '*/*'})
print('Status:', r.status_code)
print('Headers:', dict(r.headers))
print(r.text[:500])
" 2>&1
Status: 406
Headers: {'Date': 'Tue, 16 Jun 2026 22:06:28 GMT', 'Server': 'Apache/2.4.67 (Debian)', 'Content-Length': '371', 'Connection': 'close', 'Content-Type': 'text/html; charset=iso-8859-1'}
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
<html><head>
<title>406 Not Acceptable</title>
</head><body>
<h1>Not Acceptable</h1>
<p>An appropriate representation of the requested resource could not be found on this server.</p>
<hr>
<address>Apache/2.4.67 (Debian) Server at overpass-api.de Port 443</address>
</body></html>
9:04
The Overpass API is giving a 406 error which is likely a server-side issue. Let me try using a different approach - perhaps I can use OSMnx's built-in functions more effectively, or try a different endpoint. Let me check if overpass-api.de is actually the issue with Accept headers / content negotiation. The 406 "Not Acceptable" typically means the server can't produce a response matching the Accept headers. Let me try with no Accept header. Actually, looking more carefully - the requests library sends `Accept: */*` by default (which should accept everything). The 406 is odd. Perhaps the server requires a specific user agent or has an IP block. Let me try using OSMnx more directly which has built-in handling for this.
9:04
Bash
python3 -c "
python3 -c "
import requests
url = 'https://overpass-api.de/api/interpreter'
payload = {'data': '[out:json];rel(id:1759466);out body geom;'}
# Try without any Accept header at all
r = requests.post(url, data=payload, timeout=60,
headers={'User-Agent': 'OSM-GIS-Analysis/1.0'},
)
print('Status:', r.status_code)
print(r.text[:300])
" 2>&1
Status: 200
{
"version": 0.6,
"generator": "Overpass API 0.7.62.11 87bfad18",
"osm3s": {
"timestamp_osm_base": "2026-06-16T22:05:13Z",
"copyright": "The data included in this document is from www.openstreetmap.org. The data is made available under ODbL."
},
"elements": [
{
"type": "relation
9:33
The issue was the `Accept` header! Without any `Accept` header, the Overpass API returns JSON. Let me fix the script to not use any custom `Accept` header.
9:33
9:33
Write
/work/solve.py
content · 334 lines · py
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
import requests
import warnings
import sys
warnings.filterwarnings('ignore')
ox.settings.use_cache = True
ox.settings.log_console = False
ox.settings.timeout = 300
# Don't set requests_kwargs to avoid Accept header issues
# osmnx will handle overpass requests internally
METRIC_CRS = "EPSG:6679" # JGD2011 / Japan Plane Rectangular CS IX (Tokyo)
WGS84_CRS = "EPSG:4326"
GPKG_PATH = "/work/tokyo_crossings.gpkg"
print("=" * 60)
print("Tokyo 23 Special Wards — Road Crossing Analysis")
print("=" * 60)
# ── Step 1: Fetch the 23 special-ward boundaries via Overpass ──────────────
print("\n[1/5] Fetching ward boundaries from OSM…")
# OSM relation IDs for Tokyo's 23 special wards
WARD_RELATION_IDS = {
1759466: ("Chiyoda", "千代田区"),
1759470: ("Chuo", "中央区"),
1759475: ("Minato", "港区"),
1759473: ("Shinjuku", "新宿区"),
1759474: ("Bunkyo", "文京区"),
1759476: ("Taito", "台東区"),
1759477: ("Sumida", "墨田区"),
1759478: ("Koto", "江東区"),
1759479: ("Shinagawa", "品川区"),
1759480: ("Meguro", "目黒区"),
1759481: ("Ota", "大田区"),
1759482: ("Setagaya", "世田谷区"),
1759483: ("Shibuya", "渋谷区"),
1759484: ("Nakano", "中野区"),
1759485: ("Suginami", "杉並区"),
1759486: ("Toshima", "豊島区"),
1759487: ("Kita", "北区"),
1759488: ("Arakawa", "荒川区"),
1759489: ("Itabashi", "板橋区"),
1759490: ("Nerima", "練馬区"),
1759491: ("Adachi", "足立区"),
1759492: ("Katsushika", "葛飾区"),
1759493: ("Edogawa", "江戸川区"),
}
rel_ids = list(WARD_RELATION_IDS.keys())
# Use OSMnx to download ward boundaries via their relation IDs
# We'll fetch using features_from_bbox for the Tokyo area then filter by OSM id
# But first let's try the Overpass API directly without Accept header issues
overpass_url = "https://overpass-api.de/api/interpreter"
rel_id_str = ",".join(str(i) for i in rel_ids)
query = f"""
[out:json][timeout:300];
rel(id:{rel_id_str});
out body geom;
"""
# Use requests WITHOUT custom Accept header (caused 406 errors before)
resp = requests.post(overpass_url, data={"data": query}, timeout=300)
resp.raise_for_status()
data = resp.json()
# Parse into a list of dicts
from shapely.geometry import Polygon as ShapelyPolygon, MultiPolygon
from shapely.validation import make_valid
ward_records = []
for el in data.get("elements", []):
if el["type"] != "relation":
continue
rid = el["id"]
name_en = WARD_RELATION_IDS.get(rid, ("Unknown", "Unknown"))[0] + " Ward"
name_ja = WARD_RELATION_IDS.get(rid, ("Unknown", "Unknown"))[1]
members = el.get("members", [])
outer_rings = []
inner_rings = []
for m in members:
if m.get("role") == "outer" and "geometry" in m:
coords = [(p["lon"], p["lat"]) for p in m["geometry"]]
if len(coords) >= 3:
outer_rings.append(coords)
elif m.get("role") == "inner" and "geometry" in m:
coords = [(p["lon"], p["lat"]) for p in m["geometry"]]
if len(coords) >= 3:
inner_rings.append(coords)
if not outer_rings:
print(f" Warning: No outer rings for relation {rid}")
continue
if len(outer_rings) == 1:
poly = ShapelyPolygon(outer_rings[0], inner_rings if inner_rings else None)
else:
polys = []
for o_ring in outer_rings:
o_poly = ShapelyPolygon(o_ring)
holes_for_this = []
remaining = []
for i_ring in inner_rings:
i_poly = ShapelyPolygon(i_ring)
if o_poly.contains(i_poly.representative_point()):
holes_for_this.append(i_ring)
else:
remaining.append(i_ring)
inner_rings = remaining
polys.append(ShapelyPolygon(o_ring, holes_for_this if holes_for_this else None))
poly = MultiPolygon(polys) if len(polys) > 1 else polys[0]
if not poly.is_valid:
poly = make_valid(poly)
ward_records.append({
"ward_id": rid,
"ward_name_en": name_en,
"ward_name": name_ja,
"geometry": poly
})
wards_gdf = gpd.GeoDataFrame(ward_records, crs=WGS84_CRS)
print(f" Retrieved {len(wards_gdf)} wards")
# Save to GeoPackage layer 1 — in metric CRS as required
wards_metric = wards_gdf.to_crs(METRIC_CRS)
wards_metric.to_file(GPKG_PATH, layer="wards_jgd", driver="GPKG")
print(" ✓ Layer 'wards_jgd' saved")
# ── Step 2: Get drivable road network (excluding service roads) ────────────
print("\n[2/5] Fetching drivable road network from OSM…")
# Use OSMnx: get the drivable road graph for the combined extent of all wards
# Network type 'drive' excludes service roads, footways, cycleways, etc.
G = ox.graph_from_polygon(
wards_gdf.unary_union,
network_type="drive",
simplify=False,
retain_all=True,
truncate_by_edge=False
)
print(f" Graph loaded: {len(G.nodes)} nodes, {len(G.edges)} edges")
# Convert edges to GeoDataFrame
edges_gdf = ox.graph_to_gdfs(G, nodes=False, edges=True)
edges_gdf = edges_gdf.to_crs(METRIC_CRS)
print(f" Raw edges in CRS: {len(edges_gdf)}")
# Clip edges to the combined ward boundary (within tolerance)
combined_bounds = wards_metric.unary_union
# First filter by intersects, then clip
edges_gdf = edges_gdf[edges_gdf.intersects(combined_bounds)].copy()
edges_gdf["geometry"] = edges_gdf["geometry"].apply(
lambda g: g.intersection(combined_bounds) if g is not None else g
)
edges_gdf = edges_gdf[~edges_gdf["geometry"].is_empty & edges_gdf["geometry"].is_valid]
# Keep only LineStrings / MultiLineStrings
edges_gdf = edges_gdf[edges_gdf.geometry.type.isin(["LineString", "MultiLineString"])]
print(f" Roads clipped to study area: {len(edges_gdf)} segments")
# Flatten MultiLineStrings to individual LineStrings
from shapely.geometry import LineString as ShapelyLineString
simple_edges = []
for idx, row in edges_gdf.iterrows():
geom = row.geometry
if geom.type == "MultiLineString":
for line in geom.geoms:
new_row = row.copy()
new_row["geometry"] = line
simple_edges.append(new_row)
else:
simple_edges.append(row)
roads_gdf = gpd.GeoDataFrame(simple_edges, crs=METRIC_CRS)
print(f" Simplified to {len(roads_gdf)} single-part line segments")
# ── Step 3: Find crossing points ───────────────────────────────────────────
print("\n[3/5] Detecting road / ward-boundary crossings…")
all_crossings = []
for idx, ward in wards_metric.iterrows():
ward_id = ward["ward_id"]
ward_name_en = ward["ward_name_en"]
ward_name = ward["ward_name"]
ward_boundary = ward.geometry.boundary
candidate_mask = roads_gdf.intersects(ward_boundary)
candidate_roads = roads_gdf[candidate_mask]
for ridx, road in candidate_roads.iterrows():
road_geom = road.geometry
inter = road_geom.intersection(ward_boundary)
if inter.is_empty:
continue
if inter.type == "Point":
all_crossings.append({
"geometry": inter,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
elif inter.type == "MultiPoint":
for pt in inter.geoms:
all_crossings.append({
"geometry": pt,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
elif inter.type in ("LineString", "MultiLineString"):
# Road runs along boundary — skip
continue
print(f" Found {len(all_crossings)} crossing points")
if len(all_crossings) > 0:
crossing_pts = gpd.GeoDataFrame(all_crossings, crs=METRIC_CRS)
else:
crossing_pts = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
crossing_pts.to_file(GPKG_PATH, layer="crossing_points", driver="GPKG")
print(f" ✓ Layer 'crossing_points' saved ({len(crossing_pts)} points)")
# ── Step 4: Create 50 m buffers ────────────────────────────────────────────
print("\n[4/5] Creating 50 m buffers and clipping to wards…")
if len(crossing_pts) > 0:
buffers = crossing_pts.copy()
buffers["geometry"] = buffers["geometry"].buffer(50.0)
buffers.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
print(f" ✓ Layer 'crossing_buffers_50m' saved ({len(buffers)} buffers)")
# Clip each buffer to the ward whose boundary it crossed
clipped_parts = []
for bidx, buf in buffers.iterrows():
ward_id = buf["ward_id"]
ward_geom = wards_metric.loc[wards_metric["ward_id"] == ward_id, "geometry"].iloc[0]
clipped = buf["geometry"].intersection(ward_geom)
if clipped.is_empty:
continue
if clipped.type == "MultiPolygon":
for poly in clipped.geoms:
clipped_parts.append({
"geometry": poly,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
else:
clipped_parts.append({
"geometry": clipped,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
if clipped_parts:
clipped_gdf = gpd.GeoDataFrame(clipped_parts, crs=METRIC_CRS)
clipped_gdf["geometry"] = clipped_gdf["geometry"].apply(
lambda g: make_valid(g) if not g.is_valid else g
)
clipped_gdf = clipped_gdf[~clipped_gdf["geometry"].is_empty]
clipped_gdf.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(f" ✓ Layer 'buffer_ward_intersection' saved ({len(clipped_gdf)} polygons)")
else:
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layer 'buffer_ward_intersection' saved (empty)")
else:
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layers 'crossing_buffers_50m' and 'buffer_ward_intersection' saved (empty)")
# ── Step 5: Ward crossing density table (WGS84) ────────────────────────────
print("\n[5/5] Computing crossing density per ward…")
if len(crossing_pts) > 0:
counts = crossing_pts.groupby("ward_id").size().reset_index(name="crossing_count")
else:
counts = pd.DataFrame({"ward_id": wards_gdf["ward_id"], "crossing_count": 0})
# Area from metric CRS
wards_metric["ward_area_km2"] = wards_metric.geometry.area / 1_000_000.0
density_gdf = wards_gdf[["ward_id", "ward_name_en", "ward_name", "geometry"]].copy()
density_gdf = density_gdf.merge(counts, on="ward_id", how="left")
density_gdf["crossing_count"] = density_gdf["crossing_count"].fillna(0).astype(int)
area_map = dict(zip(wards_metric["ward_id"], wards_metric["ward_area_km2"]))
density_gdf["ward_area_km2"] = density_gdf["ward_id"].map(area_map).round(4)
density_gdf["crossings_per_km2"] = np.where(
density_gdf["ward_area_km2"] > 0,
(density_gdf["crossing_count"] / density_gdf["ward_area_km2"]).round(4),
0.0
)
# Already in WGS84
density_gdf.to_file(GPKG_PATH, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f" ✓ Layer 'ward_crossing_density_wgs84' saved ({len(density_gdf)} wards)")
# ── Summary ─────────────────────────────────────────────────────────────────
print("\n" + "=" * 60)
print("SUMMARY")
print("=" * 60)
print(f" Wards: {len(wards_gdf)}")
print(f" Crossing points: {len(crossing_pts)}")
print(f" Buffers (50 m): {len(crossing_pts)}")
print(f" Density table: {len(density_gdf)} rows")
print(f"\nOutput: {GPKG_PATH}")
print("Done.")Wrote /work/solve.py (12603 chars).
9:59
9:59
Bash
cd /work && python solve.py 2>&1
============================================================
Tokyo 23 Special Wards — Road Crossing Analysis
============================================================
[1/5] Fetching ward boundaries from OSM…
Traceback (most recent call last):
File "/work/solve.py", line 70, in <module>
resp.raise_for_status()
~~~~~~~~~~~~~~~~~~~~~^^
File "/usr/local/lib/python3.14/site-packages/requests/models.py", line 1167, in raise_for_status
raise HTTPError(http_error_msg, response=self)
requests.exceptions.HTTPError: 406 Client Error: Not Acceptable for url: https://overpass-api.de/api/interpreter
10:12
The requests library is adding some default Accept header. Let me try with urllib instead, or use raw sockets, or set the Accept header to the right value.
Actually, looking at the successful test earlier, the successful request was:
```python
r = requests.post(url, data=payload, timeout=60,
headers={'User-Agent': 'OSM-GIS-Analysis/1.0'},
)
```
It worked without any Accept header. But the default requests library adds `Accept: */*`. Let me check if we can suppress that.
Actually, the issue might be that `requests` adds `Accept: */*` by default, which is being rejected. Let me try using `urllib` directly since that worked earlier.
Wait, the successful test was with `urllib.request` not `requests`. Let me check:
In test 5 (the successful one): `r = requests.post(url, data=payload, timeout=60, headers={'User-Agent': 'OSM-GIS-Analysis/1.0'})` returned 200.
But in test 4: `r = requests.post(url, data=payload, timeout=60, headers={'User-Agent': 'Mozilla/5.0', 'Accept': '*/*'})` returned 406.
So the issue is the Accept header. If we don't set it, requests adds `Accept: */*` by default. But when we explicitly set `{'User-Agent': 'OSM-GIS-Analysis/1.0'}`, it worked...
Wait, but the script also failed. Let me check if the requests library is somehow adding the Accept header. Let me just use urllib directly instead.
10:12
Bash
python3 -c "
python3 -c "
import requests
# Check what headers requests sends by default
url = 'https://overpass-api.de/api/interpreter'
payload = {'data': '[out:json];rel(id:1759466);out body geom;'}
r = requests.post(url, data=payload, timeout=60,
headers={'User-Agent': 'OSM-GIS-Analysis/1.0'})
print('Status:', r.status_code)
print(r.text[:200])
" 2>&1
Status: 200
{
"version": 0.6,
"generator": "Overpass API 0.7.62.11 87bfad18",
"osm3s": {
"timestamp_osm_base": "2026-06-16T22:06:00Z",
"copyright": "The data included in this document is from www.op
10:19
So the issue is that `requests` library in newer versions adds an `Accept` header by default. When we explicitly set `User-Agent` but don't set `Accept`, it might still send the default `Accept` header. Let me check what headers are being sent.
Actually, looking at my code: I removed the header entirely in the first test and it worked (with just User-Agent). But in the script I used `resp = requests.post(overpass_url, data={"data": query}, timeout=300)` without any headers. Let me check what happens.
Let me check what the requests library sends by default vs. when User-Agent is set.
10:19
Bash
python3 -c "
python3 -c " import requests # Check default headers url = 'https://httpbin.org/headers' r = requests.get(url) print(r.json()) " 2>&1
{'headers': {'Accept': '*/*', 'Accept-Encoding': 'gzip, deflate, zstd', 'Host': 'httpbin.org', 'User-Agent': 'python-requests/2.34.2', 'X-Amzn-Trace-Id': 'Root=1-6a31c93a-76e7d20421ad20ba4c58bdea'}}
10:25
Interesting - the requests library sends `Accept: */*` by default, and that works fine from the terminal. But in the script it fails. Let me try to figure out what's different. Maybe it's the User-Agent? Let me add a User-Agent header like the working test.
10:25
Edit
/work/solve.py
Edited /work/solve.py.
10:31
Wait, earlier the test with Mozilla User-Agent failed (406). But then without User-Agent it worked (200). Let me check... Actually the test that worked was:
```python
r = requests.post(url, data=payload, timeout=60,
headers={'User-Agent': 'OSM-GIS-Analysis/1.0'})
```
So the key is the User-Agent string. `python-requests/2.34.2` causes 406, but `OSM-GIS-Analysis/1.0` works. Let me use that.
10:31
Edit
/work/solve.py
Edited /work/solve.py.
10:37
Actually, the issue is that `python-requests/2.34.2` sends `Accept: */*` by default, which should be fine. But somehow it still gets 406. Let me check more carefully. The previous test that worked had `User-Agent: 'OSM-GIS-Analysis/1.0'` but no explicit Accept header. Let me test more. Actually the solution is simpler - just set a specific User-Agent that isn't "python-requests".
10:37
Bash
cd /work && python solve.py 2>&1
============================================================
Tokyo 23 Special Wards — Road Crossing Analysis
============================================================
[1/5] Fetching ward boundaries from OSM…
Retrieved 21 wards
✓ Layer 'wards_jgd' saved
[2/5] Fetching drivable road network from OSM…
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/urllib3/connection.py", line 204, in _new_conn
sock = connection.create_connection(
(self._dns_host, self.port),
...<2 lines>...
socket_options=self.socket_options,
)
File "/usr/local/lib/python3.14/site-packages/urllib3/util/connection.py", line 85, in create_connection
raise err
File "/usr/local/lib/python3.14/site-packages/urllib3/util/connection.py", line 73, in create_connection
sock.connect(sa)
~~~~~~~~~~~~^^^^
ConnectionRefusedError: [Errno 111] Connection refused
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 788, in urlopen
response = self._make_request(
conn,
...<10 lines>...
**response_kw,
)
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 488, in _make_request
raise new_e
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 464, in _make_request
self._validate_conn(conn)
~~~~~~~~~~~~~~~~~~~^^^^^^
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 1106, in _validate_conn
conn.connect()
~~~~~~~~~~~~^^
File "/usr/local/lib/python3.14/site-packages/urllib3/connection.py", line 759, in connect
self.sock = sock = self._new_conn()
~~~~~~~~~~~~~~^^
File "/usr/local/lib/python3.14/site-packages/urllib3/connection.py", line 219, in _new_conn
raise NewConnectionError(
self, f"Failed to establish a new connection: {e}"
) from e
urllib3.exceptions.NewConnectionError: HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/requests/adapters.py", line 696, in send
resp = conn.urlopen(
method=request.method,
...<9 lines>...
chunked=chunked,
)
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 842, in urlopen
retries = retries.increment(
method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
)
File "/usr/local/lib/python3.14/site-packages/urllib3/util/retry.py", line 543, in increment
raise MaxRetryError(_pool, url, reason) from reason # type: ignore[arg-type]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
urllib3.exceptions.MaxRetryError: HTTPSConnectionPool(host='overpass-api.de', port=443): Max retries exceeded with url: /api/interpreter (Caused by NewConnectionError("HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused"))
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/work/solve.py", line 144, in <module>
G = ox.graph_from_polygon(
wards_gdf.unary_union,
...<3 lines>...
truncate_by_edge=False
)
File "/usr/local/lib/python3.14/site-packages/osmnx/graph.py", line 496, in graph_from_polygon
G_buff = _create_graph(response_jsons, bidirectional)
File "/usr/local/lib/python3.14/site-packages/osmnx/graph.py", line 625, in _create_graph
for response_json in response_jsons:
^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/osmnx/_overpass.py", line 402, in _download_overpass_network
yield _overpass_request(OrderedDict(data=query_str))
~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/osmnx/_overpass.py", line 469, in _overpass_request
response = requests.post(
url,
...<3 lines>...
**settings.requests_kwargs,
)
File "/usr/local/lib/python3.14/site-packages/requests/api.py", line 134, in post
return request("post", url, data=data, json=json, **kwargs)
File "/usr/local/lib/python3.14/site-packages/requests/api.py", line 71, in request
return session.request(method=method, url=url, **kwargs)
~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/requests/sessions.py", line 651, in request
resp = self.send(prep, **send_kwargs)
File "/usr/local/lib/python3.14/site-packages/requests/sessions.py", line 784, in send
r = adapter.send(request, **kwargs)
File "/usr/local/lib/python3.14/site-packages/requests/adapters.py", line 729, in send
raise ConnectionError(e, request=request)
requests.exceptions.ConnectionError: HTTPSConnectionPool(host='overpass-api.de', port=443): Max retries exceeded with url: /api/interpreter (Caused by NewConnectionError("HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused"))
12:06
The OSMnx library is also hitting the overpass-api.de endpoint with its own requests. The overpass-api.de seems to be rate-limiting or blocking connections after our first successful request. Let me try a different approach: instead of using OSMnx for the road network, I'll make direct Overpass API queries to get the roads, but I'll handle the query and parsing myself. I'll also avoid overpass-api.de for the second query by using a different approach. Actually, the issue is that overpass-api.de is rate-limiting us. The first query for 23 wards works fine (21 retrieved, 2 failed - maybe Setagaya and another ward have complex boundaries). But the second query for roads gets connection refused. Let me use a different Overpass endpoint or use OSMnx differently. Actually, looking at this more carefully, the overpass-api.de seems to have intermittent connectivity issues. Let me try using overpass.kumi.systems instead, or build the queries directly. Wait, let me re-think. The issue might be that OSMnx internally connects to overpass-api.de (the default) and doesn't use the custom User-Agent. Let me check if there's a way to configure OSMnx to use the working endpoint with proper headers. Actually, looking at the osmnx source code, it uses `settings.requests_kwargs` for extra kwargs. Let me try setting the User-Agent there. But first, let me also understand why 2 of the 23 wards were not retrieved. Let me check which ones are missing.
12:06
Bash
python3 -c "
python3 -c "
import requests
# Check OSMnx default settings
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# First query: all 23 wards
rel_ids = [1759466, 1759470, 1759475, 1759473, 1759474, 1759476, 1759477, 1759478, 1759479, 1759480, 1759481, 1759482, 1759483, 1759484, 1759485, 1759486, 1759487, 1759488, 1759489, 1759490, 1759491, 1759492, 1759493]
rel_id_str = ','.join(str(i) for i in rel_ids)
query = f'[out:json][timeout:300];rel(id:{rel_id_str});out body geom;'
resp = requests.post(url, data={'data': query}, timeout=300, headers=headers)
print(f'Status: {resp.status_code}')
data = resp.json()
print(f'Elements: {len(data.get(\"elements\", []))}')
for el in data.get('elements', []):
if el['type'] == 'relation':
print(f' Relation {el[\"id\"]} - members: {len(el.get(\"members\", []))}')
" 2>&1
Status: 200 Elements: 21 Relation 1759466 - members: 4 Relation 1759470 - members: 3 Relation 1759473 - members: 5 Relation 1759474 - members: 192 Relation 1759475 - members: 18 Relation 1759477 - members: 112 Relation 1759479 - members: 19 Relation 1759480 - members: 69 Relation 1759481 - members: 2 Relation 1759482 - members: 9 Relation 1759483 - members: 5 Relation 1759484 - members: 2 Relation 1759485 - members: 2 Relation 1759486 - members: 2 Relation 1759487 - members: 2 Relation 1759488 - members: 2 Relation 1759489 - members: 2 Relation 1759490 - members: 2 Relation 1759491 - members: 2 Relation 1759492 - members: 2 Relation 1759493 - members: 2
12:16
We got 21 out of 23 wards. Let me check which are missing: 1759476 (Taito), 1759478 (Koto). Let me query them individually.
12:16
Bash
python3 -c "
python3 -c "
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Check missing wards
for rid in [1759476, 1759478]:
query = f'[out:json][timeout:30];rel(id:{rid});out body geom;'
resp = requests.post(url, data={'data': query}, timeout=60, headers=headers)
print(f'Relation {rid}: status={resp.status_code}')
data = resp.json()
els = data.get('elements', [])
print(f' Elements: {len(els)}')
for el in els:
print(f' Type: {el[\"type\"]}, tags: {el.get(\"tags\", {})}')
if el['type'] == 'relation':
print(f' Members: {len(el.get(\"members\", []))}')
" 2>&1
Relation 1759476: status=200 Elements: 0 Relation 1759478: status=200 Elements: 0
12:25
Two wards are missing - Taito (1759476) and Koto (1759478). These might have different relation IDs or might need a different approach. Let me check.
12:25
(no output)
Bash
python3 -c "
python3 -c "
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Search for Taito and Koto wards
query = '''
[out:json][timeout:30];
area[id=\"3600175946\"]; // Tokyo prefecture
rel(area)[boundary=administrative][admin_level=9];
out tags;
'''
resp = requests.post(url, data={'data': query}, timeout=60, headers=headers)
data = resp.json()
for el in data.get('elements', []):
name = el.get('tags', {}).get('name', '')
name_en = el.get('tags', {}).get('name:en', '')
rid = el['id']
print(f' {rid}: {name} / {name_en}')
" 2>&1
12:38
The area query didn't work. Let me try a different approach to find the missing ward relations.
12:38
Bash
python3 -c "
python3 -c "
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Check the two missing relations with 'out body' (no geom, just to check if they exist)
query = '[out:json];rel(id:1759476,1759478);out body;'
resp = requests.post(url, data={'data': query}, timeout=60, headers=headers)
data = resp.json()
for el in data.get('elements', []):
print(f'id={el[\"id\"]}, type={el[\"type\"]}, tags={el.get(\"tags\", {})}')
print(f'Total elements: {len(data.get(\"elements\", []))}')
" 2>&1
Total elements: 0
12:49
The Taito (1759476) and Koto (1759478) relations are returning no results. These might be old/deleted relations. Let me search for them by name.
12:49
Bash
python3 -c "
python3 -c "
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Search by name
query = '''
[out:json][timeout:30];
(
rel[\"name:en\"~\"Taito|Koto|台東|江東\"][boundary=administrative]({{bbox}});
);
out body;
'''
# Use a bbox that covers roughly Tokyo
bbox = '35.5,139.5,35.8,139.95'
query2 = f'''
[out:json][timeout:30];
(
rel[\"name:en\"~\"Taito|Koto\"][boundary=administrative]({bbox});
rel[\"name\"~\"台東区|江東区\"][boundary=administrative]({bbox});
);
out body;
'''
resp = requests.post(url, data={'data': query2}, timeout=60, headers=headers)
data = resp.json()
for el in data.get('elements', []):
tags = el.get('tags', {})
print(f'id={el[\"id\"]}, name={tags.get(\"name\",\"\")}, name:en={tags.get(\"name:en\",\"\")}, admin_level={tags.get(\"admin_level\",\"\")}')
" 2>&1
id=1758888, name=台東区, name:en=Taito, admin_level=7 id=3554015, name=江東区, name:en=Koto, admin_level=7 id=15920350, name=台東四丁目, name:en=Taito 4-chome, admin_level=10 id=15920351, name=台東三丁目, name:en=Taito 3-chome, admin_level=10 id=15920352, name=台東二丁目, name:en=Taito 2-chome, admin_level=10 id=16929798, name=江東橋二丁目, name:en=Kotobashi 2, admin_level=10 id=16929801, name=江東橋, name:en=Kotobashi, admin_level=9 id=18158501, name=寿四丁目, name:en=Kotobuki 4-chome, admin_level=10 id=18158502, name=寿三丁目, name:en=Kotobuki 3-chome, admin_level=10 id=18158503, name=寿二丁目, name:en=Kotobuki 2-chome, admin_level=10 id=18158504, name=寿一丁目, name:en=Kotobuki 1-chome, admin_level=10 id=18158505, name=寿, name:en=Kotobuki, admin_level=9 id=18158950, name=台東一丁目, name:en=Taito 1-chome, admin_level=10 id=18158951, name=台東, name:en=Taito, admin_level=9
12:58
So the relation IDs I used for Taito (1759476) and Koto (1759478) seem to be wrong. The correct IDs are: - Taito: 1758888 (admin_level=7) - Koto: 3554015 (admin_level=7) Note that admin_level=7 for special wards in Japan, not 9. Let me check all the others.
12:58
Bash
python3 -c "
python3 -c "
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Find all 23 special wards of Tokyo
bbox = '35.5,139.5,35.8,139.95'
query = f'''
[out:json][timeout:60];
(
rel[\"admin_level\"=\"7\"][boundary=administrative]({bbox});
rel[\"admin_level\"=\"9\"][boundary=administrative]({bbox});
);
out tags;
'''
resp = requests.post(url, data={'data': query}, timeout=120, headers=headers)
data = resp.json()
wards = {}
for el in data.get('elements', []):
tags = el.get('tags', {})
name = tags.get('name', '')
name_en = tags.get('name:en', '')
admin = tags.get('admin_level', '')
if name_en and 'ward' in name_en.lower():
wards[el['id']] = (name, name_en, admin)
elif name and '区' in name and admin in ('7', '9'):
wards[el['id']] = (name, name_en, admin)
for rid, (name, name_en, admin) in sorted(wards.items()):
print(f'{rid}: {name} ({name_en}) admin_level={admin}')
" 2>&1
1543055: 杉並区 (Suginami) admin_level=7 1543056: 中野区 (Nakano) admin_level=7 1758858: 新宿区 (Shinjuku) admin_level=7 1758878: 文京区 (Bunkyō) admin_level=7 1758888: 台東区 (Taito) admin_level=7 1758891: 墨田区 (Sumida) admin_level=7 1758897: 中央区 (Chuo) admin_level=7 1758936: 目黒区 (Meguro) admin_level=7 1758947: 大田区 (Ota) admin_level=7 1759474: 世田谷区 (Setagaya) admin_level=7 1759477: 渋谷区 (Shibuya) admin_level=7 1759506: 豊島区 (Toshima) admin_level=7 1760038: 北区 (Kita) admin_level=7 1760040: 荒川区 (Arakawa) admin_level=7 1760078: 板橋区 (Itabashi) admin_level=7 1760119: 練馬区 (Nerima) admin_level=7 1760124: 足立区 (Adachi) admin_level=7 1761717: 港区 (Minato) admin_level=7 1761718: 葛飾区 (Katsushika) admin_level=7 1761742: 千代田区 (Chiyoda) admin_level=7 1761743: 江戸川区 (Edogawa) admin_level=7 3554015: 江東区 (Koto) admin_level=7 3554304: 品川区 (Shinagawa) admin_level=7
13:14
So the Tokyo special wards have admin_level=7 in OSM (not 9). The original relation IDs I had were a mix. Let me get the correct relation IDs for all 23 wards and their geometries. Let me also get the actual geometries. Key relations: - 1758888 - Taito (台東区) - 3554015 - Koto (江東区) - 1761742 - Chiyoda (千代田区) - 1758897 - Chuo (中央区) - 1761717 - Minato (港区) - 1758858 - Shinjuku (新宿区) - 1758878 - Bunkyo (文京区) - 1758891 - Sumida (墨田区) - 1759477 - Shibuya (渋谷区) - 1543056 - Nakano (中野区) - 1543055 - Suginami (杉並区) - 1759506 - Toshima (豊島区) - 1760038 - Kita (北区) - 1760040 - Arakawa (荒川区) - 1760078 - Itabashi (板橋区) - 1760119 - Nerima (練馬区) - 1760124 - Adachi (足立区) - 1761718 - Katsushika (葛飾区) - 1761743 - Edogawa (江戸川区) - 1758936 - Meguro (目黒区) - 1759474 - Setagaya (世田谷区) <-- wait, 1759474 is Setagaya but before I had it as Bunkyo - 1758947 - Ota (大田区) - 3554304 - Shinagawa (品川区) Let me compile the correct mapping and update the script.
13:14
Bash
python3 -c "
python3 -c "
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Get all 23 wards with their geometries
rel_ids = [1761742, 1758897, 1761717, 1758858, 1758878, 1758888, 1758891, 3554015, 3554304, 1758936, 1758947, 1759474, 1759477, 1543056, 1543055, 1759506, 1760038, 1760040, 1760078, 1760119, 1760124, 1761718, 1761743]
rel_id_str = ','.join(str(i) for i in rel_ids)
# First just get tags
query = f'[out:json][timeout:120];rel(id:{rel_id_str});out tags;'
resp = requests.post(url, data={'data': query}, timeout=120, headers=headers)
data = resp.json()
for el in data.get('elements', []):
tags = el.get('tags', {})
print(f'{el[\"id\"]}: {tags.get(\"name\",\"\")} / {tags.get(\"name:en\",\"\")} / {tags.get(\"name:ja\",\"\")}')
" 2>&1
1543055: 杉並区 / Suginami / 杉並区 1543056: 中野区 / Nakano / 中野区 1758858: 新宿区 / Shinjuku / 新宿区 1758878: 文京区 / Bunkyō / 文京区 1758888: 台東区 / Taito / 台東区 1758891: 墨田区 / Sumida / 墨田区 1758897: 中央区 / Chuo / 中央区 1758936: 目黒区 / Meguro / 目黒区 1758947: 大田区 / Ota / 大田区 1759474: 世田谷区 / Setagaya / 世田谷区 1759477: 渋谷区 / Shibuya / 渋谷区 1759506: 豊島区 / Toshima / 豊島区 1760038: 北区 / Kita / 北区 1760040: 荒川区 / Arakawa / 荒川区 1760078: 板橋区 / Itabashi / 板橋区 1760119: 練馬区 / Nerima / 練馬区 1760124: 足立区 / Adachi / 足立区 1761717: 港区 / Minato / 港区 1761718: 葛飾区 / Katsushika / 葛飾区 1761742: 千代田区 / Chiyoda / 千代田区 1761743: 江戸川区 / Edogawa / 江戸川区 3554015: 江東区 / Koto / 江東区 3554304: 品川区 / Shinagawa / 品川区
13:27
So the actual OSM relation IDs for Tokyo's 23 special wards (admin_level=7) are different from what I initially used. Let me also check relation 1759476 (Taito) and 1759478 (Koto) - these were wrong. The correct IDs are: From the search: - 1761742: Chiyoda (千代田区) - 1758897: Chuo (中央区) - 1761717: Minato (港区) - 1758858: Shinjuku (新宿区) - 1758878: Bunkyo (文京区) - 1758888: Taito (台東区) - 1758891: Sumida (墨田区) - 3554015: Koto (江東区) - 3554304: Shinagawa (品川区) - 1758936: Meguro (目黒区) - 1758947: Ota (大田区) - 1759474: Setagaya (世田谷区) - 1759477: Shibuya (渋谷区) - 1543056: Nakano (中野区) - 1543055: Suginami (杉並区) - 1759506: Toshima (豊島区) - 1760038: Kita (北区) - 1760040: Arakawa (荒川区) - 1760078: Itabashi (板橋区) - 1760119: Nerima (練馬区) - 1760124: Adachi (足立区) - 1761718: Katsushika (葛飾区) - 1761743: Edogawa (江戸川区) My original IDs like 1759466 (Chiyoda), 1759470 (Chuo), etc. were wrong. Let me update the script with the correct IDs. Now for the OSMnx step: the problem is that OSMnx uses its own requests settings which might trigger the same 406/connection issues. Let me check if I can configure OSMnx to use a different User-Agent or override the overpass endpoint.
13:27
13:27
Bash
python3 -c "
python3 -c " import osmnx as ox print(ox.settings.requests_kwargs) print(ox.settings.overpass_endpoint) print(ox.settings.default_user_agent) " 2>&1
{}
Traceback (most recent call last):
File "<string>", line 4, in <module>
print(ox.settings.overpass_endpoint)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AttributeError: module 'osmnx.settings' has no attribute 'overpass_endpoint'
13:37
Bash
python3 -c "
python3 -c " import osmnx as ox print(dir(ox.settings)) " 2>&1 | tr ',' '\n'
['Any' 'TYPE_CHECKING' '__annotations__' '__builtins__' '__cached__' '__conditional_annotations__' '__doc__' '__file__' '__loader__' '__name__' '__package__' '__spec__' 'all_oneway' 'annotations' 'bidirectional_network_types' 'cache_folder' 'cache_only_mode' 'data_folder' 'default_access' 'default_crs' 'doh_url_template' 'elevation_url_template' 'http_accept_language' 'http_referer' 'http_user_agent' 'imgs_folder' 'lg' 'log_console' 'log_file' 'log_filename' 'log_level' 'log_name' 'logs_folder' 'max_query_area_size' 'nominatim_key' 'nominatim_url' 'overpass_memory' 'overpass_rate_limit' 'overpass_settings' 'overpass_url' 'requests_kwargs' 'requests_timeout' 'use_cache' 'useful_tags_node' 'useful_tags_way']
13:46
So the settings are `overpass_url` and `requests_kwargs`. Let me set the requests_kwargs to include a proper User-Agent, or set the http_user_agent.
13:46
Bash
python3 -c "
python3 -c "
import osmnx as ox
print('default_user_agent:', ox.settings.http_user_agent)
print('overpass_url:', ox.settings.overpass_url)
print('requests_kwargs:', ox.settings.requests_kwargs)
" 2>&1
default_user_agent: OSMnx Python package (https://github.com/gboeing/osmnx)
overpass_url: https://overpass-api.de/api
requests_kwargs: {}
14:24
The issue is that when osmnx uses requests, it doesn't set a custom User-Agent, so the Overpass API returns 406. Let me override the user-agent setting in osmnx.
14:24
Bash
python3 -c "
python3 -c "
import osmnx as ox
ox.settings.http_user_agent = 'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36'
# Test with small area
try:
G = ox.graph_from_bbox(35.71, 35.69, 139.78, 139.76, network_type='drive')
print(f'Nodes: {len(G.nodes)}, Edges: {len(G.edges)}')
except Exception as e:
print(f'Error: {e}')
" 2>&1
Error: graph_from_bbox() takes 1 positional argument but 4 positional arguments (and 1 keyword-only argument) were given
14:42
Bash
python3 -c "
python3 -c "
import osmnx as ox
ox.settings.http_user_agent = 'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36'
# Test with small area - correct API
bbox = (139.76, 35.69, 139.78, 35.71) # left, bottom, right, top
try:
G = ox.graph_from_bbox(bbox, network_type='drive')
print(f'Nodes: {len(G.nodes)}, Edges: {len(G.edges)}')
except Exception as e:
print(f'Error: {type(e).__name__}: {e}')
" 2>&1
Error: ConnectionError: HTTPSConnectionPool(host='overpass-api.de', port=443): Max retries exceeded with url: /api/interpreter (Caused by NewConnectionError("HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused"))
15:53
The overpass-api.de seems to be refusing connections intermittently. However, I was able to connect directly via Python requests with a custom User-Agent. The issue with OSMnx is that it's using the python-requests default User-Agent which might be blocked. Let me try to set the requests_kwargs in OSMnx settings to use a different User-Agent.
15:53
Bash
python3 -c "
python3 -c "
import osmnx as ox
ox.settings.http_user_agent = 'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36'
ox.settings.requests_kwargs = {'headers': {'User-Agent': 'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36'}}
bbox = (139.76, 35.69, 139.78, 35.71)
try:
G = ox.graph_from_bbox(bbox, network_type='drive')
print(f'Nodes: {len(G.nodes)}, Edges: {len(G.edges)}')
except Exception as e:
print(f'Error: {type(e).__name__}: {e}')
" 2>&1
Error: TypeError: requests.api.get() got multiple values for keyword argument 'headers'
16:04
response
16:04
Bash
python3 -c "
python3 -c "
import osmnx as ox
ox.settings.http_user_agent = 'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36'
# Check what osmnx passes to requests
print('user_agent:', repr(ox.settings.http_user_agent))
# Try with a simple request first to see if network works
import requests
r = requests.post('https://overpass-api.de/api/interpreter',
data={'data': '[out:json];node(35.69,139.76,35.71,139.78);out;'},
timeout=30,
headers={'User-Agent': 'Mozilla/5.0'})
print(f'Status: {r.status_code}')
print(r.text[:200])
" 2>&1
user_agent: 'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36' Status: 406 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> <html><head> <title>406 Not Acceptable</title> </head><body> <h1>Not Acceptable</h1> <p>An appropriate repres
16:16
So the 406 error is back! It seems like when we set a User-Agent in the requests headers, the overpass API sends a 406. But earlier, when I used `{'User-Agent': 'OSM-GIS-Analysis/1.0'}` it worked. Let me check more carefully...
Actually wait, the 406 error was happening specifically with the requests library's default headers. Let me re-check: when I tried with `{'User-Agent': 'OSM-GIS-Analysis/1.0'}` for the single relation, it returned 200. But now when I use `{'User-Agent': 'Mozilla/5.0'}` it returns 406? That's inconsistent...
Let me re-test more carefully.
16:16
Bash
python3 -c "
python3 -c "
import requests
# Test exact same query with different user agents
url = 'https://overpass-api.de/api/interpreter'
query = '[out:json];node(35.69,139.76,35.71,139.78);out;'
for ua in ['Mozilla/5.0', 'OSM-GIS-Analysis/1.0', 'python-requests/2.34.2', 'wget']:
try:
r = requests.post(url, data={'data': query}, timeout=30, headers={'User-Agent': ua})
print(f'UA=\"{ua}\": Status={r.status_code}, len={len(r.text)}')
except Exception as e:
print(f'UA=\"{ua}\": Error={e}')
" 2>&1
UA="Mozilla/5.0": Status=406, len=371 UA="OSM-GIS-Analysis/1.0": Status=200, len=10409187 UA="python-requests/2.34.2": Status=406, len=371 UA="wget": Status=200, len=10409187
16:33
Interesting! So the User-Agent "OSM-GIS-Analysis/1.0" and "wget" work, but "Mozilla/5.0" and "python-requests/2.34.2" cause 406 errors. The issue with using OSMnx is that it uses its own User-Agent ("OSMnx Python package...") which I need to check what status that gets.
Let me check:
16:33
Bash
python3 -c "
python3 -c "
import requests
url = 'https://overpass-api.de/api/interpreter'
query = '[out:json];node(35.69,139.76,35.71,139.78);out;'
for ua in ['OSMnx Python package (https://github.com/gboeing/osmnx)']:
r = requests.post(url, data={'data': query}, timeout=30, headers={'User-Agent': ua})
print(f'UA=\"{ua}\": Status={r.status_code}')
" 2>&1
UA="OSMnx Python package (https://github.com/gboeing/osmnx)": Status=429
16:49
429 means too many requests. Let me set the user agent for osmnx to something that works. The issue is that the environment is rate-limited by overpass-api.de when using certain UA strings. Let me set the user agent to something that works. Actually, let me check the osmnx source to see exactly how it sends requests.
16:49
Bash
python3 -c "
python3 -c " # Let's look at how osmnx sends requests import inspect import osmnx._overpass as op src = inspect.getsource(op._overpass_request) print(src[:1500]) " 2>&1
def _overpass_request(data: OrderedDict[str, Any]) -> dict[str, Any]:
"""
Send a HTTP POST request to the Overpass API and return response.
Parameters
----------
data
Key-value pairs of parameters.
Returns
-------
response_json
The Overpass API's response.
"""
# resolve url to same IP even if there is server round-robin redirecting
_http._config_dns(settings.overpass_url)
# prepare the Overpass API URL and see if request already exists in cache
url = settings.overpass_url.rstrip("/") + "/interpreter"
prepared_url = str(requests.Request("GET", url, params=data).prepare().url)
cached_response_json = _http._retrieve_from_cache(prepared_url)
if isinstance(cached_response_json, dict):
return cached_response_json
# pause then request this URL
pause = _get_overpass_pause(settings.overpass_url)
hostname = _http._hostname_from_url(url)
msg = f"Pausing {pause} second(s) before making HTTP POST request to {hostname!r}"
utils.log(msg, level=lg.INFO)
time.sleep(pause)
# transmit the HTTP POST request
msg = f"Post {prepared_url} with timeout={settings.requests_timeout}"
utils.log(msg, level=lg.INFO)
response = requests.post(
url,
data=data,
timeout=settings.requests_timeout,
headers=_http._get_http_headers(),
**settings.requests_kwargs,
)
# handle 429 and 504 errors by pausing then recursively re-trying request
16:57
Bash
python3 -c "
python3 -c " import osmnx._http as _http src = inspect.getsource(_http._get_http_headers) print(src) " 2>&1
Traceback (most recent call last):
File "<string>", line 3, in <module>
src = inspect.getsource(_http._get_http_headers)
^^^^^^^
NameError: name 'inspect' is not defined. Did you forget to import 'inspect'?
17:04
Bash
python3 -c "
python3 -c " import inspect, osmnx._http as _http src = inspect.getsource(_http._get_http_headers) print(src) " 2>&1
def _get_http_headers(
*,
user_agent: str | None = None,
referer: str | None = None,
accept_language: str | None = None,
) -> dict[str, str]:
"""
Update the default requests HTTP headers with OSMnx information.
Parameters
----------
user_agent
The user agent. If None, use `settings.http_user_agent` value.
referer
The referer. If None, use `settings.http_referer` value.
accept_language
The accept language. If None, use `settings.http_accept_language`
value.
Returns
-------
headers
The updated HTTP headers.
"""
if user_agent is None:
user_agent = settings.http_user_agent
if referer is None:
referer = settings.http_referer
if accept_language is None:
accept_language = settings.http_accept_language
info = {"User-Agent": user_agent, "referer": referer, "Accept-Language": accept_language}
headers = dict(requests.utils.default_headers())
headers.update(info)
return headers
17:14
The issue is that `requests.utils.default_headers()` includes `Accept: */*` which the overpass-api.de server apparently doesn't like (406 error). When we use `User-Agent: 'OSM-GIS-Analysis/1.0'` or `wget`, it works fine.
Looking at the code, `_get_http_headers()` calls `requests.utils.default_headers()` which returns `{'User-Agent': 'python-requests/2.34.2', 'Accept-Encoding': 'gzip, deflate, zstd', 'Accept': '*/*', 'Connection': 'keep-alive'}` and then updates with OSMnx's user agent.
So the Accept: */* is causing the 406. We need to override the headers to not include Accept, or set it to something the server accepts.
We can use `ox.settings.requests_kwargs` to override headers... but earlier we saw it caused a conflict. Let me look at how to fix this properly.
Actually, looking at the code:
```python
response = requests.post(
url,
data=data,
timeout=settings.requests_timeout,
headers=_http._get_http_headers(),
**settings.requests_kwargs,
)
```
The headers dict from `_get_http_headers()` is passed as `headers=`, then `**settings.requests_kwargs` is unpacked. If `requests_kwargs` contains `headers`, there would be a conflict.
But we could monkey-patch `_get_http_headers` or set `http_user_agent` to something that works and also modify the default headers behavior.
Actually, looking more carefully: when User-Agent is 'OSM-GIS-Analysis/1.0', the requests.default_headers() still includes Accept: */* but it works fine. So the issue is specifically with the default User-Agent from osmnx (which is 'OSMnx Python package...') that's rate-limited or blocked.
Wait, earlier we tested:
- UA="OSM-GIS-Analysis/1.0": Status=200
- UA="OSMnx Python package (https://github.com/gboeing/osmnx)": Status=429
So the OSMnx default UA is getting rate-limited! Let's just change the user agent.
But wait, there's also the `Accept: */*` issue - when we tested with `Mozilla/5.0` it gave 406. Let me re-test with just changing the user agent in osmnx settings.
17:14
Bash
python3 -c "
python3 -c "
import osmnx as ox
ox.settings.http_user_agent = 'OSM-GIS-Analysis/1.0'
bbox = (139.76, 35.69, 139.78, 35.71)
try:
G = ox.graph_from_bbox(bbox, network_type='drive')
print(f'Nodes: {len(G.nodes)}, Edges: {len(G.edges)}')
except Exception as e:
print(f'Error: {type(e).__name__}: {e}')
" 2>&1
Nodes: 1296, Edges: 2660
19:02
The key insight: OSMnx was using the default User-Agent "OSMnx Python package..." which got rate-limited (429). We need to set a custom User-Agent. Let me rewrite the solution with: 1. The correct 23 ward relation IDs (I found them from the Overpass query) 2. Setting the OSMnx user agent to something that works 3. Using OSMnx to download road network (which handles the overpass queries internally)
19:02
19:02
Write
/work/solve.py
content · 323 lines · py
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
import requests
import warnings
warnings.filterwarnings('ignore')
ox.settings.use_cache = True
ox.settings.log_console = False
ox.settings.timeout = 300
ox.settings.http_user_agent = "OSM-GIS-Analysis/1.0"
METRIC_CRS = "EPSG:6679" # JGD2011 / Japan Plane Rectangular CS IX (Tokyo)
WGS84_CRS = "EPSG:4326"
GPKG_PATH = "/work/tokyo_crossings.gpkg"
print("=" * 60)
print("Tokyo 23 Special Wards — Road Crossing Analysis")
print("=" * 60)
# ── Step 1: Fetch the 23 special-ward boundaries from Overpass ─────────────
print("\n[1/5] Fetching ward boundaries from OSM…")
# Correct OSM relation IDs for Tokyo's 23 special wards (admin_level=7)
WARD_INFO = {
1761742: ("Chiyoda Ward", "千代田区"),
1758897: ("Chuo Ward", "中央区"),
1761717: ("Minato Ward", "港区"),
1758858: ("Shinjuku Ward", "新宿区"),
1758878: ("Bunkyo Ward", "文京区"),
1758888: ("Taito Ward", "台東区"),
1758891: ("Sumida Ward", "墨田区"),
3554015: ("Koto Ward", "江東区"),
3554304: ("Shinagawa Ward", "品川区"),
1758936: ("Meguro Ward", "目黒区"),
1758947: ("Ota Ward", "大田区"),
1759474: ("Setagaya Ward", "世田谷区"),
1759477: ("Shibuya Ward", "渋谷区"),
1543056: ("Nakano Ward", "中野区"),
1543055: ("Suginami Ward", "杉並区"),
1759506: ("Toshima Ward", "豊島区"),
1760038: ("Kita Ward", "北区"),
1760040: ("Arakawa Ward", "荒川区"),
1760078: ("Itabashi Ward", "板橋区"),
1760119: ("Nerima Ward", "練馬区"),
1760124: ("Adachi Ward", "足立区"),
1761718: ("Katsushika Ward","葛飾区"),
1761743: ("Edogawa Ward", "江戸川区"),
}
rel_ids = list(WARD_INFO.keys())
overpass_url = "https://overpass-api.de/api/interpreter"
rel_id_str = ",".join(str(i) for i in rel_ids)
query = f"""
[out:json][timeout:300];
rel(id:{rel_id_str});
out body geom;
"""
resp = requests.post(overpass_url, data={"data": query}, timeout=300,
headers={"User-Agent": "OSM-GIS-Analysis/1.0"})
resp.raise_for_status()
data = resp.json()
from shapely.geometry import Polygon as ShapelyPolygon, MultiPolygon
from shapely.validation import make_valid
ward_records = []
for el in data.get("elements", []):
if el["type"] != "relation":
continue
rid = el["id"]
name_en, name_ja = WARD_INFO[rid]
members = el.get("members", [])
outer_rings = []
inner_rings = []
for m in members:
if m.get("role") == "outer" and "geometry" in m:
coords = [(p["lon"], p["lat"]) for p in m["geometry"]]
if len(coords) >= 3:
outer_rings.append(coords)
elif m.get("role") == "inner" and "geometry" in m:
coords = [(p["lon"], p["lat"]) for p in m["geometry"]]
if len(coords) >= 3:
inner_rings.append(coords)
if not outer_rings:
print(f" ⚠ No outer rings for relation {rid}, skipping")
continue
if len(outer_rings) == 1:
poly = ShapelyPolygon(outer_rings[0], inner_rings if inner_rings else None)
else:
polys = []
for o_ring in outer_rings:
o_poly = ShapelyPolygon(o_ring)
holes_for_this = []
keep = []
for i_ring in inner_rings:
i_poly = ShapelyPolygon(i_ring)
if o_poly.contains(i_poly.representative_point()):
holes_for_this.append(i_ring)
else:
keep.append(i_ring)
inner_rings = keep
polys.append(ShapelyPolygon(o_ring, holes_for_this if holes_for_this else None))
poly = MultiPolygon(polys) if len(polys) > 1 else polys[0]
if not poly.is_valid:
poly = make_valid(poly)
ward_records.append({
"ward_id": rid,
"ward_name_en": name_en,
"ward_name": name_ja,
"geometry": poly
})
wards_gdf = gpd.GeoDataFrame(ward_records, crs=WGS84_CRS)
print(f" Retrieved {len(wards_gdf)} wards")
# Project to metric CRS and save layer 1
wards_metric = wards_gdf.to_crs(METRIC_CRS)
wards_metric.to_file(GPKG_PATH, layer="wards_jgd", driver="GPKG")
print(" ✓ Layer 'wards_jgd' saved")
# ── Step 2: Get drivable road network (excluding service roads) ────────────
print("\n[2/5] Fetching drivable road network from OSM…")
# Use OSMnx with the combined ward polygon
combined_bounds_wgs84 = wards_gdf.unary_union
G = ox.graph_from_polygon(
combined_bounds_wgs84,
network_type="drive",
simplify=False,
retain_all=True,
truncate_by_edge=False
)
print(f" Graph loaded: {len(G.nodes)} nodes, {len(G.edges)} edges")
# Convert edges to GeoDataFrame
edges_gdf = ox.graph_to_gdfs(G, nodes=False, edges=True)
edges_gdf = edges_gdf.to_crs(METRIC_CRS)
print(f" Raw edges in metric CRS: {len(edges_gdf)}")
# Clip to combined ward boundary
combined_bounds = wards_metric.unary_union
edges_gdf = edges_gdf[edges_gdf.intersects(combined_bounds)].copy()
edges_gdf["geometry"] = edges_gdf["geometry"].apply(
lambda g: g.intersection(combined_bounds) if g is not None else g
)
edges_gdf = edges_gdf[~edges_gdf["geometry"].is_empty & edges_gdf["geometry"].is_valid]
edges_gdf = edges_gdf[edges_gdf.geometry.type.isin(["LineString", "MultiLineString"])]
print(f" Roads clipped to study area: {len(edges_gdf)} segments")
# Flatten MultiLineStrings
from shapely.geometry import LineString as ShapelyLineString
simple_edges = []
for idx, row in edges_gdf.iterrows():
geom = row.geometry
if geom.type == "MultiLineString":
for line in geom.geoms:
new_row = row.copy()
new_row["geometry"] = line
simple_edges.append(new_row)
else:
simple_edges.append(row)
roads_gdf = gpd.GeoDataFrame(simple_edges, crs=METRIC_CRS)
print(f" Simplified to {len(roads_gdf)} single-part line segments")
# ── Step 3: Find crossing points ───────────────────────────────────────────
print("\n[3/5] Detecting road / ward-boundary crossings…")
all_crossings = []
for idx, ward in wards_metric.iterrows():
ward_id = ward["ward_id"]
ward_name_en = ward["ward_name_en"]
ward_name = ward["ward_name"]
ward_boundary = ward.geometry.boundary
candidate_mask = roads_gdf.intersects(ward_boundary)
candidate_roads = roads_gdf[candidate_mask]
for ridx, road in candidate_roads.iterrows():
road_geom = road.geometry
inter = road_geom.intersection(ward_boundary)
if inter.is_empty:
continue
if inter.type == "Point":
all_crossings.append({
"geometry": inter,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
elif inter.type == "MultiPoint":
for pt in inter.geoms:
all_crossings.append({
"geometry": pt,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
# Skip LineString/MultiLineString (road runs along boundary)
print(f" Found {len(all_crossings)} crossing points")
if len(all_crossings) > 0:
crossing_pts = gpd.GeoDataFrame(all_crossings, crs=METRIC_CRS)
else:
crossing_pts = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
crossing_pts.to_file(GPKG_PATH, layer="crossing_points", driver="GPKG")
print(f" ✓ Layer 'crossing_points' saved ({len(crossing_pts)} points)")
# ── Step 4: Create 50 m buffers ────────────────────────────────────────────
print("\n[4/5] Creating 50 m buffers and clipping to wards…")
if len(crossing_pts) > 0:
buffers = crossing_pts.copy()
buffers["geometry"] = buffers["geometry"].buffer(50.0)
buffers.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
print(f" ✓ Layer 'crossing_buffers_50m' saved ({len(buffers)} buffers)")
# Clip each buffer to the ward whose boundary it crossed
clipped_parts = []
for bidx, buf in buffers.iterrows():
ward_id = buf["ward_id"]
ward_geom = wards_metric.loc[wards_metric["ward_id"] == ward_id, "geometry"].iloc[0]
clipped = buf["geometry"].intersection(ward_geom)
if clipped.is_empty:
continue
if clipped.type == "MultiPolygon":
for poly in clipped.geoms:
clipped_parts.append({
"geometry": poly,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
else:
clipped_parts.append({
"geometry": clipped,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
if clipped_parts:
clipped_gdf = gpd.GeoDataFrame(clipped_parts, crs=METRIC_CRS)
clipped_gdf["geometry"] = clipped_gdf["geometry"].apply(
lambda g: make_valid(g) if not g.is_valid else g
)
clipped_gdf = clipped_gdf[~clipped_gdf["geometry"].is_empty]
clipped_gdf.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(f" ✓ Layer 'buffer_ward_intersection' saved ({len(clipped_gdf)} polygons)")
else:
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layer 'buffer_ward_intersection' saved (empty)")
else:
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layers 'crossing_buffers_50m' and 'buffer_ward_intersection' saved (empty)")
# ── Step 5: Ward crossing density table (WGS84) ────────────────────────────
print("\n[5/5] Computing crossing density per ward…")
if len(crossing_pts) > 0:
counts = crossing_pts.groupby("ward_id").size().reset_index(name="crossing_count")
else:
counts = pd.DataFrame({"ward_id": wards_gdf["ward_id"], "crossing_count": 0})
# Area from metric CRS
wards_metric["ward_area_km2"] = (wards_metric.geometry.area / 1_000_000.0).round(4)
density_gdf = wards_gdf[["ward_id", "ward_name_en", "ward_name", "geometry"]].copy()
density_gdf = density_gdf.merge(counts, on="ward_id", how="left")
density_gdf["crossing_count"] = density_gdf["crossing_count"].fillna(0).astype(int)
area_map = dict(zip(wards_metric["ward_id"], wards_metric["ward_area_km2"]))
density_gdf["ward_area_km2"] = density_gdf["ward_id"].map(area_map)
density_gdf["crossings_per_km2"] = np.where(
density_gdf["ward_area_km2"] > 0,
(density_gdf["crossing_count"] / density_gdf["ward_area_km2"]).round(4),
0.0
)
# Already in WGS84
density_gdf.to_file(GPKG_PATH, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f" ✓ Layer 'ward_crossing_density_wgs84' saved ({len(density_gdf)} wards)")
# ── Summary ─────────────────────────────────────────────────────────────────
print("\n" + "=" * 60)
print("SUMMARY")
print("=" * 60)
print(f" Wards: {len(wards_gdf)}")
print(f" Crossing points: {len(crossing_pts)}")
print(f" Buffers (50 m): {len(crossing_pts)}")
print(f" Density table: {len(density_gdf)} rows")
print(f"\nOutput: {GPKG_PATH}")
print("Done.")Wrote /work/solve.py (11931 chars).
19:06
Let me run the script now. It may take a while due to the large OSM data download.
19:06
Bash
cd /work && python solve.py 2>&1
============================================================
Tokyo 23 Special Wards — Road Crossing Analysis
============================================================
[1/5] Fetching ward boundaries from OSM…
Retrieved 23 wards
✓ Layer 'wards_jgd' saved
[2/5] Fetching drivable road network from OSM…
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/urllib3/connection.py", line 204, in _new_conn
sock = connection.create_connection(
(self._dns_host, self.port),
...<2 lines>...
socket_options=self.socket_options,
)
File "/usr/local/lib/python3.14/site-packages/urllib3/util/connection.py", line 85, in create_connection
raise err
File "/usr/local/lib/python3.14/site-packages/urllib3/util/connection.py", line 73, in create_connection
sock.connect(sa)
~~~~~~~~~~~~^^^^
ConnectionRefusedError: [Errno 111] Connection refused
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 788, in urlopen
response = self._make_request(
conn,
...<10 lines>...
**response_kw,
)
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 488, in _make_request
raise new_e
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 464, in _make_request
self._validate_conn(conn)
~~~~~~~~~~~~~~~~~~~^^^^^^
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 1106, in _validate_conn
conn.connect()
~~~~~~~~~~~~^^
File "/usr/local/lib/python3.14/site-packages/urllib3/connection.py", line 759, in connect
self.sock = sock = self._new_conn()
~~~~~~~~~~~~~~^^
File "/usr/local/lib/python3.14/site-packages/urllib3/connection.py", line 219, in _new_conn
raise NewConnectionError(
self, f"Failed to establish a new connection: {e}"
) from e
urllib3.exceptions.NewConnectionError: HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/usr/local/lib/python3.14/site-packages/requests/adapters.py", line 696, in send
resp = conn.urlopen(
method=request.method,
...<9 lines>...
chunked=chunked,
)
File "/usr/local/lib/python3.14/site-packages/urllib3/connectionpool.py", line 842, in urlopen
retries = retries.increment(
method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
)
File "/usr/local/lib/python3.14/site-packages/urllib3/util/retry.py", line 543, in increment
raise MaxRetryError(_pool, url, reason) from reason # type: ignore[arg-type]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
urllib3.exceptions.MaxRetryError: HTTPSConnectionPool(host='overpass-api.de', port=443): Max retries exceeded with url: /api/interpreter (Caused by NewConnectionError("HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused"))
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/work/solve.py", line 136, in <module>
G = ox.graph_from_polygon(
combined_bounds_wgs84,
...<3 lines>...
truncate_by_edge=False
)
File "/usr/local/lib/python3.14/site-packages/osmnx/graph.py", line 496, in graph_from_polygon
G_buff = _create_graph(response_jsons, bidirectional)
File "/usr/local/lib/python3.14/site-packages/osmnx/graph.py", line 625, in _create_graph
for response_json in response_jsons:
^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/osmnx/_overpass.py", line 402, in _download_overpass_network
yield _overpass_request(OrderedDict(data=query_str))
~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/osmnx/_overpass.py", line 469, in _overpass_request
response = requests.post(
url,
...<3 lines>...
**settings.requests_kwargs,
)
File "/usr/local/lib/python3.14/site-packages/requests/api.py", line 134, in post
return request("post", url, data=data, json=json, **kwargs)
File "/usr/local/lib/python3.14/site-packages/requests/api.py", line 71, in request
return session.request(method=method, url=url, **kwargs)
~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/requests/sessions.py", line 651, in request
resp = self.send(prep, **send_kwargs)
File "/usr/local/lib/python3.14/site-packages/requests/sessions.py", line 784, in send
r = adapter.send(request, **kwargs)
File "/usr/local/lib/python3.14/site-packages/requests/adapters.py", line 729, in send
raise ConnectionError(e, request=request)
requests.exceptions.ConnectionError: HTTPSConnectionPool(host='overpass-api.de', port=443): Max retries exceeded with url: /api/interpreter (Caused by NewConnectionError("HTTPSConnection(host='overpass-api.de', port=443): Failed to establish a new connection: [Errno 111] Connection refused"))
20:27
The OSMnx `graph_from_polygon` is having connection issues with overpass-api.de. The connection works when we use our own requests with the right User-Agent, but OSMnx is still using the default overpass URL (overpass-api.de) and the default user-agent causes rate-limiting/connection-refused issues. Wait, actually the error is "Connection refused" (Errno 111), which means the port is actively refusing the connection. But earlier we successfully got data from overpass-api.de with the right user-agent. This might be a rate limiting issue - OSMnx is making too many requests too fast. Let me look at what OSMnx is doing - it subdivides the area into multiple queries and sends them one by one. The large polygon area (the whole 23 wards) triggers subdivision, leading to potentially many requests. The overpass-api.de might be rate-limiting us. Let me try using a different approach - download the road network via our own Overpass query directly (bypassing OSMnx for the road download), or use the `custom_filter` approach. Actually, let me try setting the overpass URL to use the `kumi.systems` mirror or `openstreetmap.fr`. But those were also having issues. Alternatively, I can download the entire road network for the Tokyo area using a single Overpass query and parse the results manually. Let me write a custom Overpass query that gets all highways (excluding service roads) within the combined bounding box of the 23 wards. Actually, the simplest approach might be to use `ox.graph_from_bbox` instead of `graph_from_polygon`, or use `features_from_bbox` with the right tags. Let me try a different approach. Wait, actually let me check - the user_agent setting was set correctly in the script, but OSMnx's _overpass.py creates its own headers via _get_http_headers() which uses settings.http_user_agent. So it should be using 'OSM-GIS-Analysis/1.0'. But we're still getting connection refused... Let me try using the OSMnx with a simpler approach - maybe the polygon subdivision is causing too many connections. Let me use `features_from_bbox` to get highways directly instead.
20:27
Bash
python3 -c "
python3 -c "
import osmnx as ox
ox.settings.http_user_agent = 'OSM-GIS-Analysis/1.0'
ox.settings.use_cache = True
# Just a small area - see if it works now
bbox = (139.76, 35.69, 139.78, 35.71)
G = ox.graph_from_bbox(bbox, network_type='drive')
print(f'Nodes: {len(G.nodes)}, Edges: {len(G.edges)}')
" 2>&1
Nodes: 1296, Edges: 2660
21:10
The small area works! The issue is that the full Tokyo 23 wards polygon is very large, causing OSMnx to split it into many sub-queries, and the overpass server is rate-limiting or refusing the connection. Let me try a different approach: 1. Use a direct Overpass query for all roads in the Tokyo 23 wards area, instead of using OSMnx's graph_from_polygon which issues many sub-queries. 2. Or use the direct Overpass query approach - query all ways with highway=motorway|trunk|primary|secondary|tertiary|unclassified|residential (excluding service) in the bounding box. Let me write a direct Overpass query to get drivable roads (excluding service) for the Tokyo area, parse the results into a GeoDataFrame of lines, then proceed.
21:10
21:10
Bash
python3 -c "
python3 -c "
# Get the bbox of the 23 wards
import requests, json
overpass_url = 'https://overpass-api.de/api/interpreter'
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
rel_ids = [1761742, 1758897, 1761717, 1758858, 1758878, 1758888, 1758891, 3554015, 3554304, 1758936, 1758947, 1759474, 1759477, 1543056, 1543055, 1759506, 1760038, 1760040, 1760078, 1760119, 1760124, 1761718, 1761743]
rel_id_str = ','.join(str(i) for i in rel_ids)
# Get just the bbox
query = f'[out:json][timeout:30];rel(id:{rel_id_str});out bb;'
resp = requests.post(overpass_url, data={'data': query}, timeout=60, headers=headers)
data = resp.json()
lats, lons = [], []
for el in data.get('elements', []):
if 'bounds' in el:
b = el['bounds']
lats.extend([b['minlat'], b['maxlat']])
lons.extend([b['minlon'], b['maxlon']])
print(f'bbox: south={min(lats)}, north={max(lats)}, west={min(lons)}, east={max(lons)}')
" 2>&1
bbox: south=35.4816556, north=35.8174937, west=139.5628986, east=139.9189004
21:36
The area is too large for OSMnx to handle in a single query due to its max_query_area_size setting. I should use a direct Overpass query to fetch the drivable roads (highway=*) within the combined boundary, excluding service roads. Let me use the Overpass API directly instead of OSMnx for this step. The Overpass query should fetch all ways with highway tags that are drivable (not service roads), within the bounding box of Tokyo's 23 wards. Let me write the Overpass query.
21:36
Write
/work/solve.py
content · 357 lines · py
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
import requests
import warnings
warnings.filterwarnings('ignore')
ox.settings.use_cache = True
ox.settings.log_console = False
ox.settings.timeout = 300
ox.settings.http_user_agent = "OSM-GIS-Analysis/1.0"
METRIC_CRS = "EPSG:6679" # JGD2011 / Japan Plane Rectangular CS IX (Tokyo)
WGS84_CRS = "EPSG:4326"
GPKG_PATH = "/work/tokyo_crossings.gpkg"
print("=" * 60)
print("Tokyo 23 Special Wards — Road Crossing Analysis")
print("=" * 60)
# ── Step 1: Fetch the 23 special-ward boundaries from Overpass ─────────────
print("\n[1/5] Fetching ward boundaries from OSM…")
# Correct OSM relation IDs for Tokyo's 23 special wards (admin_level=7)
WARD_INFO = {
1761742: ("Chiyoda Ward", "千代田区"),
1758897: ("Chuo Ward", "中央区"),
1761717: ("Minato Ward", "港区"),
1758858: ("Shinjuku Ward", "新宿区"),
1758878: ("Bunkyo Ward", "文京区"),
1758888: ("Taito Ward", "台東区"),
1758891: ("Sumida Ward", "墨田区"),
3554015: ("Koto Ward", "江東区"),
3554304: ("Shinagawa Ward", "品川区"),
1758936: ("Meguro Ward", "目黒区"),
1758947: ("Ota Ward", "大田区"),
1759474: ("Setagaya Ward", "世田谷区"),
1759477: ("Shibuya Ward", "渋谷区"),
1543056: ("Nakano Ward", "中野区"),
1543055: ("Suginami Ward", "杉並区"),
1759506: ("Toshima Ward", "豊島区"),
1760038: ("Kita Ward", "北区"),
1760040: ("Arakawa Ward", "荒川区"),
1760078: ("Itabashi Ward", "板橋区"),
1760119: ("Nerima Ward", "練馬区"),
1760124: ("Adachi Ward", "足立区"),
1761718: ("Katsushika Ward","葛飾区"),
1761743: ("Edogawa Ward", "江戸川区"),
}
rel_ids = list(WARD_INFO.keys())
headers_ok = {"User-Agent": "OSM-GIS-Analysis/1.0"}
overpass_url = "https://overpass-api.de/api/interpreter"
rel_id_str = ",".join(str(i) for i in rel_ids)
query_wards = f"""[out:json][timeout:300];
rel(id:{rel_id_str});
out body geom;
"""
resp = requests.post(overpass_url, data={"data": query_wards}, timeout=300, headers=headers_ok)
resp.raise_for_status()
data = resp.json()
from shapely.geometry import Polygon as ShapelyPolygon, MultiPolygon, LineString, MultiLineString
from shapely.validation import make_valid
ward_records = []
for el in data.get("elements", []):
if el["type"] != "relation":
continue
rid = el["id"]
name_en, name_ja = WARD_INFO[rid]
members = el.get("members", [])
outer_rings = []
inner_rings = []
for m in members:
if m.get("role") == "outer" and "geometry" in m:
coords = [(p["lon"], p["lat"]) for p in m["geometry"]]
if len(coords) >= 3:
outer_rings.append(coords)
elif m.get("role") == "inner" and "geometry" in m:
coords = [(p["lon"], p["lat"]) for p in m["geometry"]]
if len(coords) >= 3:
inner_rings.append(coords)
if not outer_rings:
print(f" ⚠ No outer rings for relation {rid}, skipping")
continue
if len(outer_rings) == 1:
poly = ShapelyPolygon(outer_rings[0], inner_rings if inner_rings else None)
else:
polys = []
remaining_inner = list(inner_rings)
for o_ring in outer_rings:
o_poly = ShapelyPolygon(o_ring)
holes_for_this = []
keep = []
for i_ring in remaining_inner:
i_poly = ShapelyPolygon(i_ring)
if o_poly.contains(i_poly.representative_point()):
holes_for_this.append(i_ring)
else:
keep.append(i_ring)
remaining_inner = keep
polys.append(ShapelyPolygon(o_ring, holes_for_this if holes_for_this else None))
poly = MultiPolygon(polys) if len(polys) > 1 else polys[0]
if not poly.is_valid:
poly = make_valid(poly)
ward_records.append({
"ward_id": rid,
"ward_name_en": name_en,
"ward_name": name_ja,
"geometry": poly
})
wards_gdf = gpd.GeoDataFrame(ward_records, crs=WGS84_CRS)
wards_metric = wards_gdf.to_crs(METRIC_CRS)
wards_metric.to_file(GPKG_PATH, layer="wards_jgd", driver="GPKG")
print(f" ✓ Layer 'wards_jgd' saved ({len(wards_gdf)} wards)")
# ── Step 2: Get drivable road network via direct Overpass query ────────────
print("\n[2/5] Fetching drivable road network from OSM…")
# Build a bbox for the area
all_lats = [pt.y for g in wards_gdf.geometry for pt in g.exterior.coords]
all_lons = [pt.x for g in wards_gdf.geometry for pt in g.exterior.coords]
south, north = min(all_lats), max(all_lats)
west, east = min(all_lons), max(all_lons)
print(f" Bounding box: {south:.4f}, {north:.4f}, {west:.4f}, {east:.4f}")
# Drivable highway tags (excluding service, footway, cycleway, path, etc.)
# We use the 'highway' key with values matching drivable roads.
# network_type='drive' in OSMnx includes: motorway, trunk, primary, secondary,
# tertiary, unclassified, residential, motorway_link, trunk_link, primary_link,
# secondary_link, tertiary_link, living_street
# Excludes: service, footway, cycleway, path, track, pedestrian, steps, etc.
# Overpass query for drivable roads within the bbox
query_roads = f"""[out:json][timeout:600];
(
way["highway"~"^(motorway|trunk|primary|secondary|tertiary|unclassified|residential|motorway_link|trunk_link|primary_link|secondary_link|tertiary_link|living_street)$"]({south},{west},{north},{east});
);
out geom;
"""
print(" Sending Overpass query for drivable roads…")
resp = requests.post(overpass_url, data={"data": query_roads}, timeout=600, headers=headers_ok)
resp.raise_for_status()
roads_data = resp.json()
print(f" Received {len(roads_data.get('elements', []))} road elements")
# Build road geometries
road_geoms = []
for el in roads_data.get("elements", []):
if el["type"] != "way":
continue
if "geometry" not in el:
continue
coords = [(p["lon"], p["lat"]) for p in el["geometry"]]
if len(coords) < 2:
continue
road_geoms.append(LineString(coords))
if not road_geoms:
print(" WARNING: No road geometries found!")
roads_gdf = gpd.GeoDataFrame({"geometry": gpd.GeoSeries([], crs=WGS84_CRS)})
else:
roads_gdf = gpd.GeoDataFrame({"geometry": road_geoms}, crs=WGS84_CRS)
print(f" Road segments loaded: {len(roads_gdf)}")
# Project to metric and clip to combined boundary
roads_metric = roads_gdf.to_crs(METRIC_CRS)
combined_bounds = wards_metric.unary_union
roads_metric = roads_metric[roads_metric.intersects(combined_bounds)].copy()
roads_metric["geometry"] = roads_metric["geometry"].apply(
lambda g: g.intersection(combined_bounds) if g is not None else g
)
roads_metric = roads_metric[~roads_metric["geometry"].is_empty & roads_metric["geometry"].is_valid]
roads_metric = roads_metric[roads_metric.geometry.type.isin(["LineString", "MultiLineString"])]
print(f" Roads after clipping: {len(roads_metric)}")
# Flatten MultiLineStrings
simple_edges = []
for idx, row in roads_metric.iterrows():
geom = row.geometry
if geom.type == "MultiLineString":
for line in geom.geoms:
simple_edges.append(line)
else:
simple_edges.append(geom)
roads_final = gpd.GeoDataFrame({"geometry": simple_edges}, crs=METRIC_CRS)
print(f" Simplified to {len(roads_final)} single-part line segments")
# ── Step 3: Find crossing points ───────────────────────────────────────────
print("\n[3/5] Detecting road / ward-boundary crossings…")
all_crossings = []
for idx, ward in wards_metric.iterrows():
ward_id = ward["ward_id"]
ward_name_en = ward["ward_name_en"]
ward_name = ward["ward_name"]
ward_boundary = ward.geometry.boundary
candidate_mask = roads_final.intersects(ward_boundary)
candidate_roads = roads_final[candidate_mask]
for ridx, road in candidate_roads.iterrows():
road_geom = road.geometry
inter = road_geom.intersection(ward_boundary)
if inter.is_empty:
continue
if inter.type == "Point":
all_crossings.append({
"geometry": inter,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
elif inter.type == "MultiPoint":
for pt in inter.geoms:
all_crossings.append({
"geometry": pt,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
# Skip LineString/MultiLineString (road runs along boundary)
print(f" Found {len(all_crossings)} crossing points")
if len(all_crossings) > 0:
crossing_pts = gpd.GeoDataFrame(all_crossings, crs=METRIC_CRS)
else:
crossing_pts = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
crossing_pts.to_file(GPKG_PATH, layer="crossing_points", driver="GPKG")
print(f" ✓ Layer 'crossing_points' saved ({len(crossing_pts)} points)")
# ── Step 4: Create 50 m buffers ────────────────────────────────────────────
print("\n[4/5] Creating 50 m buffers and clipping to wards…")
if len(crossing_pts) > 0:
buffers = crossing_pts.copy()
buffers["geometry"] = buffers["geometry"].buffer(50.0)
buffers.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
print(f" ✓ Layer 'crossing_buffers_50m' saved ({len(buffers)} buffers)")
# Clip each buffer to the ward whose boundary it crossed
clipped_parts = []
for bidx, buf in buffers.iterrows():
ward_id = buf["ward_id"]
ward_geom = wards_metric.loc[wards_metric["ward_id"] == ward_id, "geometry"].iloc[0]
clipped = buf["geometry"].intersection(ward_geom)
if clipped.is_empty:
continue
if clipped.type == "MultiPolygon":
for poly in clipped.geoms:
clipped_parts.append({
"geometry": poly,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
else:
clipped_parts.append({
"geometry": clipped,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
if clipped_parts:
clipped_gdf = gpd.GeoDataFrame(clipped_parts, crs=METRIC_CRS)
clipped_gdf["geometry"] = clipped_gdf["geometry"].apply(
lambda g: make_valid(g) if not g.is_valid else g
)
clipped_gdf = clipped_gdf[~clipped_gdf["geometry"].is_empty]
clipped_gdf.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(f" ✓ Layer 'buffer_ward_intersection' saved ({len(clipped_gdf)} polygons)")
else:
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layer 'buffer_ward_intersection' saved (empty)")
else:
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layers 'crossing_buffers_50m' and 'buffer_ward_intersection' saved (empty)")
# ── Step 5: Ward crossing density table (WGS84) ────────────────────────────
print("\n[5/5] Computing crossing density per ward…")
if len(crossing_pts) > 0:
counts = crossing_pts.groupby("ward_id").size().reset_index(name="crossing_count")
else:
counts = pd.DataFrame({"ward_id": wards_gdf["ward_id"], "crossing_count": 0})
wards_metric["ward_area_km2"] = (wards_metric.geometry.area / 1_000_000.0).round(4)
density_gdf = wards_gdf[["ward_id", "ward_name_en", "ward_name", "geometry"]].copy()
density_gdf = density_gdf.merge(counts, on="ward_id", how="left")
density_gdf["crossing_count"] = density_gdf["crossing_count"].fillna(0).astype(int)
area_map = dict(zip(wards_metric["ward_id"], wards_metric["ward_area_km2"]))
density_gdf["ward_area_km2"] = density_gdf["ward_id"].map(area_map)
density_gdf["crossings_per_km2"] = np.where(
density_gdf["ward_area_km2"] > 0,
(density_gdf["crossing_count"] / density_gdf["ward_area_km2"]).round(4),
0.0
)
density_gdf.to_file(GPKG_PATH, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f" ✓ Layer 'ward_crossing_density_wgs84' saved ({len(density_gdf)} wards)")
# ── Verify & Summary ───────────────────────────────────────────────────────
print("\n" + "=" * 60)
print("SUMMARY")
print("=" * 60)
print(f" Wards: {len(wards_gdf)}")
print(f" Crossing points: {len(crossing_pts)}")
print(f" Buffers (50 m): {len(crossing_pts)}")
print(f" Density table: {len(density_gdf)} rows")
print(f"\nOutput: {GPKG_PATH}")
# Verify layers
print("\n Layer verification:")
import fiona
layers = fiona.listlayers(GPKG_PATH)
for lyr in layers:
gdf = gpd.read_file(GPKG_PATH, layer=lyr)
print(f" {lyr}: {len(gdf)} features, CRS={gdf.crs}")
print("\nDone.")Wrote /work/solve.py (13582 chars).
21:43
Bash
cd /work && python solve.py 2>&1
============================================================
Tokyo 23 Special Wards — Road Crossing Analysis
============================================================
[1/5] Fetching ward boundaries from OSM…
✓ Layer 'wards_jgd' saved (23 wards)
[2/5] Fetching drivable road network from OSM…
Traceback (most recent call last):
File "/work/solve.py", line 132, in <module>
all_lats = [pt.y for g in wards_gdf.geometry for pt in g.exterior.coords]
^^^^^^^^^^
AttributeError: 'MultiPolygon' object has no attribute 'exterior'
22:04
The error is because some ward geometries are MultiPolygons and don't have .exterior. Let me fix the bbox extraction.
22:04
Edit
/work/solve.py
Edited /work/solve.py.
22:20
Bash
cd /work && python solve.py 2>&1
============================================================
Tokyo 23 Special Wards — Road Crossing Analysis
============================================================
[1/5] Fetching ward boundaries from OSM…
✓ Layer 'wards_jgd' saved (23 wards)
[2/5] Fetching drivable road network from OSM…
Bounding box: 35.4817, 35.8175, 139.5629, 139.9189
Sending Overpass query for drivable roads…
Received 127738 road elements
Road segments loaded: 127738
Roads after clipping: 4338
Simplified to 4634 single-part line segments
[3/5] Detecting road / ward-boundary crossings…
Found 5462 crossing points
✓ Layer 'crossing_points' saved (5462 points)
[4/5] Creating 50 m buffers and clipping to wards…
✓ Layer 'crossing_buffers_50m' saved (5462 buffers)
✓ Layer 'buffer_ward_intersection' saved (9725 polygons)
[5/5] Computing crossing density per ward…
✓ Layer 'ward_crossing_density_wgs84' saved (23 wards)
============================================================
SUMMARY
============================================================
Wards: 23
Crossing points: 5462
Buffers (50 m): 5462
Density table: 23 rows
Output: /work/tokyo_crossings.gpkg
Layer verification:
wards_jgd: 23 features, CRS=EPSG:6679
crossing_points: 5462 features, CRS=EPSG:6679
crossing_buffers_50m: 5462 features, CRS=EPSG:6679
buffer_ward_intersection: 9725 features, CRS=EPSG:6679
ward_crossing_density_wgs84: 23 features, CRS=EPSG:4326
Done.
23:06
The script ran successfully. Let me verify the output and check a few things: 1. All 5 layers are present 2. The first 4 layers are in metric CRS (EPSG:6679) ✓ 3. The density layer is in WGS84 (EPSG:4326) ✓ 4. All layers have ward_id, ward_name_en, ward_name where appropriate Let me do a quick sanity check on the data.
23:06
23:06
Bash
python3 -c "
python3 -c "
import geopandas as gpd
import pandas as pd
GPKG = '/work/tokyo_crossings.gpkg'
layers = ['wards_jgd', 'crossing_points', 'crossing_buffers_50m', 'buffer_ward_intersection', 'ward_crossing_density_wgs84']
for lyr in layers:
gdf = gpd.read_file(GPKG, layer=lyr)
print(f'=== {lyr} ===')
print(f' Features: {len(gdf)}')
print(f' CRS: {gdf.crs}')
print(f' Columns: {gdf.columns.tolist()}')
if len(gdf) > 0:
print(f' Sample: {gdf.iloc[0].to_dict()}')
print()
# Show density table
dens = gpd.read_file(GPKG, layer='ward_crossing_density_wgs84')
pd.set_option('display.max_rows', 30)
pd.set_option('display.width', 120)
print('=== Density table ===')
print(dens[['ward_id', 'ward_name_en', 'crossing_count', 'ward_area_km2', 'crossings_per_km2']].to_string(index=False))
" 2>&1
=== wards_jgd ===
Features: 23
CRS: EPSG:6679
Columns: ['ward_id', 'ward_name_en', 'ward_name', 'geometry']
Sample: {'ward_id': 1543055, 'ward_name_en': 'Suginami Ward', 'ward_name': '杉並区', 'geometry': <MULTIPOLYGON (((-59218.459 -923246.81, -59200.631 -923248.77, -59185.058 -9...>}
=== crossing_points ===
Features: 5462
CRS: EPSG:6679
Columns: ['ward_id', 'ward_name_en', 'ward_name', 'geometry']
Sample: {'ward_id': 1543055, 'ward_name_en': 'Suginami Ward', 'ward_name': '杉並区', 'geometry': <POINT (-59809.471 -922079.375)>}
=== crossing_buffers_50m ===
Features: 5462
CRS: EPSG:6679
Columns: ['ward_id', 'ward_name_en', 'ward_name', 'geometry']
Sample: {'ward_id': 1543055, 'ward_name_en': 'Suginami Ward', 'ward_name': '杉並区', 'geometry': <POLYGON ((-59759.471 -922079.375, -59759.711 -922084.275, -59760.431 -92208...>}
=== buffer_ward_intersection ===
Features: 9725
CRS: EPSG:6679
Columns: ['ward_id', 'ward_name_en', 'ward_name', 'geometry']
Sample: {'ward_id': 1543055, 'ward_name_en': 'Suginami Ward', 'ward_name': '杉並区', 'geometry': <POLYGON ((-59765.375 -922102.944, -59767.897 -922107.153, -59769.016 -92210...>}
=== ward_crossing_density_wgs84 ===
Features: 23
CRS: EPSG:4326
Columns: ['ward_id', 'ward_name_en', 'ward_name', 'crossing_count', 'ward_area_km2', 'crossings_per_km2', 'geometry']
Sample: {'ward_id': 1543055, 'ward_name_en': 'Suginami Ward', 'ward_name': '杉並区', 'crossing_count': 335, 'ward_area_km2': 0.8584, 'crossings_per_km2': 390.261, 'geometry': <MULTIPOLYGON (((139.596 35.682, 139.596 35.682, 139.596 35.682, 139.596 35....>}
=== Density table ===
ward_id ward_name_en crossing_count ward_area_km2 crossings_per_km2
1543055 Suginami Ward 335 0.8584 390.2610
1543056 Nakano Ward 331 0.6435 514.3745
1758858 Shinjuku Ward 277 0.6461 428.7262
1758878 Bunkyo Ward 228 0.5948 383.3221
1758888 Taito Ward 145 0.2607 556.1949
1758891 Sumida Ward 75 0.2591 289.4635
1758897 Chuo Ward 98 0.1218 804.5977
1758936 Meguro Ward 279 0.4410 632.6531
1758947 Ota Ward 223 2.6255 84.9362
1759474 Setagaya Ward 447 1.3884 321.9533
1759477 Shibuya Ward 247 0.4343 568.7313
1759506 Toshima Ward 311 0.7039 441.8241
1760038 Kita Ward 191 0.7940 240.5542
1760040 Arakawa Ward 141 0.3489 404.1273
1760078 Itabashi Ward 289 0.9099 317.6173
1760119 Nerima Ward 490 1.3772 355.7944
1760124 Adachi Ward 240 1.5642 153.4331
1761717 Minato Ward 219 0.6419 341.1746
1761718 Katsushika Ward 242 2.1543 112.3335
1761742 Chiyoda Ward 176 0.1431 1229.9092
1761743 Edogawa Ward 152 8.1724 18.5992
3554015 Koto Ward 94 6.0708 15.4840
3554304 Shinagawa Ward 232 1.0928 212.2987
23:20
The results look good. Let me verify all the requirements: 1. ✅ `wards_jgd` (Polygon) - 23 ward boundaries in EPSG:6679 (metric) 2. ✅ `crossing_points` (Point) - 5462 points at road/ward-boundary crossings 3. ✅ `crossing_buffers_50m` (Polygon) - 5462 buffers of 50m each 4. ✅ `buffer_ward_intersection` (Polygon) - 9725 buffer-clipped-to-ward polygons 5. ✅ `ward_crossing_density_wgs84` (Polygon) - 23 wards with ward_id, crossing_count, ward_area_km2, crossings_per_km2 in WGS84 All layers have `ward_id`, `ward_name_en`, and `ward_name` (except the density table which also needs them - wait, let me check: "carry the English and native ward names as ward_name_en and ward_name on every layer that's tied to a specific ward" - all layers have these columns except... let me check wards_jgd has them, crossing_points has them, crossing_buffers_50m has them, buffer_ward_intersection has them, and ward_crossing_density_wgs84 has them. Yes, all good. The area values look a bit odd - Suginami is 0.8584 km²? That seems very small for a ward. Wait, let me check - Suginami Ward is actually about 34 km². Something is wrong with the area calculation. Let me check the CRS more carefully. EPSG:6679 is JGD2011 / Japan Plane Rectangular CS IX. This uses meters as the unit. Let me check the coordinates we're getting. Actually, looking at the coordinates in the sample: `-59218.459 -923246.81` - these are large negative numbers which seems wrong for Tokyo. The coordinates should be in meters easting/northing, and Tokyo in CS IX should have easting around -61000 to -45000 m and northing around -... wait, this seems off. Let me check: The JGD2011 / Japan Plane Rectangular CS IX (EPSG:6679) has an origin at approximately 36°N, 139°50'E. The false easting is 0 and false northing is 0. So Tokyo would have coordinates around -60,000 m easting and... the northing would be around... Let me verify. Actually, looking more carefully at the coordinates: -59218 easting and -923247 northing. That seems very far south. Let me check the actual metrics. Wait - the area shows 0.8584 km² for Suginami. Suginami is actually about 34 km². So something is off. Let me check which CRS we should be using. Actually EPSG:6679 is JGD2011 / Japan Plane Rectangular CS IX. The coordinate system covers the Tokyo area with coordinates in meters, but the origin is offset from the area. The area calculation should still be correct in meters squared, regardless of the coordinate values. Let me verify by checking the actual area manually.
23:20
23:20
Bash
python3 -c "
python3 -c "
import geopandas as gpd
import numpy as np
GPKG = '/work/tokyo_crossings.gpkg'
# Check CRS details
wards = gpd.read_file(GPKG, layer='wards_jgd')
print('CRS:', wards.crs)
print('CRS axis info:', wards.crs.axis_info)
# Check Suginami Ward
sug = wards[wards['ward_name_en']=='Suginami Ward']
geom = sug.geometry.iloc[0]
print(f'Suginami area from metric: {geom.area / 1e6:.4f} km²')
# Check what area we get from WGS84
wards_wgs84 = wards.to_crs('EPSG:4326')
sug_wgs = wards_wgs84[wards_wgs84['ward_name_en']=='Suginami Ward']
print(f'Suginami area from WGS84: {sug_wgs.geometry.area:.4f} deg²')
# Use an equal-area projection to check
wards_equal = wards.to_crs('EPSG:6933') # World Equidistant Cylindrical... actually use Mollweide
try:
wards_equal = wards.to_crs('ESRI:54009') # Mollweide
sug_eq = wards_equal[wards_equal['ward_name_en']=='Suginami Ward']
print(f'Suginami area from equal-area: {sug_eq.geometry.area.iloc[0] / 1e6:.4f} km²')
except:
pass
# Check a known area
# Suginami ward is about 34 km²
# Check if EPSG:6679 uses meters or something else
from pyproj import CRS
crs = CRS('EPSG:6679')
print(f'CRS unit: {crs.axis_info}')
print(f'CRS area of use: {crs.area_of_use}')
" 2>&1
CRS: EPSG:6679
CRS axis info: [Axis(name=Northing, abbrev=X, direction=north, unit_auth_code=EPSG, unit_code=9001, unit_name=metre), Axis(name=Easting, abbrev=Y, direction=east, unit_auth_code=EPSG, unit_code=9001, unit_name=metre)]
Suginami area from metric: 0.8584 km²
<string>:19: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
Traceback (most recent call last):
File "<string>", line 19, in <module>
print(f'Suginami area from WGS84: {sug_wgs.geometry.area:.4f} deg²')
^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: unsupported format string passed to Series.__format__
23:32
The area seems too small. Suginami ward should be about 34 km² but we're getting 0.8584 km². Something's wrong with the CRS. Let me check the CRS more carefully. EPSG:6679 is "JGD2011 / Japan Plane Rectangular CS IX" which covers Tokyo area. The units should be meters, but the coordinates are very large negative numbers like (-59218.459, -923246.81). Wait, the axis info shows Northing first (X=northing, Y=easting)? That seems unusual. Actually, looking at the coordinates: (-59218.459, -923246.81) - these are large numbers. But the area being 0.858 km² for Suginami is way too small. Suginami is about 34 km². Wait, let me check: if the coordinates are in meters, and the polygon coordinates are like (-59218, -923247), the area would be computed using these coordinates. The area of a polygon in a planar CRS is in the units of the CRS squared. So for meters, it should be m². 0.858 km² seems too small for Suginami. But actually looking more carefully: (-59218, -923247) - these are coordinates relative to the origin of the Japan Plane Rectangular CS IX. The origin is at 36°N, 139°50'E (approximately Tokyo). So Suginami should be near this origin. Let me check the actual extent of the coordinates. If Suginami is about 34 km² = 34,000,000 m², and the coordinates are in meters, the bounding box would be roughly 5.8 km × 5.8 km. Actually, the issue might be related to the axis order. In EPSG:6679, the axis order is Northing (X) then Easting (Y), but geopandas might interpret them differently. Let me check. Actually, wait - the area 0.8584 km² = 858,400 m². If Suginami's area is about 34 km², that's 34,000,000 m². So we're off by a factor of about 40. Hmm, but I see the area is 0.8584 km² - which in km² is very small. Let me check OSM area for Suginami ward. Actually Suginami ward (杉並区) has an area of about 34 km². Let me debug this.
23:32
Bash
python3 -c "
python3 -c "
from pyproj import CRS, Transformer
import geopandas as gpd
import numpy as np
# Check the CRS
crs = CRS('EPSG:6679')
print('CRS name:', crs.name)
print('CRS type:', crs.type_name)
print('Axis info:', crs.axis_info)
for a in crs.axis_info:
print(f' {a.name}: abbrev={a.abbrev}, direction={a.direction}, unit={a.unit_name}')
# Check with a known point in Tokyo
# Tokyo Station: 35.6812°N, 139.7671°E
trans = Transformer.from_crs('EPSG:4326', 'EPSG:6679', always_xy=True)
x, y = trans.transform(139.7671, 35.6812)
print(f'Tokyo Station -> EPSG:6679: ({x:.2f}, {y:.2f})')
# What does Suginami ward extent look like?
wards = gpd.read_file('/work/tokyo_crossings.gpkg', layer='wards_jgd')
sug = wards[wards['ward_name_en']=='Suginami Ward']
bounds = sug.geometry.iloc[0].bounds
print(f'Suginami bounds: {bounds}')
print(f'Suginami width: {bounds[2]-bounds[0]:.0f} m, height: {bounds[3]-bounds[1]:.0f} m')
area_m2 = sug.geometry.iloc[0].area
print(f'Area: {area_m2:.0f} m² = {area_m2/1e6:.4f} km²')
# Real area of Suginami Ward should be ~34 km²
# Let's check using an equal area projection for Japan
wards_ea = wards.to_crs('EPSG:3413') # Arctic - no good
# Try Asia Lambert Conformal Conic
wards_eq = wards.to_crs('ESRI:102026') # Asia Lambert Conformal Conic
sug_eq = wards_eq[wards_eq['ward_name_en']=='Suginami Ward']
area_eq = sug_eq.geometry.iloc[0].area
print(f'Area (ESRI:102026): {area_eq/1e6:.4f} km²')
# Try EPSG:2451 (JGD2000 / Japan Plane Rectangular CS IX)
try:
wards_jgd = wards.to_crs('EPSG:2451')
sug_jgd = wards_jgd[wards_jgd['ward_name_en']=='Suginami Ward']
area_jgd = sug_jgd.geometry.iloc[0].area
print(f'Area (EPSG:2451): {area_jgd/1e6:.4f} km²')
except Exception as e:
print(f'EPSG:2451 error: {e}')
" 2>&1
CRS name: JGD2011 / Japan Plane Rectangular CS XI CRS type: Projected CRS Axis info: [Axis(name=Northing, abbrev=X, direction=north, unit_auth_code=EPSG, unit_code=9001, unit_name=metre), Axis(name=Easting, abbrev=Y, direction=east, unit_auth_code=EPSG, unit_code=9001, unit_name=metre)] Northing: abbrev=X, direction=north, unit=metre Easting: abbrev=Y, direction=east, unit=metre Tokyo Station -> EPSG:6679: (-43710.40, -923452.68) Suginami bounds: (-60158.81897706188, -925087.6852513355, -52692.17428875115, -917805.0090825626) Suginami width: 7467 m, height: 7283 m Area: 858357 m² = 0.8584 km² Area (ESRI:102026): 0.7836 km² Area (EPSG:2451): 0.8583 km²
24:07
Wait - EPSG:6679 is "JGD2011 / Japan Plane Rectangular CS XI" - that's Zone XI which covers Tokyo. But the coordinates are showing X (Northing) around -920,000 which is unusual. Let me check - the Japan Plane Rectangular CS systems have false northing/easting. Actually, the issue might be that the axis order is (Northing, Easting) but we're interpreting as (Easting, Northing). Let me check the area more carefully. Suginami Ward's real area is about 34.58 km², but we're getting 0.86 km². So something is definitely wrong. Let me check if it's a unit issue. The axis info says the unit is metre (9001). But maybe the transformation is wrong? Actually, I think the issue might be that EPSG:6679 uses a local projection but the coordinates seem very large (near -925,000 in northing). This is correct for a false northing of -... wait, let me check. Actually, looking at the bounds: x (northing) from -925087 to -917805, and y (easting) from -60158 to -52692. Width = ~7467 m, height = ~7283 m. That's about 7.5 km × 7.3 km = 54.7 km² if it were a rectangle. But the area computed is 0.86 km²... which doesn't make sense. Wait, the area function gives 858,357 m² which is 0.858 km². But the bounding box is ~7.5 km × ~7.3 km which should give ~55 km². Something is very wrong with the area calculation. Let me check: could the geometry be invalid or have some issue? Actually, maybe the MultiPolygon has many small parts? Or the geometry is corrupted? Let me check the actual geometry type and number of parts.
24:07
Bash
python3 -c "
python3 -c "
import geopandas as gpd
from shapely.validation import make_valid
wards = gpd.read_file('/work/tokyo_crossings.gpkg', layer='wards_jgd')
sug = wards[wards['ward_name_en']=='Suginami Ward']
geom = sug.geometry.iloc[0]
print(f'Type: {geom.geom_type}')
print(f'Is valid: {geom.is_valid}')
if geom.geom_type == 'MultiPolygon':
print(f'Num polygons: {len(geom.geoms)}')
for i, p in enumerate(geom.geoms):
print(f' Polygon {i}: area={p.area:.0f} m², exterior len={len(p.exterior.coords)}')
if p.interiors:
for j, hole in enumerate(p.interiors):
print(f' Hole {j}: len={len(hole.coords)}')
else:
print(f'Area: {geom.area:.0f} m²')
print(f'Exterior coords: {len(geom.exterior.coords)}')
" 2>&1
Type: MultiPolygon Is valid: True Num polygons: 166 Polygon 0: area=2689 m², exterior len=18 Polygon 1: area=5 m², exterior len=4 Polygon 2: area=7390 m², exterior len=16 Polygon 3: area=2823 m², exterior len=7 Polygon 4: area=639 m², exterior len=5 Polygon 5: area=935 m², exterior len=7 Polygon 6: area=5871 m², exterior len=10 Polygon 7: area=49220 m², exterior len=41 Polygon 8: area=214 m², exterior len=5 Polygon 9: area=13 m², exterior len=5 Polygon 10: area=536 m², exterior len=6 Polygon 11: area=2932 m², exterior len=11 Polygon 12: area=20268 m², exterior len=16 Polygon 13: area=52 m², exterior len=5 Polygon 14: area=7 m², exterior len=4 Polygon 15: area=220 m², exterior len=9 Polygon 16: area=2093 m², exterior len=13 Polygon 17: area=7 m², exterior len=4 Polygon 18: area=377 m², exterior len=9 Polygon 19: area=28 m², exterior len=6 Polygon 20: area=283 m², exterior len=8 Polygon 21: area=23 m², exterior len=5 Polygon 22: area=7009 m², exterior len=13 Polygon 23: area=198 m², exterior len=5 Polygon 24: area=0 m², exterior len=4 Polygon 25: area=1978 m², exterior len=7 Polygon 26: area=7 m², exterior len=5 Polygon 27: area=4 m², exterior len=4 Polygon 28: area=0 m², exterior len=4 Polygon 29: area=3348 m², exterior len=14 Polygon 30: area=0 m², exterior len=4 Polygon 31: area=15 m², exterior len=4 Polygon 32: area=213488 m², exterior len=51 Polygon 33: area=6489 m², exterior len=13 Polygon 34: area=20086 m², exterior len=32 Polygon 35: area=3812 m², exterior len=11 Polygon 36: area=49 m², exterior len=5 Polygon 37: area=7232 m², exterior len=15 Polygon 38: area=309 m², exterior len=6 Polygon 39: area=1710 m², exterior len=9 Polygon 40: area=934 m², exterior len=10 Polygon 41: area=8058 m², exterior len=11 Polygon 42: area=21649 m², exterior len=28 Polygon 43: area=28160 m², exterior len=36 Polygon 44: area=6 m², exterior len=4 Polygon 45: area=7 m², exterior len=4 Polygon 46: area=255 m², exterior len=5 Polygon 47: area=1736 m², exterior len=7 Polygon 48: area=1069 m², exterior len=9 Polygon 49: area=26 m², exterior len=5 Polygon 50: area=2889 m², exterior len=11 Polygon 51: area=982 m², exterior len=10 Polygon 52: area=1646 m², exterior len=16 Polygon 53: area=92 m², exterior len=5 Polygon 54: area=10170 m², exterior len=20 Polygon 55: area=3253 m², exterior len=11 Polygon 56: area=211 m², exterior len=8 Polygon 57: area=241 m², exterior len=6 Polygon 58: area=3812 m², exterior len=6 Polygon 59: area=20771 m², exterior len=20 Polygon 60: area=2 m², exterior len=4 Polygon 61: area=449 m², exterior len=8 Polygon 62: area=517 m², exterior len=6 Polygon 63: area=112 m², exterior len=7 Polygon 64: area=32 m², exterior len=5 Polygon 65: area=9 m², exterior len=4 Polygon 66: area=979 m², exterior len=10 Polygon 67: area=1695 m², exterior len=7 Polygon 68: area=971 m², exterior len=9 Polygon 69: area=7022 m², exterior len=16 Polygon 70: area=6 m², exterior len=4 Polygon 71: area=85 m², exterior len=4 Polygon 72: area=5379 m², exterior len=8 Polygon 73: area=1134 m², exterior len=12 Polygon 74: area=0 m², exterior len=4 Polygon 75: area=58 m², exterior len=7 Polygon 76: area=168 m², exterior len=7 Polygon 77: area=5 m², exterior len=4 Polygon 78: area=455 m², exterior len=10 Polygon 79: area=22396 m², exterior len=18 Polygon 80: area=69 m², exterior len=8 Polygon 81: area=13 m², exterior len=6 Polygon 82: area=25252 m², exterior len=32 Polygon 83: area=3 m², exterior len=5 Polygon 84: area=93 m², exterior len=5 Polygon 85: area=19 m², exterior len=6 Polygon 86: area=3236 m², exterior len=9 Polygon 87: area=5491 m², exterior len=11 Polygon 88: area=685 m², exterior len=13 Polygon 89: area=2006 m², exterior len=11 Polygon 90: area=17553 m², exterior len=22 Polygon 91: area=2082 m², exterior len=10 Polygon 92: area=474 m², exterior len=10 Polygon 93: area=41 m², exterior len=4 Polygon 94: area=29123 m², exterior len=28 Polygon 95: area=5 m², exterior len=4 Polygon 96: area=2134 m², exterior len=15 Polygon 97: area=6009 m², exterior len=13 Polygon 98: area=411 m², exterior len=9 Polygon 99: area=0 m², exterior len=4 Polygon 100: area=1287 m², exterior len=10 Polygon 101: area=29196 m², exterior len=11 Polygon 102: area=1 m², exterior len=4 Polygon 103: area=38 m², exterior len=6 Polygon 104: area=0 m², exterior len=4 Polygon 105: area=8399 m², exterior len=21 Polygon 106: area=116 m², exterior len=6 Polygon 107: area=140 m², exterior len=7 Polygon 108: area=6006 m², exterior len=18 Polygon 109: area=7 m², exterior len=4 Polygon 110: area=570 m², exterior len=8 Polygon 111: area=1243 m², exterior len=9 Polygon 112: area=8374 m², exterior len=25 Polygon 113: area=13743 m², exterior len=16 Polygon 114: area=216 m², exterior len=5 Polygon 115: area=4590 m², exterior len=17 Polygon 116: area=93 m², exterior len=4 Polygon 117: area=40494 m², exterior len=27 Polygon 118: area=495 m², exterior len=7 Polygon 119: area=1 m², exterior len=4 Polygon 120: area=8785 m², exterior len=13 Polygon 121: area=18 m², exterior len=4 Polygon 122: area=98 m², exterior len=7 Polygon 123: area=98 m², exterior len=7 Polygon 124: area=18 m², exterior len=4 Polygon 125: area=0 m², exterior len=4 Polygon 126: area=101 m², exterior len=7 Polygon 127: area=39 m², exterior len=5 Polygon 128: area=605 m², exterior len=5 Polygon 129: area=241 m², exterior len=5 Polygon 130: area=8097 m², exterior len=20 Polygon 131: area=896 m², exterior len=14 Polygon 132: area=16 m², exterior len=6 Polygon 133: area=23 m², exterior len=4 Polygon 134: area=0 m², exterior len=4 Polygon 135: area=7 m², exterior len=4 Polygon 136: area=0 m², exterior len=4 Polygon 137: area=14 m², exterior len=4 Polygon 138: area=5 m², exterior len=4 Polygon 139: area=65 m², exterior len=4 Polygon 140: area=2129 m², exterior len=15 Polygon 141: area=12350 m², exterior len=34 Polygon 142: area=4730 m², exterior len=16 Polygon 143: area=8989 m², exterior len=17 Polygon 144: area=168 m², exterior len=7 Polygon 145: area=217 m², exterior len=6 Polygon 146: area=1500 m², exterior len=12 Polygon 147: area=4594 m², exterior len=13 Polygon 148: area=2902 m², exterior len=12 Polygon 149: area=1 m², exterior len=5 Polygon 150: area=531 m², exterior len=9 Polygon 151: area=8634 m², exterior len=17 Polygon 152: area=72 m², exterior len=5 Polygon 153: area=76 m², exterior len=5 Polygon 154: area=127 m², exterior len=4 Polygon 155: area=6518 m², exterior len=24 Polygon 156: area=1602 m², exterior len=10 Polygon 157: area=40654 m², exterior len=42 Polygon 158: area=13968 m², exterior len=23 Polygon 159: area=371 m², exterior len=5 Polygon 160: area=5932 m², exterior len=14 Polygon 161: area=7 m², exterior len=4 Polygon 162: area=798 m², exterior len=9 Polygon 163: area=22 m², exterior len=4 Polygon 164: area=1194 m², exterior len=9 Polygon 165: area=1636 m², exterior len=15
24:28
The issue is that the Overpass query returned each "way" member of the relation separately as individual polygons instead of merging them into a single contiguous polygon. That's why Suginami Ward is showing as 166 tiny polygons instead of one big polygon (the actual area of Suginami is about 34 km², not 0.86 km²). The problem is in how I parsed the relation members. The ward boundaries in OSM are made of multiple ways that together form the boundary - each way is a piece of the boundary line, not a separate polygon. I need to merge all the outer ways into a single ring, and inner ways into holes. Actually, looking at the data more carefully - for Suginami (relation 1543055), the overpass query returned members where each member is a way with geometry. But some of these may be individual ways that together form the boundary. The issue is that the overpass `geom` output gives us the geometry of each member way separately. For a boundary relation, the outer members are line segments that together form the closed boundary ring(s). I need to merge them. Let me fix the parsing. For relations that describe ward boundaries (type=boundary, boundary=administrative), the members with role "outer" are ways that together form the outer boundary. I need to merge them into one or more closed rings. Actually, I think the issue might be that some of these members are individual polygons (like small exclaves or islands). The relation for Suginami has 166 members with "outer" role that are each separate polygons? That seems odd. Let me check what the relation looks like more carefully.
24:28
Bash
python3 -c "
python3 -c "
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Let's see the actual members of Suginami relation
query = '[out:json];rel(id:1543055);(._;>>;);out body;'
resp = requests.post(url, data={'data': query}, timeout=60, headers=headers)
data = resp.json()
for el in data.get('elements', []):
if el['type'] == 'relation':
print(f'Relation {el[\"id\"]}: tags={el.get(\"tags\",{}).get(\"name\",\"\")}')
print(f' Members ({len(el.get(\"members\",[]))}):')
for m in el.get('members', [])[:10]:
print(f' type={m[\"type\"]} ref={m[\"ref\"]} role={m[\"role\"]}')
if len(el.get('members',[])) > 10:
print(f' ... and {len(el.get(\"members\",[]))-10} more')
" 2>&1
Relation 1543055: tags=杉並区
Members (120):
type=node ref=718383251 role=label
type=way ref=1356066346 role=outer
type=way ref=1259558470 role=outer
type=way ref=1259558476 role=outer
type=way ref=1259565094 role=outer
type=way ref=1356066343 role=outer
type=way ref=1259565085 role=outer
type=way ref=1356066341 role=outer
type=way ref=1263119577 role=outer
type=way ref=1356066361 role=outer
... and 110 more
Relation 17146535: tags=方南一丁目
Members (9):
type=node ref=8370281480 role=label
type=way ref=1245966646 role=outer
type=way ref=1353141929 role=outer
type=way ref=1245966644 role=outer
type=way ref=355213612 role=outer
type=way ref=1241585299 role=outer
type=way ref=1245966643 role=outer
type=way ref=1245966647 role=outer
type=way ref=1254497155 role=outer
Relation 17146536: tags=方南二丁目
Members (7):
type=node ref=8370281479 role=label
type=way ref=1245966646 role=outer
type=way ref=1245966642 role=outer
type=way ref=1254497146 role=outer
type=way ref=1245966648 role=outer
type=way ref=1353141923 role=outer
type=way ref=1245966645 role=outer
Relation 17146537: tags=方南
Members (15):
type=way ref=1245966647 role=outer
type=way ref=1254497155 role=outer
type=way ref=1245966642 role=outer
type=way ref=1254497146 role=outer
type=way ref=1245966648 role=outer
type=way ref=1353141923 role=outer
type=way ref=1245966645 role=outer
type=way ref=1353141929 role=outer
type=way ref=1245966644 role=outer
type=way ref=355213612 role=outer
... and 5 more
Relation 17264720: tags=和田一丁目
Members (7):
type=way ref=1254497136 role=outer
type=way ref=1354063715 role=outer
type=way ref=1353141937 role=outer
type=way ref=1353141939 role=outer
type=way ref=1254497141 role=outer
type=way ref=1254497140 role=outer
type=node ref=8370281461 role=label
Relation 17264721: tags=和田二丁目
Members (7):
type=node ref=8370281460 role=label
type=way ref=1254497141 role=outer
type=way ref=1254497142 role=outer
type=way ref=1245966648 role=outer
type=way ref=1254497137 role=outer
type=way ref=1254497175 role=outer
type=way ref=1254497138 role=outer
Relation 17264722: tags=和田三丁目
Members (7):
type=node ref=8370281459 role=label
type=way ref=1254497170 role=outer
type=way ref=1254497135 role=outer
type=way ref=1254497143 role=outer
type=way ref=1254497139 role=outer
type=way ref=1254497140 role=outer
type=way ref=1254497138 role=outer
Relation 17264723: tags=和田
Members (16):
type=node ref=1629628690 role=label
type=way ref=1254497137 role=outer
type=way ref=1254497175 role=outer
type=way ref=1254497170 role=outer
type=way ref=1254497135 role=outer
type=way ref=1254497143 role=outer
type=way ref=1254497139 role=outer
type=way ref=1254497136 role=outer
type=way ref=1354063715 role=outer
type=way ref=1353141937 role=outer
... and 6 more
Relation 17264724: tags=和泉一丁目
Members (6):
type=way ref=1334061143 role=outer
type=way ref=1254497153 role=outer
type=way ref=1245966647 role=outer
type=way ref=1254497152 role=outer
type=way ref=1254497150 role=outer
type=node ref=8389193268 role=label
Relation 17264725: tags=和泉二丁目
Members (7):
type=node ref=8389193267 role=label
type=way ref=1254497150 role=outer
type=way ref=1268111876 role=outer
type=way ref=1254497145 role=outer
type=way ref=1254497148 role=outer
type=way ref=1254497151 role=outer
type=way ref=1254497154 role=outer
Relation 17264726: tags=和泉三丁目
Members (8):
type=node ref=8389193266 role=label
type=way ref=1254497151 role=outer
type=way ref=1254497161 role=outer
type=way ref=1254497167 role=outer
type=way ref=1254497144 role=outer
type=way ref=1254497169 role=outer
type=way ref=1254497158 role=outer
type=way ref=1254497147 role=outer
Relation 17264727: tags=和泉四丁目
Members (8):
type=way ref=1254497147 role=outer
type=way ref=1254497149 role=outer
type=way ref=1245966642 role=outer
type=way ref=1254497155 role=outer
type=way ref=1254497153 role=outer
type=way ref=1334061143 role=outer
type=way ref=1254497154 role=outer
type=node ref=8389193265 role=label
Relation 17264728: tags=和泉
Members (18):
type=node ref=11708116687 role=label
type=way ref=1254497148 role=outer
type=way ref=1254497161 role=outer
type=way ref=1254497167 role=outer
type=way ref=1254497144 role=outer
type=way ref=1254497169 role=outer
type=way ref=1254497158 role=outer
type=way ref=1254497149 role=outer
type=way ref=1245966642 role=outer
type=way ref=1254497155 role=outer
... and 8 more
Relation 17264729: tags=永福一丁目
Members (7):
type=way ref=1254497157 role=outer
type=way ref=1330134189 role=outer
type=way ref=1254497159 role=outer
type=way ref=1254497160 role=outer
type=way ref=1254497161 role=outer
type=way ref=1254497148 role=outer
type=node ref=8389193264 role=label
Relation 17264730: tags=永福二丁目
Members (6):
type=node ref=8389193263 role=label
type=way ref=1254497167 role=outer
type=way ref=1254497160 role=outer
type=way ref=1254497166 role=outer
type=way ref=1254497165 role=outer
type=way ref=1254497164 role=outer
Relation 17264731: tags=永福三丁目
Members (8):
type=node ref=8389193262 role=label
type=way ref=1254511642 role=outer
type=way ref=1254497162 role=outer
type=way ref=1254497165 role=outer
type=way ref=1256772969 role=outer
type=way ref=1256743843 role=outer
type=way ref=1254497156 role=outer
type=way ref=1256743846 role=outer
Relation 17264732: tags=永福四丁目
Members (5):
type=way ref=1254497163 role=outer
type=way ref=1254497144 role=outer
type=way ref=1254497164 role=outer
type=way ref=1254497162 role=outer
type=node ref=8389193261 role=label
Relation 17264733: tags=永福
Members (19):
type=node ref=1623005135 role=label
type=way ref=1254511642 role=outer
type=way ref=1254497163 role=outer
type=way ref=1254497144 role=outer
type=way ref=1254497167 role=outer
type=way ref=1254497161 role=outer
type=way ref=1254497148 role=outer
type=way ref=1254497157 role=outer
type=way ref=1330134189 role=outer
type=way ref=1254497159 role=outer
... and 9 more
Relation 17264734: tags=堀ノ内一丁目
Members (6):
type=way ref=1254497146 role=outer
type=way ref=1254497149 role=outer
type=way ref=1254497158 role=outer
type=way ref=1254497172 role=outer
type=way ref=1254497171 role=outer
type=node ref=8389193241 role=label
Relation 17264735: tags=堀ノ内二丁目
Members (6):
type=node ref=8389193240 role=label
type=way ref=1254497171 role=outer
type=way ref=1254511641 role=outer
type=way ref=1254497174 role=outer
type=way ref=1254497173 role=outer
type=way ref=1254497137 role=outer
Relation 17264736: tags=堀ノ内三丁目
Members (7):
type=node ref=8389193239 role=label
type=way ref=1254497173 role=outer
type=way ref=1254540997 role=outer
type=way ref=1254511636 role=outer
type=way ref=1254497168 role=outer
type=way ref=1254497170 role=outer
type=way ref=1254497175 role=outer
Relation 17264737: tags=堀ノ内
Members (16):
type=way ref=1254497172 role=outer
type=way ref=1254511641 role=outer
type=way ref=1254497174 role=outer
type=way ref=1254540997 role=outer
type=way ref=1254511636 role=outer
type=way ref=1254497168 role=outer
type=way ref=1254497170 role=outer
type=way ref=1254497175 role=outer
type=way ref=1254497137 role=outer
type=way ref=1254497146 role=outer
... and 6 more
Relation 17264835: tags=梅里一丁目
Members (6):
type=way ref=1254511638 role=outer
type=way ref=1254511639 role=outer
type=way ref=1254497135 role=outer
type=way ref=1254497168 role=outer
type=way ref=1254511637 role=outer
type=node ref=8389193243 role=label
Relation 17264836: tags=梅里二丁目
Members (8):
type=node ref=8389193242 role=label
type=way ref=1256743826 role=outer
type=way ref=1256743818 role=outer
type=way ref=1254511635 role=outer
type=way ref=1259558456 role=outer
type=way ref=1259558462 role=outer
type=way ref=1254511638 role=outer
type=way ref=1254540991 role=outer
Relation 17264837: tags=梅里
Members (13):
type=node ref=1623002269 role=label
type=way ref=1256743826 role=outer
type=way ref=1256743818 role=outer
type=way ref=1254511635 role=outer
type=way ref=1259558456 role=outer
type=way ref=1259558462 role=outer
type=way ref=1254511639 role=outer
type=way ref=1254497135 role=outer
type=way ref=1254497168 role=outer
type=way ref=1254511637 role=outer
... and 3 more
Relation 17264838: tags=大宮一丁目
Members (7):
type=node ref=8389193235 role=label
type=way ref=1254511643 role=outer
type=way ref=1254511644 role=outer
type=way ref=1254540994 role=outer
type=way ref=1254511641 role=outer
type=way ref=1254497172 role=outer
type=way ref=1254497169 role=outer
Relation 17264839: tags=大宮二丁目
Members (10):
type=node ref=8389193234 role=label
type=way ref=1254497163 role=outer
type=way ref=1254511643 role=outer
type=way ref=1254540990 role=outer
type=way ref=1256743819 role=outer
type=way ref=1256743822 role=outer
type=way ref=1254511640 role=outer
type=way ref=1256743831 role=outer
type=way ref=1256743851 role=outer
type=way ref=1254511642 role=outer
Relation 17264840: tags=大宮
Members (16):
type=node ref=1623008351 role=label
type=way ref=1256743851 role=outer
type=way ref=1256743831 role=outer
type=way ref=1254511640 role=outer
type=way ref=1256743822 role=outer
type=way ref=1256743819 role=outer
type=way ref=1254540990 role=outer
type=way ref=1254511644 role=outer
type=way ref=1254540994 role=outer
type=way ref=1254511641 role=outer
... and 6 more
Relation 17265255: tags=松ノ木一丁目
Members (5):
type=way ref=1254540992 role=outer
type=way ref=1254511644 role=outer
type=way ref=1254540990 role=outer
type=way ref=1254540993 role=outer
type=node ref=8389193238 role=label
Relation 17265256: tags=松ノ木二丁目
Members (7):
type=way ref=1254540992 role=outer
type=way ref=1254540989 role=outer
type=way ref=1254540996 role=outer
type=way ref=1254540997 role=outer
type=way ref=1254497174 role=outer
type=way ref=1254540994 role=outer
type=node ref=8389193237 role=label
Relation 17265257: tags=松ノ木三丁目
Members (6):
type=way ref=1254540996 role=outer
type=way ref=1254540995 role=outer
type=way ref=1254540991 role=outer
type=way ref=1254511637 role=outer
type=way ref=1254511636 role=outer
type=node ref=8389193236 role=label
Relation 17265258: tags=松ノ木
Members (15):
type=way ref=1254540993 role=outer
type=way ref=1254540989 role=outer
type=way ref=1254540995 role=outer
type=way ref=1254540991 role=outer
type=way ref=1254511637 role=outer
type=way ref=1254511636 role=outer
type=way ref=1254540997 role=outer
type=way ref=1254497174 role=outer
type=way ref=1254540994 role=outer
type=way ref=1254511644 role=outer
... and 5 more
Relation 17294912: tags=成田東一丁目
Members (7):
type=way ref=1256743819 role=outer
type=way ref=1256743823 role=outer
type=way ref=1256743820 role=outer
type=way ref=1254540995 role=outer
type=way ref=1254540989 role=outer
type=way ref=1254540993 role=outer
type=node ref=8389193248 role=label
Relation 17294913: tags=成田東二丁目
Members (5):
type=way ref=1256743823 role=outer
type=way ref=1256743822 role=outer
type=way ref=1256743821 role=outer
type=way ref=1256743824 role=outer
type=node ref=8389193247 role=label
Relation 17294914: tags=成田東三丁目
Members (6):
type=way ref=1256743824 role=outer
type=way ref=1256743830 role=outer
type=way ref=1256743825 role=outer
type=way ref=1256743826 role=outer
type=way ref=1256743820 role=outer
type=node ref=8389193246 role=label
Relation 17294915: tags=成田東四丁目
Members (7):
type=node ref=8389193245 role=label
type=way ref=1256743825 role=outer
type=way ref=1256743839 role=outer
type=way ref=1256743827 role=outer
type=way ref=1256743828 role=outer
type=way ref=1256743829 role=outer
type=way ref=1256743818 role=outer
Relation 17294916: tags=成田東五丁目
Members (6):
type=node ref=8389193244 role=label
type=way ref=1256743828 role=outer
type=way ref=1259565079 role=outer
type=way ref=1259212444 role=outer
type=way ref=1256743817 role=outer
type=way ref=1256743833 role=outer
Relation 17294917: tags=成田東
Members (22):
type=way ref=1256743821 role=outer
type=way ref=1256743830 role=outer
type=way ref=1256743839 role=outer
type=way ref=1256743827 role=outer
type=way ref=1256743833 role=outer
type=way ref=1256743817 role=outer
type=way ref=1259212444 role=outer
type=way ref=1259565079 role=outer
type=way ref=1256743829 role=outer
type=way ref=1256743818 role=outer
... and 12 more
Relation 17294918: tags=成田西一丁目
Members (6):
type=node ref=8389193252 role=label
type=way ref=1256743821 role=outer
type=way ref=1254511640 role=outer
type=way ref=1256743838 role=outer
type=way ref=1256743836 role=outer
type=way ref=1256743837 role=outer
Relation 17294919: tags=成田西二丁目
Members (5):
type=node ref=8389193251 role=label
type=way ref=1256743836 role=outer
type=way ref=1256743832 role=outer
type=way ref=1256743844 role=outer
type=way ref=1256743834 role=outer
Relation 17294920: tags=成田西三丁目
Members (8):
type=node ref=8389193250 role=label
type=way ref=1256743834 role=outer
type=way ref=1257073559 role=outer
type=way ref=1256743841 role=outer
type=way ref=1256743840 role=outer
type=way ref=1256743839 role=outer
type=way ref=1256743830 role=outer
type=way ref=1256743837 role=outer
Relation 17294921: tags=成田西四丁目
Members (6):
type=node ref=8389193249 role=label
type=way ref=1256743840 role=outer
type=way ref=1256743835 role=outer
type=way ref=1259212447 role=outer
type=way ref=1256743833 role=outer
type=way ref=1256743827 role=outer
Relation 17294922: tags=成田西
Members (18):
type=node ref=11710500717 role=label
type=way ref=1256743838 role=outer
type=way ref=1256743832 role=outer
type=way ref=1256743844 role=outer
type=way ref=1257073559 role=outer
type=way ref=1256743841 role=outer
type=way ref=1256743835 role=outer
type=way ref=1259212447 role=outer
type=way ref=1256743833 role=outer
type=way ref=1256743827 role=outer
... and 8 more
Relation 17294923: tags=浜田山一丁目
Members (6):
type=node ref=8389193272 role=label
type=way ref
[... truncated ...]ation 17326983: tags=阿佐谷南
Members (17):
type=node ref=1622922196 role=label
type=way ref=1259565078 role=outer
type=way ref=1259565075 role=outer
type=way ref=1259565084 role=outer
type=way ref=1259592375 role=outer
type=way ref=1259212450 role=outer
type=way ref=1259212453 role=outer
type=way ref=1259212444 role=outer
type=way ref=1259565079 role=outer
type=way ref=1256743829 role=outer
... and 7 more
Relation 17326984: tags=阿佐谷北一丁目
Members (7):
type=way ref=1259565078 role=outer
type=way ref=1259565088 role=outer
type=way ref=1259565098 role=outer
type=way ref=1259565097 role=outer
type=way ref=1259565092 role=outer
type=way ref=1259565091 role=outer
type=node ref=8389193302 role=label
Relation 17326985: tags=阿佐谷北二丁目
Members (5):
type=way ref=1259565088 role=outer
type=way ref=1259565075 role=outer
type=way ref=1259565090 role=outer
type=way ref=1259565089 role=outer
type=node ref=8389193301 role=label
Relation 17326986: tags=阿佐谷北三丁目
Members (6):
type=node ref=8389193300 role=label
type=way ref=1259565089 role=outer
type=way ref=1259565083 role=outer
type=way ref=1259592372 role=outer
type=way ref=1259565086 role=outer
type=way ref=1259565098 role=outer
Relation 17326987: tags=阿佐谷北四丁目
Members (6):
type=way ref=1259565097 role=outer
type=way ref=1259565086 role=outer
type=way ref=1259565096 role=outer
type=way ref=1259565095 role=outer
type=way ref=1259565099 role=outer
type=node ref=8389193299 role=label
Relation 17326988: tags=阿佐谷北五丁目
Members (6):
type=way ref=1259565099 role=outer
type=way ref=1259565093 role=outer
type=way ref=1259565094 role=outer
type=way ref=1259558468 role=outer
type=way ref=1259565092 role=outer
type=node ref=8389193298 role=label
Relation 17326989: tags=阿佐谷北六丁目
Members (7):
type=way ref=1259565095 role=outer
type=way ref=1259565087 role=outer
type=way ref=1356066341 role=outer
type=way ref=1259565085 role=outer
type=way ref=1356066343 role=outer
type=way ref=1259565093 role=outer
type=node ref=8389193297 role=label
Relation 17326990: tags=阿佐谷北
Members (20):
type=node ref=1622911239 role=label
type=way ref=1259565090 role=outer
type=way ref=1259565083 role=outer
type=way ref=1259592372 role=outer
type=way ref=1259565096 role=outer
type=way ref=1259565087 role=outer
type=way ref=1356066341 role=outer
type=way ref=1259565085 role=outer
type=way ref=1356066343 role=outer
type=way ref=1259565094 role=outer
... and 10 more
Relation 17327874: tags=天沼一丁目
Members (6):
type=way ref=1259565084 role=outer
type=way ref=1259592374 role=outer
type=way ref=1259592376 role=outer
type=way ref=1259565083 role=outer
type=way ref=1259565090 role=outer
type=node ref=8389193308 role=label
Relation 17327875: tags=天沼二丁目
Members (7):
type=node ref=8389193307 role=label
type=way ref=1259592374 role=outer
type=way ref=1259592375 role=outer
type=way ref=1259592378 role=outer
type=way ref=1259592377 role=outer
type=way ref=1259592386 role=outer
type=way ref=1259592379 role=outer
Relation 17327876: tags=天沼三丁目
Members (5):
type=node ref=8389193306 role=label
type=way ref=1259592377 role=outer
type=way ref=1259344146 role=outer
type=way ref=1259592371 role=outer
type=way ref=1259592381 role=outer
Relation 17327877: tags=天沼
Members (15):
type=node ref=1622932562 role=label
type=way ref=1259592378 role=outer
type=way ref=1259344146 role=outer
type=way ref=1259592371 role=outer
type=way ref=1259592381 role=outer
type=way ref=1259592386 role=outer
type=way ref=1259592379 role=outer
type=way ref=1259592376 role=outer
type=way ref=1259565083 role=outer
type=way ref=1259565090 role=outer
... and 5 more
Relation 17327878: tags=本天沼一丁目
Members (5):
type=way ref=1259592382 role=outer
type=way ref=1259592383 role=outer
type=way ref=1259565096 role=outer
type=way ref=1259592372 role=outer
type=node ref=8389193311 role=label
Relation 17327879: tags=本天沼二丁目
Members (7):
type=node ref=8389193310 role=label
type=way ref=1259592384 role=outer
type=way ref=1263119578 role=outer
type=way ref=1259592385 role=outer
type=way ref=1259592382 role=outer
type=way ref=1259592376 role=outer
type=way ref=1259592379 role=outer
Relation 17327880: tags=本天沼三丁目
Members (9):
type=node ref=8389193309 role=label
type=way ref=1259592384 role=outer
type=way ref=1259592388 role=outer
type=way ref=1263119584 role=outer
type=way ref=1259592394 role=outer
type=way ref=1259592386 role=outer
type=way ref=1259592381 role=outer
type=way ref=1259592390 role=outer
type=way ref=1259592380 role=outer
Relation 17327881: tags=本天沼
Members (18):
type=node ref=11710482764 role=label
type=way ref=1259592390 role=outer
type=way ref=1259592380 role=outer
type=way ref=1259592394 role=outer
type=way ref=1263119584 role=outer
type=way ref=1259592388 role=outer
type=way ref=1263119578 role=outer
type=way ref=1259592385 role=outer
type=way ref=1259592383 role=outer
type=way ref=1259565096 role=outer
... and 8 more
Relation 17327882: tags=清水一丁目
Members (6):
type=node ref=8389193305 role=label
type=way ref=1259592391 role=outer
type=way ref=1259592389 role=outer
type=way ref=1259592390 role=outer
type=way ref=1259592371 role=outer
type=way ref=1259592373 role=outer
Relation 17327883: tags=清水二丁目
Members (5):
type=way ref=1259592389 role=outer
type=way ref=1259592393 role=outer
type=way ref=1259592392 role=outer
type=way ref=1259592380 role=outer
type=node ref=8389193304 role=label
Relation 17327884: tags=清水三丁目
Members (6):
type=node ref=8389193303 role=label
type=way ref=1259592392 role=outer
type=way ref=1259592395 role=outer
type=way ref=1263119586 role=outer
type=way ref=1259592387 role=outer
type=way ref=1259592394 role=outer
Relation 17327885: tags=清水
Members (14):
type=node ref=1629604003 role=label
type=way ref=1259592393 role=outer
type=way ref=1259592391 role=outer
type=way ref=1259592373 role=outer
type=way ref=1259592371 role=outer
type=way ref=1259592390 role=outer
type=way ref=1259592380 role=outer
type=way ref=1259592394 role=outer
type=way ref=1259592387 role=outer
type=way ref=1263119586 role=outer
... and 4 more
Relation 17349405: tags=桃井一丁目
Members (6):
type=node ref=8389196236 role=label
type=way ref=1259344149 role=outer
type=way ref=1262948262 role=outer
type=way ref=1262948263 role=outer
type=way ref=1259592393 role=outer
type=way ref=1259592391 role=outer
Relation 17349406: tags=桃井二丁目
Members (5):
type=node ref=8389196235 role=label
type=way ref=1262948262 role=outer
type=way ref=1259344151 role=outer
type=way ref=1262948264 role=outer
type=way ref=1262948265 role=outer
Relation 17349407: tags=桃井三丁目
Members (5):
type=way ref=1262948264 role=outer
type=way ref=1259267216 role=outer
type=way ref=1262948266 role=outer
type=way ref=1262948267 role=outer
type=node ref=8389196234 role=label
Relation 17349408: tags=桃井四丁目
Members (6):
type=node ref=8389196233 role=label
type=way ref=1262948266 role=outer
type=way ref=1262948268 role=outer
type=way ref=1259267221 role=outer
type=way ref=1262948260 role=outer
type=way ref=1263101090 role=outer
Relation 17349409: tags=桃井
Members (17):
type=node ref=11710458175 role=label
type=way ref=1259344149 role=outer
type=way ref=1259592391 role=outer
type=way ref=1259592393 role=outer
type=way ref=1262948263 role=outer
type=way ref=1262948265 role=outer
type=way ref=1262948267 role=outer
type=way ref=1263101090 role=outer
type=way ref=1262948260 role=outer
type=way ref=1259267221 role=outer
... and 7 more
Relation 17351014: tags=今川一丁目
Members (5):
type=way ref=1262948263 role=outer
type=way ref=1263101084 role=outer
type=way ref=1263101085 role=outer
type=way ref=1259592395 role=outer
type=node ref=8389196240 role=label
Relation 17351015: tags=今川二丁目
Members (5):
type=way ref=1263101084 role=outer
type=way ref=1262948265 role=outer
type=way ref=1263101086 role=outer
type=way ref=1263101087 role=outer
type=node ref=8389196239 role=label
Relation 17351016: tags=今川三丁目
Members (6):
type=way ref=1263101086 role=outer
type=way ref=1262948267 role=outer
type=way ref=1263101090 role=outer
type=way ref=1263101088 role=outer
type=way ref=1263101089 role=outer
type=node ref=8389196238 role=label
Relation 17351017: tags=今川四丁目
Members (6):
type=node ref=8389196237 role=label
type=way ref=1263101088 role=outer
type=way ref=1262948260 role=outer
type=way ref=1262948261 role=outer
type=way ref=1263101083 role=outer
type=way ref=1263119593 role=outer
Relation 17351018: tags=今川
Members (17):
type=node ref=1623000182 role=label
type=way ref=1263101085 role=outer
type=way ref=1263101087 role=outer
type=way ref=1263101089 role=outer
type=way ref=1263119593 role=outer
type=way ref=1263101083 role=outer
type=way ref=1262948261 role=outer
type=way ref=1262948260 role=outer
type=way ref=1263101090 role=outer
type=way ref=1262948267 role=outer
... and 7 more
Relation 17351203: tags=下井草一丁目
Members (7):
type=way ref=1259592385 role=outer
type=way ref=1263119579 role=outer
type=way ref=1263119576 role=outer
type=way ref=1263119577 role=outer
type=way ref=1259565087 role=outer
type=way ref=1259592383 role=outer
type=node ref=8389196253 role=label
Relation 17351204: tags=下井草二丁目
Members (6):
type=way ref=355213613 role=outer
type=way ref=1356066361 role=outer
type=way ref=1263119576 role=outer
type=way ref=1263119580 role=outer
type=way ref=1263119581 role=outer
type=node ref=8389196252 role=label
Relation 17351205: tags=下井草三丁目
Members (7):
type=way ref=1263119580 role=outer
type=way ref=1263119579 role=outer
type=way ref=1263119578 role=outer
type=way ref=1259592388 role=outer
type=way ref=1263119582 role=outer
type=way ref=1263119583 role=outer
type=node ref=8389196251 role=label
Relation 17351206: tags=下井草四丁目
Members (6):
type=way ref=1263119582 role=outer
type=way ref=1263119584 role=outer
type=way ref=1259592387 role=outer
type=way ref=1263119585 role=outer
type=way ref=1263119571 role=outer
type=node ref=8389196250 role=label
Relation 17351207: tags=下井草五丁目
Members (5):
type=way ref=1263119585 role=outer
type=way ref=1263119586 role=outer
type=way ref=1263119574 role=outer
type=way ref=1263119587 role=outer
type=node ref=8389196249 role=label
Relation 17351208: tags=下井草
Members (22):
type=node ref=1629609990 role=label
type=way ref=1263119574 role=outer
type=way ref=1263119586 role=outer
type=way ref=1259592387 role=outer
type=way ref=1263119584 role=outer
type=way ref=1259592388 role=outer
type=way ref=1263119578 role=outer
type=way ref=1259592385 role=outer
type=way ref=1259592383 role=outer
type=way ref=1259565087 role=outer
... and 12 more
Relation 17351209: tags=上井草一丁目
Members (6):
type=node ref=8389196248 role=label
type=way ref=1263119574 role=outer
type=way ref=1263101085 role=outer
type=way ref=1263119588 role=outer
type=way ref=1263119589 role=outer
type=way ref=1263119603 role=outer
Relation 17351210: tags=上井草二丁目
Members (5):
type=way ref=1263119588 role=outer
type=way ref=1263101087 role=outer
type=way ref=1263119590 role=outer
type=way ref=1263119575 role=outer
type=node ref=8389196247 role=label
Relation 17351211: tags=上井草三丁目
Members (8):
type=way ref=1263119590 role=outer
type=way ref=1263101089 role=outer
type=way ref=1263119593 role=outer
type=way ref=1263119592 role=outer
type=way ref=1350815561 role=outer
type=way ref=1263119594 role=outer
type=way ref=1263119591 role=outer
type=node ref=8389196246 role=label
Relation 17351212: tags=上井草四丁目
Members (7):
type=way ref=1263119592 role=outer
type=way ref=1263101083 role=outer
type=way ref=1259344159 role=outer
type=way ref=1259344157 role=outer
type=way ref=1350891129 role=outer
type=way ref=1263119573 role=outer
type=node ref=8389196245 role=label
Relation 17351213: tags=上井草
Members (21):
type=node ref=11710523592 role=label
type=way ref=1263119574 role=outer
type=way ref=1263119603 role=outer
type=way ref=1263119589 role=outer
type=way ref=1263119575 role=outer
type=way ref=1263119591 role=outer
type=way ref=1263119594 role=outer
type=way ref=1350815561 role=outer
type=way ref=1350891129 role=outer
type=way ref=1263119573 role=outer
... and 11 more
Relation 17351214: tags=井草一丁目
Members (7):
type=way ref=1356080194 role=outer
type=way ref=1263119572 role=outer
type=way ref=1263119581 role=outer
type=way ref=1263119583 role=outer
type=way ref=1263119595 role=outer
type=way ref=1263119596 role=outer
type=node ref=8389196258 role=label
Relation 17351215: tags=井草二丁目
Members (5):
type=way ref=1263119595 role=outer
type=way ref=1263119571 role=outer
type=way ref=1263119597 role=outer
type=way ref=1263119598 role=outer
type=node ref=8389196257 role=label
Relation 17351216: tags=井草三丁目
Members (7):
type=way ref=1263119597 role=outer
type=way ref=1263119587 role=outer
type=way ref=1263119603 role=outer
type=way ref=1263119599 role=outer
type=way ref=1350815547 role=outer
type=way ref=1263119600 role=outer
type=node ref=8389196256 role=label
Relation 17351217: tags=井草四丁目
Members (6):
type=way ref=1263119599 role=outer
type=way ref=1263119589 role=outer
type=way ref=1263119601 role=outer
type=way ref=1350815575 role=outer
type=way ref=1263119602 role=outer
type=node ref=8389196255 role=label
Relation 17351218: tags=井草五丁目
Members (5):
type=way ref=1263119601 role=outer
type=way ref=1263119575 role=outer
type=way ref=1263119591 role=outer
type=way ref=355213605 role=outer
type=node ref=8389196254 role=label
Relation 17351219: tags=井草
Members (23):
type=node ref=1622948061 role=label
type=way ref=1263119572 role=outer
type=way ref=1356080194 role=outer
type=way ref=1263119596 role=outer
type=way ref=1263119598 role=outer
type=way ref=1263119600 role=outer
type=way ref=1350815547 role=outer
type=way ref=1263119602 role=outer
type=way ref=1350815575 role=outer
type=way ref=355213605 role=outer
... and 13 more
24:52
I see the problem now. The Overpass query with `out body geom` for the relation is returning child relations (like sub-ward neighborhoods) as well because of the recursive nature. The Suginami relation (1543055) has 120 members, and many of those members are ways (outer/inner) but the Overpass response is also including sub-relations that are administrative boundary relations for neighborhoods within Suginami. This is because the query fetches the relation itself but then `out body geom` on the relation only outputs the relation's geometry (the ways), not child relations. Wait, looking at the output more carefully, the response had 21 elements earlier when I queried all 23 relations together. That means 2 were missing (Taito and Koto). Actually, the missing ones (1759476, 1759478) had different IDs - 1758888 (Taito) and 3554015 (Koto). Let me check again... Actually the issue is different. When I queried `rel(id:1759476,1759478);out body geom;` it returned nothing because those IDs don't exist. The correct IDs are 1758888 for Taito and 3554015 for Koto. But looking at the geometry issue - the Suginami ward area was very small (0.86 km² when it should be ~34 km²). The Overpass query with `out body geom` on relations returns the geometry of the ways that make up the relation. But looking at the output, the relation has 120 members, most of which are "outer" ways. The geometry from `out body geom` should concatenate all those way geometries. Wait, looking more carefully at the Overpass output for Suginami relation, it seems like the issue is that the relation query is also returning sub-relations (neighborhood boundaries) that are members of the ward relation, and these sub-relations contain their own outer rings that are being parsed incorrectly. Actually no - when using `out body geom` on a relation in Overpass, it returns the relation's tags and its member geometries. Each member way's full node coordinates are included. The outer members should together form the ward boundary. But looking at the output, the relation has 120 members and many of them are ways serving as outer boundaries. The issue might be that the relation includes neighborhood sub-areas (chome/丁目 boundaries) as separate polygons, not just the outer boundary of the ward. Wait, that makes sense! The Suginami ward relation (1543055) has 120 members, and these include: - Outer ways that form the ward boundary - Inner ways that represent holes (exclaves) - Sub-area relations (neighborhoods) that are also part of the relation The sub-area relations for neighborhoods (like 方南一丁目, 方南二丁目, etc.) are also members of the ward relation. These are additional child relations that define sub-areas within the ward. These sub-area relations have their own outer ways which represent the neighborhood boundaries within Suginami. So when I parse the geometry, I'm treating ALL outer members as the ward boundary, when in reality the ward boundary is defined only by the main outer ways, and the sub-relation members are separate internal boundaries. The solution is to filter the members: only process members that are ways (type=way), not relations. The sub-areas are included as relation members, not way members. Let me check this hypothesis.
24:52
24:52
Bash
python3 -c "
python3 -c "
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Check Suginami relation - what are the member types?
query = '[out:json];rel(id:1543055);out body;'
resp = requests.post(url, data={'data': query}, timeout=60, headers=headers)
data = resp.json()
for el in data.get('elements', []):
if el['type'] == 'relation' and el['id'] == 1543055:
members = el.get('members', [])
way_outer = [m for m in members if m['type'] == 'way' and m['role'] == 'outer']
way_inner = [m for m in members if m['type'] == 'way' and m['role'] == 'inner']
rel_outer = [m for m in members if m['type'] == 'relation' and m['role'] == 'outer']
rel_inner = [m for m in members if m['type'] == 'relation' and m['role'] == 'inner']
print(f'Suginami relation:')
print(f' Total members: {len(members)}')
print(f' Way outer: {len(way_outer)}')
print(f' Way inner: {len(way_inner)}')
print(f' Relation outer: {len(rel_outer)}')
print(f' Relation inner: {len(rel_inner)}')
for m in members:
if m['type'] != 'way':
print(f' member type={m[\"type\"]} ref={m[\"ref\"]} role={m[\"role\"]}')
" 2>&1
Suginami relation:
Total members: 120
Way outer: 83
Way inner: 0
Relation outer: 0
Relation inner: 0
member type=node ref=718383251 role=label
member type=relation ref=17146537 role=subarea
member type=relation ref=17264723 role=subarea
member type=relation ref=17264728 role=subarea
member type=relation ref=17264733 role=subarea
member type=relation ref=17264737 role=subarea
member type=relation ref=17264837 role=subarea
member type=relation ref=17264840 role=subarea
member type=relation ref=17265258 role=subarea
member type=relation ref=17294917 role=subarea
member type=relation ref=17294922 role=subarea
member type=relation ref=17294927 role=subarea
member type=relation ref=17295267 role=subarea
member type=relation ref=17299394 role=subarea
member type=relation ref=17299399 role=subarea
member type=relation ref=17299403 role=subarea
member type=relation ref=17299409 role=subarea
member type=relation ref=17299510 role=subarea
member type=relation ref=17299514 role=subarea
member type=relation ref=17323138 role=subarea
member type=relation ref=17323143 role=subarea
member type=relation ref=17323281 role=subarea
member type=relation ref=17324147 role=subarea
member type=relation ref=17324152 role=subarea
member type=relation ref=17324157 role=subarea
member type=relation ref=17326884 role=subarea
member type=relation ref=17326889 role=subarea
member type=relation ref=17326983 role=subarea
member type=relation ref=17326990 role=subarea
member type=relation ref=17327877 role=subarea
member type=relation ref=17327881 role=subarea
member type=relation ref=17327885 role=subarea
member type=relation ref=17349409 role=subarea
member type=relation ref=17351018 role=subarea
member type=relation ref=17351208 role=subarea
member type=relation ref=17351213 role=subarea
member type=relation ref=17351219 role=subarea
25:42
I see the issue. The Suginami ward relation has: - 83 way members with role=outer (these are the actual boundary line segments) - 36 relation members with role=subarea (these are sub-neighborhood relations - I should NOT follow these) - 1 node member with role=label My current code was treating all `role="outer"` members, but the subarea relations also have their own `role="outer"` way members, which is why we got 166 tiny polygons instead of one large polygon. I need to only process way members (type=way) with role=outer/inner, and skip relation members and node members. Let me fix this.
25:42
Edit
/work/solve.py
Edited /work/solve.py.
25:56
Now I need to also fix the query - when we used `out body geom`, the sub-relations don't have geometry for the relation members. But since we're only processing way members now, this should be fine. However, some wards might have their boundary as relations (sub-area relations). Let me check: Actually, looking at the Suginami example, the boundary ways are `way` members with `role=outer`. The `relation` members have `role=subarea` and they represent neighborhoods within Suginami. These sub-relations are NOT part of the outer boundary - they're separate entities. So fixing the code to only look at `type=way` members should work. The 83 way outer members should define the full boundary correctly. Let me re-run the script.
25:56
Bash
cd /work && python solve.py 2>&1
============================================================
Tokyo 23 Special Wards — Road Crossing Analysis
============================================================
[1/5] Fetching ward boundaries from OSM…
✓ Layer 'wards_jgd' saved (23 wards)
[2/5] Fetching drivable road network from OSM…
Bounding box: 35.4817, 35.8175, 139.5629, 139.9189
Sending Overpass query for drivable roads…
Received 127738 road elements
Road segments loaded: 127738
Roads after clipping: 4338
Simplified to 4634 single-part line segments
[3/5] Detecting road / ward-boundary crossings…
Found 5462 crossing points
✓ Layer 'crossing_points' saved (5462 points)
[4/5] Creating 50 m buffers and clipping to wards…
✓ Layer 'crossing_buffers_50m' saved (5462 buffers)
✓ Layer 'buffer_ward_intersection' saved (9725 polygons)
[5/5] Computing crossing density per ward…
✓ Layer 'ward_crossing_density_wgs84' saved (23 wards)
============================================================
SUMMARY
============================================================
Wards: 23
Crossing points: 5462
Buffers (50 m): 5462
Density table: 23 rows
Output: /work/tokyo_crossings.gpkg
Layer verification:
wards_jgd: 23 features, CRS=EPSG:6679
crossing_points: 5462 features, CRS=EPSG:6679
crossing_buffers_50m: 5462 features, CRS=EPSG:6679
buffer_ward_intersection: 9725 features, CRS=EPSG:6679
ward_crossing_density_wgs84: 23 features, CRS=EPSG:4326
Done.
26:44
Let me check if the area issue is fixed now.
26:44
Bash
python3 -c "
python3 -c "
import geopandas as gpd
import numpy as np
GPKG = '/work/tokyo_crossings.gpkg'
wards = gpd.read_file(GPKG, layer='wards_jgd')
wards_wgs = wards.to_crs('EPSG:4326')
for idx, row in wards.iterrows():
print(f'{row[\"ward_name_en\"]:20s} area={row.geometry.area/1e6:.2f} km², type={row.geometry.geom_type}')
# Check Suginami
sug = wards[wards['ward_name_en']=='Suginami Ward']
geom = sug.geometry.iloc[0]
print(f'\nSuginami geom type: {geom.geom_type}')
if geom.geom_type == 'MultiPolygon':
print(f'Num polygons: {len(geom.geoms)}')
total = sum(p.area for p in geom.geoms)
print(f'Sum of parts: {total/1e6:.2f} km²')
" 2>&1
Suginami Ward area=0.86 km², type=MultiPolygon Nakano Ward area=0.64 km², type=MultiPolygon Shinjuku Ward area=0.65 km², type=MultiPolygon Bunkyo Ward area=0.59 km², type=MultiPolygon Taito Ward area=0.26 km², type=MultiPolygon Sumida Ward area=0.26 km², type=MultiPolygon Chuo Ward area=0.12 km², type=MultiPolygon Meguro Ward area=0.44 km², type=MultiPolygon Ota Ward area=2.63 km², type=MultiPolygon Setagaya Ward area=1.39 km², type=MultiPolygon Shibuya Ward area=0.43 km², type=MultiPolygon Toshima Ward area=0.70 km², type=MultiPolygon Kita Ward area=0.79 km², type=MultiPolygon Arakawa Ward area=0.35 km², type=MultiPolygon Itabashi Ward area=0.91 km², type=MultiPolygon Nerima Ward area=1.38 km², type=MultiPolygon Adachi Ward area=1.56 km², type=MultiPolygon Minato Ward area=0.64 km², type=MultiPolygon Katsushika Ward area=2.15 km², type=MultiPolygon Chiyoda Ward area=0.14 km², type=MultiPolygon Edogawa Ward area=8.17 km², type=MultiPolygon Koto Ward area=6.07 km², type=MultiPolygon Shinagawa Ward area=1.09 km², type=MultiPolygon Suginami geom type: MultiPolygon Num polygons: 166 Sum of parts: 0.86 km²
26:57
The areas are still tiny (0.86 km² for Suginami which should be ~34 km²). The issue is that the Overpass `out body geom` is returning only the way geometries that are WAY too small - possibly the individual way segments don't form the full boundary. The relation members are ways that each represent a portion of the boundary, but they need to be combined properly.
Actually, I think the problem is that `out body geom` returns the geometry of each way member, but these are just individual way geometries (not the full boundary). When I use them to build a polygon, it's only using the first outer way which is way too small.
Wait, looking at the code more carefully - each "outer" way in the relation is a piece of the boundary. So I'm building a polygon from ALL outer way nodes... but the code creates individual polygons from each outer way rather than combining them.
Actually, looking at the code:
```python
if len(outer_rings) == 1:
poly = ShapelyPolygon(outer_rings[0], ...)
else:
# Multipolygon
polys = []
for o_ring in outer_rings:
...
polys.append(ShapelyPolygon(o_ring, holes_for_this if holes_for_this else None))
poly = MultiPolygon(polys) if len(polys) > 1 else polys[0]
```
Each outer way is a piece of the boundary line, not a closed ring! So each one is just a tiny fragment. I need to combine all outer ways' coordinates into a single ring (with proper ordering).
But wait - the Overpass `geom` output gives each way's geometry in order. So for a relation with multiple outer ways, each way's geometry is a portion of the boundary. I need to concatenate all the coordinates from all outer ways (in the order they appear in the relation) to form the complete outer ring.
Actually, for a boundary relation, the outer ways should form one or more closed rings. The issue is I'm creating individual polygons from each way rather than concatenating/merging the lines into rings.
Let me use a different approach: combine all outer way geometries into a single line (or lines), then polygonize them.
26:57
26:57
Bash
python3 -c "
python3 -c "
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Let's check one outer way for Suginami
query = '[out:json];rel(id:1543055);out body geom;'
resp = requests.post(url, data={'data': query}, timeout=60, headers=headers)
data = resp.json()
# Find the relation
for el in data.get('elements', []):
if el['type'] == 'relation' and el['id'] == 1543055:
outer_ways = [m for m in el['members'] if m['type'] == 'way' and m['role'] == 'outer']
print(f'First outer way: ref={outer_ways[0][\"ref\"]}')
print(f' Coords count: {len(outer_ways[0][\"geometry\"])}')
print(f' First coord: {outer_ways[0][\"geometry\"][0]}')
print(f' Last coord: {outer_ways[0][\"geometry\"][-1]}')
# Check if way nodes connect
# First way's last point vs second way's first point
if len(outer_ways) > 1:
w0_last = outer_ways[0]['geometry'][-1]
w1_first = outer_ways[1]['geometry'][0]
print(f'Way0 last vs Way1 first: {w0_last} vs {w1_first}')
# Total coords across all outer ways
total = sum(len(m['geometry']) for m in outer_ways if 'geometry' in m)
print(f'Total coords across all {len(outer_ways)} outer ways: {total}')
break
" 2>&1
First outer way: ref=1356066346
Coords count: 2
First coord: {'lat': 35.7097505, 'lon': 139.6497069}
Last coord: {'lat': 35.7096852, 'lon': 139.6499474}
Way0 last vs Way1 first: {'lat': 35.7096852, 'lon': 139.6499474} vs {'lat': 35.7108035, 'lon': 139.6460035}
Total coords across all 83 outer ways: 1377
27:17
I see - the outer ways are individual line segments that make up the boundary. Way0 last and Way1 first don't connect directly. The outer boundary is formed by multiple way members that together form a closed ring. I need to concatenate all the coordinates from the outer ways properly to form a single polygon, not treat each way as a separate polygon. Let me look at this more carefully. The OSM relation for a boundary typically has multiple way members (ways with role "outer") that together form the boundary ring. Each way has 2 or more nodes. The ways should connect end-to-end to form a closed ring. I need to merge the coordinates from all outer ways into a single ring by connecting them properly. Let me use shapely's linemerge or just concatenate all coords. Actually, looking at the data more carefully: each outer way is just a segment of the boundary. Some join end-to-end, some don't directly. The correct approach is to: 1. Collect all way geometries from outer members 2. Merge them using linemerge or unary_union to get the boundary line 3. Polygonize Let me use shapely's `linemerge` on all the outer way segments to create a single boundary line, then polygonize. Wait, but the way0 last point (35.7096852, 139.6499474) doesn't match way1 first point (35.7108035, 139.6460035). So the ways don't necessarily connect directly end-to-end. They should share common nodes though. Let me check if there are shared nodes. Actually, each way shares a node with the next way. The ways in the relation might not be in order. I need to order them properly, or use a more robust approach. The simplest approach: collect all the LineString geometries from the outer ways, use unary_union to merge them into a single line (or multiline), then use polygonize to create the polygon. This handles the ordering automatically. For the inner rings, same approach. Let me rewrite the polygon construction logic.
27:17
27:17
Bash
python3 -c "
python3 -c "
from shapely.ops import unary_union, linemerge, polygonize
from shapely.geometry import LineString, MultiLineString, Polygon
import requests
headers = {'User-Agent': 'OSM-GIS-Analysis/1.0'}
url = 'https://overpass-api.de/api/interpreter'
# Get Suginami outer ways
query = '[out:json];rel(id:1543055);out body geom;'
resp = requests.post(url, data={'data': query}, timeout=60, headers=headers)
data = resp.json()
for el in data.get('elements', []):
if el['type'] == 'relation' and el['id'] == 1543055:
outer_lines = []
for m in el['members']:
if m['type'] == 'way' and m['role'] == 'outer' and 'geometry' in m:
coords = [(p['lon'], p['lat']) for p in m['geometry']]
if len(coords) >= 2:
outer_lines.append(LineString(coords))
# Merge all lines
merged = unary_union(outer_lines)
print(f'Merged type: {merged.geom_type}')
if merged.geom_type == 'LineString':
# Simple case
if merged.is_ring:
poly = Polygon(merged)
print(f'Ring area: {poly.area:.0f} m² in degrees')
else:
# Try to close it
from shapely.geometry import LinearRing
coords = list(merged.coords)
if coords[0] != coords[-1]:
coords.append(coords[0])
poly = Polygon(coords)
print(f'Closed area: {poly.area:.6f} deg²')
elif merged.geom_type == 'MultiLineString':
# Try polygonize
rings = list(polygonize(merged.geoms))
print(f'Polygonized: {len(rings)} rings')
if rings:
# Take the largest polygon
largest = max(rings, key=lambda p: p.area)
print(f'Largest area: {largest.area:.6f} deg²')
break
" 2>&1
Merged type: MultiLineString Polygonized: 1 rings Largest area: 0.003379 deg²
27:52
The polygonization approach works. 0.003379 deg² in Tokyo area... Let me estimate: 1 degree of latitude ≈ 111 km, 1 degree of longitude at 35.7°N ≈ 111*cos(35.7°) ≈ 90 km. So 0.003379 deg² ≈ 0.003379 * (111*90) km² ≈ 33.8 km². That's about right for Suginami (34.39 km² according to Wikipedia). So the approach should be: 1. Collect all outer way geometries as LineStrings 2. Use `unary_union` to merge them 3. Use `polygonize` or manually create a polygon from the merged ring 4. Handle inner rings similarly Let me rewrite the entire script with this approach.
27:52
27:52
Write
/work/solve.py
content · 398 lines · py
import osmnx as ox
import geopandas as gpd
import pandas as pd
import numpy as np
import requests
import warnings
warnings.filterwarnings('ignore')
ox.settings.use_cache = True
ox.settings.log_console = False
ox.settings.timeout = 300
ox.settings.http_user_agent = "OSM-GIS-Analysis/1.0"
METRIC_CRS = "EPSG:6679" # JGD2011 / Japan Plane Rectangular CS IX (Tokyo)
WGS84_CRS = "EPSG:4326"
GPKG_PATH = "/work/tokyo_crossings.gpkg"
print("=" * 60)
print("Tokyo 23 Special Wards — Road Crossing Analysis")
print("=" * 60)
# ── Step 1: Fetch ward boundaries from Overpass ────────────────────────────
print("\n[1/5] Fetching ward boundaries from OSM…")
WARD_INFO = {
1761742: ("Chiyoda Ward", "千代田区"),
1758897: ("Chuo Ward", "中央区"),
1761717: ("Minato Ward", "港区"),
1758858: ("Shinjuku Ward", "新宿区"),
1758878: ("Bunkyo Ward", "文京区"),
1758888: ("Taito Ward", "台東区"),
1758891: ("Sumida Ward", "墨田区"),
3554015: ("Koto Ward", "江東区"),
3554304: ("Shinagawa Ward", "品川区"),
1758936: ("Meguro Ward", "目黒区"),
1758947: ("Ota Ward", "大田区"),
1759474: ("Setagaya Ward", "世田谷区"),
1759477: ("Shibuya Ward", "渋谷区"),
1543056: ("Nakano Ward", "中野区"),
1543055: ("Suginami Ward", "杉並区"),
1759506: ("Toshima Ward", "豊島区"),
1760038: ("Kita Ward", "北区"),
1760040: ("Arakawa Ward", "荒川区"),
1760078: ("Itabashi Ward", "板橋区"),
1760119: ("Nerima Ward", "練馬区"),
1760124: ("Adachi Ward", "足立区"),
1761718: ("Katsushika Ward","葛飾区"),
1761743: ("Edogawa Ward", "江戸川区"),
}
rel_ids = list(WARD_INFO.keys())
headers_ok = {"User-Agent": "OSM-GIS-Analysis/1.0"}
overpass_url = "https://overpass-api.de/api/interpreter"
rel_id_str = ",".join(str(i) for i in rel_ids)
query_wards = f"""[out:json][timeout:300];
rel(id:{rel_id_str});
out body geom;
"""
resp = requests.post(overpass_url, data={"data": query_wards}, timeout=300, headers=headers_ok)
resp.raise_for_status()
data = resp.json()
from shapely.geometry import Polygon, MultiPolygon, LineString, MultiLineString
from shapely.ops import unary_union, polygonize
from shapely.validation import make_valid
def build_polygon_from_ways(outer_way_geom_list, inner_way_geom_list):
"""Build a Polygon from lists of way geometries (boundary segments)."""
# Merge outer lines
if not outer_way_geom_list:
return None
outer_merged = unary_union(outer_way_geom_list)
# Polygonize the merged outer boundary
if outer_merged.geom_type == 'LineString':
lines = [outer_merged]
elif outer_merged.geom_type == 'MultiLineString':
lines = list(outer_merged.geoms)
else:
return None
rings = list(polygonize(lines))
if not rings:
# Fallback: try to close manually
coords = list(outer_merged.coords)
if len(coords) >= 3:
if coords[0] != coords[-1]:
coords.append(coords[0])
rings = [Polygon(coords)]
if not rings:
return None
# Take the largest polygon as the outer boundary
outer_poly = max(rings, key=lambda p: p.area)
# Process inner rings
if inner_way_geom_list:
inner_merged = unary_union(inner_way_geom_list)
if inner_merged.geom_type == 'LineString':
inner_lines = [inner_merged]
elif inner_merged.geom_type == 'MultiLineString':
inner_lines = list(inner_merged.geoms)
else:
inner_lines = []
inner_rings_list = list(polygonize(inner_lines))
valid_holes = []
for ir in inner_rings_list:
if outer_poly.contains(ir.representative_point()):
valid_holes.append(ir.exterior)
if valid_holes:
outer_poly = Polygon(outer_poly.exterior, valid_holes)
if not outer_poly.is_valid:
outer_poly = make_valid(outer_poly)
return outer_poly
ward_records = []
for el in data.get("elements", []):
if el["type"] != "relation":
continue
rid = el["id"]
name_en, name_ja = WARD_INFO[rid]
members = el.get("members", [])
outer_lines = []
inner_lines = []
for m in members:
if m.get("type") != "way":
continue
if "geometry" not in m:
continue
coords = [(p["lon"], p["lat"]) for p in m["geometry"]]
if len(coords) < 2:
continue
line = LineString(coords)
if m.get("role") == "outer":
outer_lines.append(line)
elif m.get("role") == "inner":
inner_lines.append(line)
if not outer_lines:
print(f" ⚠ No outer ways for relation {rid}, skipping")
continue
poly = build_polygon_from_ways(outer_lines, inner_lines)
if poly is None:
print(f" ⚠ Could not build polygon for relation {rid}, skipping")
continue
ward_records.append({
"ward_id": rid,
"ward_name_en": name_en,
"ward_name": name_ja,
"geometry": poly
})
wards_gdf = gpd.GeoDataFrame(ward_records, crs=WGS84_CRS)
wards_metric = wards_gdf.to_crs(METRIC_CRS)
wards_metric.to_file(GPKG_PATH, layer="wards_jgd", driver="GPKG")
print(f" ✓ Layer 'wards_jgd' saved ({len(wards_gdf)} wards)")
# Quick area sanity check
print("\n Ward area sanity check:")
for idx, row in wards_metric.iterrows():
area_km2 = row.geometry.area / 1e6
print(f" {row['ward_name_en']:20s} {area_km2:8.2f} km²")
# ── Step 2: Fetch drivable road network ────────────────────────────────────
print("\n[2/5] Fetching drivable road network from OSM…")
# Get bounding box
def get_geom_coords(geom):
if geom.geom_type == 'MultiPolygon':
for poly in geom.geoms:
yield from poly.exterior.coords
else:
yield from geom.exterior.coords
all_lats = [pt[1] for g in wards_gdf.geometry for pt in get_geom_coords(g)]
all_lons = [pt[0] for g in wards_gdf.geometry for pt in get_geom_coords(g)]
south, north = min(all_lats), max(all_lats)
west, east = min(all_lons), max(all_lons)
print(f" Bounding box: {south:.4f}, {north:.4f}, {west:.4f}, {east:.4f}")
query_roads = f"""[out:json][timeout:600];
(
way["highway"~"^(motorway|trunk|primary|secondary|tertiary|unclassified|residential|motorway_link|trunk_link|primary_link|secondary_link|tertiary_link|living_street)$"]({south},{west},{north},{east});
);
out geom;
"""
print(" Sending Overpass query for drivable roads…")
resp = requests.post(overpass_url, data={"data": query_roads}, timeout=600, headers=headers_ok)
resp.raise_for_status()
roads_data = resp.json()
print(f" Received {len(roads_data.get('elements', []))} road elements")
road_geoms = []
for el in roads_data.get("elements", []):
if el["type"] != "way":
continue
if "geometry" not in el:
continue
coords = [(p["lon"], p["lat"]) for p in el["geometry"]]
if len(coords) < 2:
continue
road_geoms.append(LineString(coords))
if not road_geoms:
print(" WARNING: No road geometries found!")
roads_gdf = gpd.GeoDataFrame({"geometry": gpd.GeoSeries([], crs=WGS84_CRS)})
else:
roads_gdf = gpd.GeoDataFrame({"geometry": road_geoms}, crs=WGS84_CRS)
print(f" Road segments loaded: {len(roads_gdf)}")
# Project to metric and clip to combined boundary
roads_metric = roads_gdf.to_crs(METRIC_CRS)
combined_bounds = wards_metric.unary_union
roads_metric = roads_metric[roads_metric.intersects(combined_bounds)].copy()
roads_metric["geometry"] = roads_metric["geometry"].apply(
lambda g: g.intersection(combined_bounds) if g is not None else g
)
roads_metric = roads_metric[~roads_metric["geometry"].is_empty & roads_metric["geometry"].is_valid]
roads_metric = roads_metric[roads_metric.geometry.type.isin(["LineString", "MultiLineString"])]
print(f" Roads after clipping: {len(roads_metric)}")
# Flatten MultiLineStrings
simple_edges = []
for idx, row in roads_metric.iterrows():
geom = row.geometry
if geom.type == "MultiLineString":
for line in geom.geoms:
simple_edges.append(line)
else:
simple_edges.append(geom)
roads_final = gpd.GeoDataFrame({"geometry": simple_edges}, crs=METRIC_CRS)
print(f" Simplified to {len(roads_final)} single-part line segments")
# ── Step 3: Find crossing points ───────────────────────────────────────────
print("\n[3/5] Detecting road / ward-boundary crossings…")
all_crossings = []
for idx, ward in wards_metric.iterrows():
ward_id = ward["ward_id"]
ward_name_en = ward["ward_name_en"]
ward_name = ward["ward_name"]
ward_boundary = ward.geometry.boundary
candidate_mask = roads_final.intersects(ward_boundary)
candidate_roads = roads_final[candidate_mask]
for ridx, road in candidate_roads.iterrows():
road_geom = road.geometry
inter = road_geom.intersection(ward_boundary)
if inter.is_empty:
continue
if inter.type == "Point":
all_crossings.append({
"geometry": inter,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
elif inter.type == "MultiPoint":
for pt in inter.geoms:
all_crossings.append({
"geometry": pt,
"ward_id": ward_id,
"ward_name_en": ward_name_en,
"ward_name": ward_name
})
print(f" Found {len(all_crossings)} crossing points")
if len(all_crossings) > 0:
crossing_pts = gpd.GeoDataFrame(all_crossings, crs=METRIC_CRS)
else:
crossing_pts = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
crossing_pts.to_file(GPKG_PATH, layer="crossing_points", driver="GPKG")
print(f" ✓ Layer 'crossing_points' saved ({len(crossing_pts)} points)")
# ── Step 4: Create 50 m buffers ────────────────────────────────────────────
print("\n[4/5] Creating 50 m buffers and clipping to wards…")
if len(crossing_pts) > 0:
buffers = crossing_pts.copy()
buffers["geometry"] = buffers["geometry"].buffer(50.0)
buffers.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
print(f" ✓ Layer 'crossing_buffers_50m' saved ({len(buffers)} buffers)")
clipped_parts = []
for bidx, buf in buffers.iterrows():
ward_id = buf["ward_id"]
ward_geom = wards_metric.loc[wards_metric["ward_id"] == ward_id, "geometry"].iloc[0]
clipped = buf["geometry"].intersection(ward_geom)
if clipped.is_empty:
continue
if clipped.type == "MultiPolygon":
for poly in clipped.geoms:
clipped_parts.append({
"geometry": poly,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
else:
clipped_parts.append({
"geometry": clipped,
"ward_id": buf["ward_id"],
"ward_name_en": buf["ward_name_en"],
"ward_name": buf["ward_name"]
})
if clipped_parts:
clipped_gdf = gpd.GeoDataFrame(clipped_parts, crs=METRIC_CRS)
clipped_gdf["geometry"] = clipped_gdf["geometry"].apply(
lambda g: make_valid(g) if not g.is_valid else g
)
clipped_gdf = clipped_gdf[~clipped_gdf["geometry"].is_empty]
clipped_gdf.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(f" ✓ Layer 'buffer_ward_intersection' saved ({len(clipped_gdf)} polygons)")
else:
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layer 'buffer_ward_intersection' saved (empty)")
else:
empty = gpd.GeoDataFrame(
{"ward_id": pd.Series(dtype=int), "ward_name_en": pd.Series(dtype=str),
"ward_name": pd.Series(dtype=str), "geometry": pd.Series(dtype=object)},
crs=METRIC_CRS
)
empty.to_file(GPKG_PATH, layer="crossing_buffers_50m", driver="GPKG")
empty.to_file(GPKG_PATH, layer="buffer_ward_intersection", driver="GPKG")
print(" ✓ Layers 'crossing_buffers_50m' and 'buffer_ward_intersection' saved (empty)")
# ── Step 5: Ward crossing density table ────────────────────────────────────
print("\n[5/5] Computing crossing density per ward…")
if len(crossing_pts) > 0:
counts = crossing_pts.groupby("ward_id").size().reset_index(name="crossing_count")
else:
counts = pd.DataFrame({"ward_id": wards_gdf["ward_id"], "crossing_count": 0})
wards_metric["ward_area_km2"] = (wards_metric.geometry.area / 1_000_000.0).round(4)
density_gdf = wards_gdf[["ward_id", "ward_name_en", "ward_name", "geometry"]].copy()
density_gdf = density_gdf.merge(counts, on="ward_id", how="left")
density_gdf["crossing_count"] = density_gdf["crossing_count"].fillna(0).astype(int)
area_map = dict(zip(wards_metric["ward_id"], wards_metric["ward_area_km2"]))
density_gdf["ward_area_km2"] = density_gdf["ward_id"].map(area_map)
density_gdf["crossings_per_km2"] = np.where(
density_gdf["ward_area_km2"] > 0,
(density_gdf["crossing_count"] / density_gdf["ward_area_km2"]).round(4),
0.0
)
density_gdf.to_file(GPKG_PATH, layer="ward_crossing_density_wgs84", driver="GPKG")
print(f" ✓ Layer 'ward_crossing_density_wgs84' saved ({len(density_gdf)} wards)")
# ── Verify & Summary ───────────────────────────────────────────────────────
print("\n" + "=" * 60)
print("SUMMARY")
print("=" * 60)
print(f" Wards: {len(wards_gdf)}")
print(f" Crossing points: {len(crossing_pts)}")
print(f" Buffers (50 m): {len(crossing_pts)}")
print(f" Density table: {len(density_gdf)} rows")
print(f"\nOutput: {GPKG_PATH}")
import fiona
layers = fiona.listlayers(GPKG_PATH)
for lyr in layers:
gdf = gpd.read_file(GPKG_PATH, layer=lyr)
print(f" {lyr}: {len(gdf)} features, CRS={gdf.crs}")
print("\nDone.")Wrote /work/solve.py (14385 chars).