fio-l3-vienna-geofabrik-highways

Rerun

Let me break down this task:

1. Download Austria PBF from Geofabrik
2. Find the Gürtel highway ways in Vienna (name ends in 'Gürtel')
3. Create a 500m buffer around those ways
4. Extract all highways that intersect that buffer
5. Extract all public-transport route relations (bus, tram, subway, train, trolleybus, light_rail) that cross the same buffer
6. Write both layers to vienna_network.gpkg
   - Highways as LineString with attributes: name, highway, maxspeed, lanes, surface, oneway
   - PT routes as MultiLineString with attributes: ref, name, operator, route
7. Use Austria's standard projected coordinate system (EPSG:31287 - MGI / Austria Lambert, or EPSG:31258? Let me think... Austria's standard is MGI / Austria Lambert - EPSG:31287 or sometimes EPSG:31258 (MGI / Austria GK East). Actually the standard projected coordinate system for Austria is typically EPSG:31287 (MGI / Austria Lambert) or EPSG:31284-31286 (GK zones). Let me use EPSG:31287 - MGI / Austria Lambert, which is the official Austrian projection.

Actually, let me think more carefully. The standard projected coordinate system for Austria is often EPSG:31287 (MGI / Austria Lambert). For Vienna, EPSG:31285 (MGI / Austria GK Central) also works. But EPSG:31287 is the lambert conformal conic projection covering all of Austria. Let me go with EPSG:31287.

Wait, let me verify - actually "MGI / Austria Lambert" is EPSG:31287. That's the standard.

Let me plan the Python script:

1. Download the PBF file
2. Use osmnx or geopandas with osmium/pyogrio to read and process
3. Actually, the best approach might be to use osmnx for this, or directly use osm2pgsql or pyosmium, or use geopandas with the OSM PBF.

Actually, let me think about the tools available. We have geopandas, shapely, pyproj, duckdb. We might need to install osmnx or use pyosmium or osmium.

Let me think about the most efficient approach:

Option 1: Use osmnx - it can download OSM data directly. But we need to work from a local PBF.

Option 2: Use pyosmium to read the PBF and filter.

Option 3: Use ogr2ogr (GDAL) to filter the PBF.

Option 4: Use osmium command-line tool.

Let me check what's available and what's easiest. I think using osmnx might be the cleanest approach - we can use osmnx.graph_from_place or osmnx.features_from_place.

Actually, let me use a combination approach:
1. Download the PBF
2. Use osmnx to read and process it, or use pyosmium

Let me check what's installed first.

Actually, I think the most straightforward approach would be:
1. Download the PBF
2. Use osmnx to work with OSM data
3. Or use pyogrio/pyosmium to read the PBF

Let me start by checking what's installed and then write the script.

Let me think about this more carefully:

For highways:
- Load the PBF, filter for ways that have highway tag and are in/near Vienna
- Find ways where name ends with 'Gürtel'
- Buffer those by 500m
- Find all ways with highway tag that intersect the buffer
- Keep attributes: name, highway, maxspeed, lanes, surface, oneway

For PT routes:
- Find route relations (type=route, route=bus|tram|subway|train|trolleybus|light_rail)
- That cross the same buffer
- Each route as a single MultiLineString feature
- Keep attributes: ref, name, operator, route

Let me write a comprehensive Python script. I'll use osmnx since it handles OSM data well, or I could use the osmium library.

Actually, let me check what's available first.

Good, all the needed libraries are available. Let me plan the approach:

1. Download the Austria OSM PBF file
2. Use osmium or osmnx to extract:
   - Highway ways whose name ends with "Gürtel" (this is Vienna's Gürtel)
   - Buffer those by 500m
   - Extract all highways within that buffer
   - Extract PT route relations (bus, tram, subway, train, trolleybus, light_rail) that cross the buffer
3. Write to GeoPackage with proper CRS (Austria's standard projected CRS is EPSG:31287 - MGI / Austria Lambert, or EPSG:31256 - MGI / Austria GK East? Actually, the standard projected coordinate system for Austria is EPSG:31287 (MGI / Austria Lambert) or sometimes EPSG:31256/31257/31258 for Gauss-Krüger zones. I think the most common "standard" is EPSG:31287 (MGI / Austria Lambert). Let me check.

Actually, Austria's official projected CRS for nationwide use is EPSG:31287 (MGI / Austria Lambert). For Vienna specifically, EPSG:31256 (MGI / Austria GK East) is also common, but the question says "Austria's standard projected coordinate system" which is EPSG:31287.

Let me think about how to do this efficiently with the PBF file.

Approach:
1. Download the PBF
2. Use osmium to parse and filter
3. Or use osmnx which can work with OSM data

Actually, osmnx can work with local .osm files but not directly with .osm.pbf. Let me use a combination approach.

Let me think more carefully:

Step 1: Download the PBF
Step 2: Find all ways with name ending in "Gürtel" - these are the Gürtel ways
Step 3: Get the geometries of these ways, buffer by 500m
Step 4: Extract all highway ways within that buffer
Step 5: Extract all PT route relations that intersect the buffer
Step 6: Write to GeoPackage

I'll use a Python script with Pyosmium or osmnx or geopandas.

Actually, let me think about the best tools:
- `osmium` command line tool is available
- `pyosmium` (Python bindings) is available
- `osmnx` is available
- `duckdb` with Spatial extension could also work with OSM PBF

Let me try using `osmnx` to download the data - it can fetch from the OSM API or use local files. Actually osmnx can use `.osm` XML files but PBF needs to be handled differently.

Let me use a different approach:
1. Use `osmium` command line to extract the data
2. Or use `pyosmium` in Python

Actually, the cleanest approach would be:
1. Download the PBF
2. Use osmium to filter out the relevant data
3. Process with geopandas/shapely

Let me start by downloading and then writing a comprehensive Python script.

Actually, let me think about this differently. Since we need:
- Highway ways in a buffer around Gürtel
- PT route relations crossing that buffer

This is complex. Let me use osmnx's graph_from_place or similar, but we need to work with a PBF file.

Actually, osmnx can work with local .osm XML files using the `graph_from_xml` function. But handling PBF...

Let me use `osmium` CLI tools to:
1. First get the Gürtel ways
2. Then filter everything in a buffer

Or even better, I can use the `osmium` command-line tool which is very efficient for PBF processing.

Let me write a comprehensive Python script that:
1. Downloads the PBF
2. Uses pyosmium to iterate through the data
3. Processes it

Actually, let me try a different approach. I'll use `osmnx` to:
1. Get the Gürtel from OSM by name query
2. Buffer it
3. Get all highways within the buffer
4. Get PT routes

But osmnx works with the Overpass API, which might be too large for this query. The PBF file is the proper approach.

Let me use `osmium` command line tool. Here's my plan:

1. Download the PBF
2. Use `osmium tags-filter` to extract ways with name=*Gürtel and highway=* 
3. Find the Gürtel ways, buffer them
4. Then extract everything within that buffer

Actually, let me write a Python script that uses pyosmium to handle this efficiently.

Wait, I can use duckdb's spatial extension which can read OSM PBF files! Let me check...

Actually duckdb can read .osm.pbf files with the spatial extension. Let me try that approach.

No, let me just use a straightforward approach:

1. Download the PBF
2. Use osmium to extract just Vienna area first (to keep file manageable)
3. Process with Python

Let me start coding.

Actually, the simplest approach might be:
1. Download the PBF
2. Use `osmium` to extract all objects with highway=* that are in or near Vienna
3. Use Python (geopandas/shapely) to do the spatial filtering

Let me write a script that does everything in one go using pyosmium (the Python bindings).

Actually, let me think about the most efficient approach:

Step 1: Download the PBF
Step 2: Use osmium to get all ways with name ending in "Gürtel" and highway tag
Step 3: Get their coordinates, buffer by 500m
Step 4: Use osmium to extract all ways within that bounding box (plus some margin)
Step 5: Spatially filter to the actual buffer
Step 6: Do similarly for PT route relations

Let me write a Python script that does this step by step.

Actually, I think the best approach is to write a comprehensive Python script.

Let me start:

```python
import os
import urllib.request
import geopandas as gpd
import shapely
import pyproj
from shapely.geometry import shape, Point, LineString, MultiLineString, box
from shapely.ops import unary_union
import osmium
import json
from collections import defaultdict
```

Hmm, but pyosmium's API is a bit complex for extracting relations. Let me think about a cleaner approach.

Actually, let me use a combination of tools:
1. Download the PBF
2. Use `osmium` CLI to filter the data
3. Use geopandas for processing

Let me start with downloading.

Actually, let me try a more direct approach. I'll use `osmnx` to load the graph from the PBF file... no, osmnx doesn't support PBF directly.

Let me try yet another approach: use `pyosmium` to process the PBF and extract the needed features into GeoJSON, then use geopandas to write to GeoPackage.

Here's my plan:
1. Download the PBF
2. Write a Python script that uses pyosmium to:
   a. First pass: find all ways with name ending in "Gürtel" and highway tag, collect their node references
   b. Collect all node coordinates
   c. Buffer the Gürtel geometries
   d. Second pass: collect all highway ways that intersect the buffer
   e. Also collect PT route relations
3. Write to GeoPackage

Let me code this up.

Actually, I think the simplest approach is to use `osmium` command line tools for filtering, and then process the results.

Let me try:

```bash
# Download PBF
wget https://download.geofabrik.de/europe/austria-latest.osm.pbf -O /work/austria-latest.osm.pbf

# Extract Vienna area using a bounding box
osmium extract -b 16.1,48.1,16.6,48.35 /work/austria-latest.osm.pbf -o /work/vienna.osm.pbf

# Get all ways with name=*Gürtel and highway=*
osmium tags-filter /work/vienna.osm.pbf w/highway w/name=Gürtel -o /work/guertel_ways.osm.pbf
```

Hmm, but osmium tags-filter has limitations on matching name ending with "Gürtel".

Let me just write a Python script using pyosmium. It's the most flexible.

Actually, let me try yet another approach: use `osmnx` with the `graph_from_place` function but that uses Overpass API which might be slow for large areas...

No, let me just write a comprehensive Python script using pyosmium. Let me do it.

Here's my plan for the script:

```python
#!/usr/bin/env python3
"""
Extract highways and PT routes around Vienna's Gürtel.
"""

import os
import sys
import urllib.request
import geopandas as gpd
import numpy as np
import shapely
import pyproj
from shapely.geometry import shape, Point, LineString, MultiLineString, box
from shapely.ops import unary_union
import osmium
from collections import defaultdict
import math

# ... etc
```

Let me write the full script now. I'll break it into steps.

Step 1: Download the PBF
Step 2: First pass through the PBF to find Gürtel ways and their geometries
Step 3: Buffer the Gürtel by 500m
Step 4: Second pass to find all highway ways in the buffer
Step 5: Third pass to find PT route relations crossing the buffer
Step 6: Write to GeoPackage

Actually, for efficiency, I can do multiple things in fewer passes. Since pyosmium streams through the file, I can:

Pass 1: Collect all nodes, find Gürtel ways, build Gürtel geometries
Pass 2 (or continuing): Find all highway ways, check if they intersect buffer, collect them
Pass 3: Find PT route relations

But the issue is that relations reference ways, and ways reference nodes. So I need multiple passes or I need to store everything.

Let me think about a more efficient approach:

Pass 1: 
- Collect all nodes (store id -> (lon, lat))
- Find ways with highway tag (store id -> [node_ids])
- Among those, find ones with name ending in "Gürtel"
- Build geometries for Gürtel ways
- Buffer Gürtel

Pass 2 (after buffer is known):
- Find highway ways that intersect the buffer
- Find PT route relations and their member ways

But for relations, I need the way geometries first. So maybe:

Pass 1: Collect all nodes
Pass 2: Collect all ways (build linestrings from nodes), identify Gürtel, buffer
Pass 3: Filter ways by buffer, collect relations and check if they cross buffer

This is getting complex. Let me just write it step by step.

Actually, let me simplify. I'll use two approaches:
1. Use osmium CLI to extract a bbox around Vienna
2. Use Python/pyosmium to process

Let me start coding now.

Actually, I just realized I should use `osmnx` more creatively. osmnx can download data from Overpass API and work with it locally. But for large datasets like all highways in Vienna, it might be slow.

Let me just go with the pyosmium approach. Here's my refined plan:

```python
import osmium
import geopandas as gpd
from shapely.geometry import LineString, MultiLineString, Point, box
from shapely.ops import unary_union
import pyproj
import json
from collections import defaultdict
import math

# Step 1: Download
# ...

# Step 2: Handler class for osmium
class OSMHandler(osmium.SimpleHandler):
    def __init__(self):
        super().__init__()
        self.nodes = {}  # id -> (lon, lat)
        self.ways = {}   # id -> way info
        self.relations = {}  # id -> relation info
        self.guertel_way_ids = set()
        self.guertel_geoms = []
        self.buffer_geom = None
        
    def node(self, n):
        self.nodes[n.id] = (n.location.lon, n.location.lat)
        
    def way(self, w):
        tags = dict(w.tags)
        if 'highway' in tags:
            name = tags.get('name', '')
            node_ids = [n.ref for n in w.nodes]
            self.ways[w.id] = {
                'tags': tags,
                'nodes': node_ids
            }
            if name.endswith('Gürtel') or name.endswith('Gürtel'):
                self.guertel_way_ids.add(w.id)
                
    def relation(self, r):
        tags = dict(r.tags)
        route = tags.get('route', '')
        if route in ('bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail'):
            self.relations[r.id] = {
                'tags': tags,
                'members': [(m.type, m.ref, m.role) for m in r.members]
            }
```

This is getting complex. Let me just write it all out properly.

Actually, I wonder if there's a simpler way. Let me check if `osmnx` can load from PBF... No, it can't directly.

What about using `osmium` CLI to convert the PBF to GeoJSON or another format? osmium can output GeoJSON but only for simple cases.

Let me use `osmium` CLI with the `-F pbf` flag to filter, then convert to something I can read with geopandas.

Actually, the cleanest approach might be:

1. Download the PBF
2. Use `osmium extract` with a bbox around Vienna to get a smaller PBF
3. Use `osmium export` to convert to GeoJSON
4. But this won't handle relations well...

Let me just write a comprehensive Python script. Let me do it now.

I'll write the script to /work/solve.py and run it.

Key considerations:
- The PBF is large (Austria is ~1GB)
- Need to process efficiently
- Need to handle UTF-8 (Gürtel has ü)

Let me write the script:

```python
#!/usr/bin/env python3
"""
Extract highway and PT route data around Vienna's Gürtel from OSM PBF.
"""

import os
import sys
import urllib.request
import osmium
import geopandas as gpd
import shapely
import pyproj
from shapely.geometry import LineString, MultiLineString, Point, box, shape
from shapely.ops import unary_union
from collections import defaultdict
import json

WORKDIR = '/work'
PBF_URL = 'https://download.geofabrik.de/europe/austria-latest.osm.pbf'
PBF_FILE = os.path.join(WORKDIR, 'austria-latest.osm.pbf')
GPKG_FILE = os.path.join(WORKDIR, 'vienna_network.gpkg')

# Target CRS: EPSG:31287 (MGI / Austria Lambert) - Austria's standard projected CRS
TARGET_CRS = 'EPSG:31287'

# Download PBF if not exists
if not os.path.exists(PBF_FILE):
    print("Downloading Austria PBF...")
    urllib.request.urlretrieve(PBF_URL, PBF_FILE)
    print("Download complete.")

# Pass 1: Collect all nodes
print("Pass 1: Collecting nodes...")
nodes = {}
class NodeHandler(osmium.SimpleHandler):
    def node(self, n):
        nodes[n.id] = (n.location.lon, n.location.lat)

handler = NodeHandler()
handler.apply_file(PBF_FILE)
print(f"  Collected {len(nodes)} nodes.")

# Pass 2: Find Gürtel ways and build their geometries
print("Pass 2: Finding Gürtel ways...")
guertel_ways = []  # list of (way_id, tags, linestring)
all_highway_ways = {}  # way_id -> {tags, geometry}

class WayHandler(osmium.SimpleHandler):
    def way(self, w):
        tags = dict(w.tags)
        if 'highway' not in tags:
            return
        name = tags.get('name', '')
        node_refs = [n.ref for n in w.nodes if n.ref in nodes]
        if len(node_refs) < 2:
            return
        coords = [(nodes[nid][0], nodes[nid][1]) for nid in node_refs if nid in nodes]
        if len(coords) < 2:
            return
        geom = LineString(coords)
        is_guertel = 'Gürtel' in name and name.endswith('Gürtel')
        if is_guertel:
            guertel_ways.append((w.id, tags, geom))
        all_highway_ways[w.id] = {'tags': tags, 'geometry': geom}

handler = WayHandler()
handler.apply_file(PBF_FILE)
print(f"  Found {len(guertel_ways)} Gürtel ways, {len(all_highway_ways)} total highway ways.")

# Build Gürtel buffer
print("Building Gürtel buffer...")
if guertel_ways:
    guertel_geom = unary_union([gw[2] for gw in guertel_ways])
    # Project to a local projection for buffering
    gdf_guertel = gpd.GeoDataFrame({'geometry': [guertel_geom]}, crs='EPSG:4326')
    gdf_guertel = gdf_guertel.to_crs(TARGET_CRS)
    buffer_geom = gdf_guertel.buffer(500).iloc[0]
    print(f"  Buffer area: {buffer_geom.area:.0f} m²")
else:
    print("  ERROR: No Gürtel ways found!")
    sys.exit(1)

# Reproject buffer back to WGS84 for intersection checks
gdf_buffer_4326 = gpd.GeoDataFrame({'geometry': [buffer_geom]}, crs=TARGET_CRS).to_crs('EPSG:4326')
buffer_geom_4326 = gdf_buffer_4326.iloc[0].geometry

# Filter highway ways within buffer
print("Filtering highway ways within buffer...")
highway_features = []
for wid, info in all_highway_ways.items():
    geom = info['geometry']
    if geom.intersects(buffer_geom_4326):
        tags = info['tags']
        highway_features.append({
            'geometry': geom,
            'name': tags.get('name', None),
            'highway': tags.get('highway', None),
            'maxspeed': tags.get('maxspeed', None),
            'lanes': tags.get('lanes', None),
            'surface': tags.get('surface', None),
            'oneway': tags.get('oneway', None),
        })

print(f"  Found {len(highway_features)} highway features in buffer.")

# Pass 3: Find PT route relations
# We need to find relations with route=bus|tram|subway|train|trolleybus|light_rail
# Then check if they cross the buffer (by checking their member ways)
print("Pass 3: Finding PT route relations...")

# First, collect relations
pt_route_types = {'bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail'}
pt_relations = []

class RelationHandler(osmium.SimpleHandler):
    def relation(self, r):
        tags = dict(r.tags)
        if tags.get('type') == 'route' and tags.get('route') in pt_route_types:
            members = [(m.type, m.ref, m.role) for m in r.members]
            pt_relations.append({
                'id': r.id,
                'tags': tags,
                'members': members
            })

handler = RelationHandler()
handler.apply_file(PBF_FILE)
print(f"  Found {len(pt_relations)} PT route relations.")

# Now, for each PT relation, collect the member ways that are within the buffer
# and build the route geometry
pt_features = []
for rel in pt_relations:
    member_ways = []
    for mtype, mref, mrole in rel['members']:
        if mtype == 'w' and mref in all_highway_ways:
            geom = all_highway_ways[mref]['geometry']
            if geom.intersects(buffer_geom_4326):
                member_ways.append(geom)
    
    if member_ways:
        # Merge into MultiLineString
        merged = unary_union(member_ways)
        if merged.geom_type == 'LineString':
            merged = MultiLineString([merged])
        elif merged.geom_type == 'MultiLineString':
            pass  # already good
        else:
            # Might be GeometryCollection or something else
            lines = [g for g in merged.geoms if g.geom_type in ('LineString', 'MultiLineString')]
            if lines:
                all_lines = []
                for l in lines:
                    if l.geom_type == 'LineString':
                        all_lines.append(l)
                    else:
                        all_lines.extend(l.geoms)
                merged = MultiLineString(all_lines)
            else:
                continue
        
        tags = rel['tags']
        pt_features.append({
            'geometry': merged,
            'ref': tags.get('ref', None),
            'name': tags.get('name', None),
            'operator': tags.get('operator', None),
            'route': tags.get('route', None),
        })

print(f"  Found {len(pt_features)} PT route features crossing the buffer.")

# Create GeoDataFrames and write to GeoPackage
print("Writing to GeoPackage...")

# Highways layer
gdf_highways = gpd.GeoDataFrame(highway_features, crs='EPSG:4326')
gdf_highways = gdf_highways.to_crs(TARGET_CRS)
print(f"  Highways: {len(gdf_highways)} features")

# PT routes layer
gdf_pt = gpd.GeoDataFrame(pt_features, crs='EPSG:4326')
gdf_pt = gdf_pt.to_crs(TARGET_CRS)
print(f"  PT routes: {len(gdf_pt)} features")

# Write to GeoPackage
if os.path.exists(GPKG_FILE):
    os.remove(GPKG_FILE)

gdf_highways.to_file(GPKG_FILE, layer='highways', driver='GPKG')
gdf_pt.to_file(GPKG_FILE, layer='pt_routes', driver='GPKG')

print(f"Done! Output: {GPKG_FILE}")
```

Wait, there's an issue with the unary_union of LineStrings - if they don't connect, it creates a MultiLineString automatically. But if they do connect, it merges them into a single LineString. Let me handle that properly.

Also, I need to make sure the PT routes are MultiLineString type as specified.

Let me also think about performance - this script processes the entire Austria PBF (which could be >1GB) three times. That could be slow. Let me optimize.

Actually, let me restructure to do it in fewer passes:

Pass 1: Collect nodes AND find Gürtel ways AND find all highway ways AND find PT relations
- Store nodes
- Store highway ways (with their node refs, not geometries yet)
- Store PT relations (with their member references)
- Identify Gürtel ways

Then build geometries for highway ways using the nodes.

Then check which highway ways are in the buffer.

Then for PT relations, check if any of their member ways are in the buffer.

This is 1 pass instead of 3, which is much faster.

Let me rewrite the script.

Actually, storing all highway ways in memory could be large. But let's try it - Austria's highway network should be manageable.

Let me write the optimized script.

One more thing - I need to handle the case where `name` ends with 'Gürtel'. The 'ü' character could be 'u' with umlaut or could be written as 'ue'. Let me check what's in OSM data. In OSM, street names in Vienna use 'ü' (the actual UTF-8 character). So `name.endswith('Gürtel')` or `name.endswith('ürtel')` - actually the full word is "Gürtel" so checking endswith('Gürtel') or name.stripped endswith 'Gürtel' should work.

Wait, the user says "name ends in 'Gürtel'" - so the full name ends with the word "Gürtel". Let me check: in Vienna, the street is called "Gürtel" (e.g., "Wiener Gürtel", "Innerer Gürtel", "Äußerer Gürtel"). So the name field would be something like "Wiener Gürtel" or just "Gürtel Straße"... Actually in OSM, the main ring road is tagged as "Gürtel" or with various prefixes.

Let me just check with `name.endswith('Gürtel')` - this should catch names like "Wiener Gürtel", "Äußerer Gürtel Straße"... wait, "Äußerer Gürtel Straße" ends with "Straße", not "Gürtel". Hmm.

Actually, in OSM, Vienna's Gürtel is typically tagged with names like:
- "Gürtel" (simple)
- "Wiener Gürtel" 
- "Äußerer Gürtel"
- etc.

But some segments might have "Gürtel" in the middle of the name, like "Gürtel Straße" - no, "Gürtel Straße" ends with "Straße" not "Gürtel".

Let me use a more flexible condition: `'Gürtel' in name` or `name.endswith('Gürtel')`. Actually the user specifically says "name ends in 'Gürtel'" so let me do exactly that. But I should also consider that some might have "Gürtel" as the full name (exact match). Let me check:

`name.endswith('Gürtel')` would match:
- "Gürtel" → True
- "Wiener Gürtel" → True
- "Gürtel Straße" → False

Let me use `name.endswith('Gürtel') or name == 'Gürtel'` which is equivalent to just `name.endswith('Gürtel')` since 'Gürtel'.endswith('Gürtel') is True.

OK let me write the complete script now.

I also need to think about the possibility that nodes might not exist in the nodes dict for some way refs (if nodes are outside the bounding box or something). In a full country PBF, all referenced nodes should be present.

Let me also think about the size - Austria's PBF is about 1.6GB as of recent. Processing it might take a while but should be fine.

Let me write the script:

```python
#!/usr/bin/env python3
"""
Extract highway and PT route data around Vienna's Gürtel.
"""

import os
import sys
import urllib.request
import osmium
import geopandas as gpd
from shapely.geometry import LineString, MultiLineString, shape
from shapely.ops import unary_union
from collections import defaultdict

WORKDIR = '/work'
PBF_URL = 'https://download.geofabrik.de/europe/austria-latest.osm.pbf'
PBF_FILE = os.path.join(WORKDIR, 'austria-latest.osm.pbf')
GPKG_FILE = os.path.join(WORKDIR, 'vienna_network.gpkg')
TARGET_CRS = 'EPSG:31287'  # MGI / Austria Lambert

PT_ROUTE_TYPES = {'bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail'}

# ---- Download ----
if not os.path.exists(PBF_FILE):
    print("Downloading Austria PBF...", flush=True)
    urllib.request.urlretrieve(PBF_URL, PBF_FILE)
    print("Download complete.", flush=True)

# ---- Single-pass processing ----
print("Processing PBF...", flush=True)

nodes = {}
highway_ways = {}  # way_id -> {tags, node_refs}
pt_relations = []  # list of {id, tags, members: [(type, ref, role)]}
guertel_way_ids = set()

class Handler(osmium.SimpleHandler):
    def node(self, n):
        if n.location.valid():
            nodes[n.id] = (n.location.lon, n.location.lat)
    
    def way(self, w):
        tags = dict(w.tags)
        if 'highway' not in tags:
            return
        node_refs = [n.ref for n in w.nodes]
        highway_ways[w.id] = {
            'tags': tags,
            'nodes': node_refs
        }
        name = tags.get('name', '')
        if name.endswith('Gürtel'):
            guertel_way_ids.add(w.id)
    
    def relation(self, r):
        tags = dict(r.tags)
        if tags.get('type') == 'route' and tags.get('route') in PT_ROUTE_TYPES:
            members = [(m.type, m.ref, m.role) for m in r.members]
            pt_relations.append({
                'id': r.id,
                'tags': tags,
                'members': members
            })

handler = Handler()
handler.apply_file(PBF_FILE)

print(f"Nodes: {len(nodes)}", flush=True)
print(f"Highway ways: {len(highway_ways)}", flush=True)
print(f"Gürtel ways: {len(guertel_way_ids)}", flush=True)
print(f"PT relations: {len(pt_relations)}", flush=True)

# ---- Build Gürtel geometries ----
print("Building Gürtel geometries...", flush=True)
guertel_geoms = []
for wid in guertel_way_ids:
    info = highway_ways[wid]
    coords = []
    for nid in info['nodes']:
        if nid in nodes:
            coords.append(nodes[nid])
    if len(coords) >= 2:
        guertel_geoms.append(LineString(coords))

if not guertel_geoms:
    print("ERROR: No Gürtel ways found!", flush=True)
    sys.exit(1)

guertel_union = unary_union(guertel_geoms)
gdf_guertel = gpd.GeoDataFrame({'geometry': [guertel_union]}, crs='EPSG:4326')
gdf_guertel_proj = gdf_guertel.to_crs(TARGET_CRS)
buffer_proj = gdf_guertel_proj.buffer(500).iloc[0]
print(f"Buffer area: {buffer_proj.area:.0f} m²", flush=True)

# Convert buffer back to WGS84 for intersection with WGS84 geometries
gdf_buffer_4326 = gpd.GeoDataFrame({'geometry': [buffer_proj]}, crs=TARGET_CRS).to_crs('EPSG:4326')
buffer_geom = gdf_buffer_4326.iloc[0].geometry

# ---- Filter highway ways ----
print("Filtering highway ways in buffer...", flush=True)
highway_features = []
for wid, info in highway_ways.items():
    coords = []
    for nid in info['nodes']:
        if nid in nodes:
            coords.append(nodes[nid])
    if len(coords) < 2:
        continue
    geom = LineString(coords)
    if geom.intersects(buffer_geom):
        tags = info['tags']
        highway_features.append({
            'geometry': geom,
            'name': tags.get('name'),
            'highway': tags.get('highway'),
            'maxspeed': tags.get('maxspeed'),
            'lanes': tags.get('lanes'),
            'surface': tags.get('surface'),
            'oneway': tags.get('oneway'),
        })

print(f"Highway features in buffer: {len(highway_features)}", flush=True)

# ---- Build PT route geometries ----
print("Building PT route geometries...", flush=True)
pt_features = []
for rel in pt_relations:
    member_lines = []
    for mtype, mref, mrole in rel['members']:
        if mtype == 'w' and mref in highway_ways:
            info = highway_ways[mref]
            coords = []
            for nid in info['nodes']:
                if nid in nodes:
                    coords.append(nodes[nid])
            if len(coords) >= 2:
                geom = LineString(coords)
                if geom.intersects(buffer_geom):
                    member_lines.append(geom)
    
    if member_lines:
        merged = unary_union(member_lines)
        if merged.geom_type == 'LineString':
            merged = MultiLineString([merged])
        elif merged.geom_type == 'MultiLineString':
            pass  # already correct
        else:
            # Skip unexpected geometry types
            continue
        
        tags = rel['tags']
        pt_features.append({
            'geometry': merged,
            'ref': tags.get('ref'),
            'name': tags.get('name'),
            'operator': tags.get('operator'),
            'route': tags.get('route'),
        })

print(f"PT route features: {len(pt_features)}", flush=True)

# ---- Write to GeoPackage ----
print("Writing to GeoPackage...", flush=True)

gdf_highways = gpd.GeoDataFrame(highway_features, crs='EPSG:4326')
gdf_highways = gdf_highways.to_crs(TARGET_CRS)

gdf_pt = gpd.GeoDataFrame(pt_features, crs='EPSG:4326')
gdf_pt = gdf_pt.to_crs(TARGET_CRS)

if os.path.exists(GPKG_FILE):
    os.remove(GPKG_FILE)

gdf_highways.to_file(GPKG_FILE, layer='highways', driver='GPKG')
gdf_pt.to_file(GPKG_FILE, layer='pt_routes', driver='GPKG')

print("Done!", flush=True)
print(f"Output: {GPKG_FILE}", flush=True)
print(f"  Highways layer: {len(gdf_highways)} LineString features", flush=True)
print(f"  PT Routes layer: {len(gdf_pt)} MultiLineString features", flush=True)
```

This looks good. Let me write it and run it.

Wait, one issue: the script might run out of memory since we're storing all highway ways and nodes for the entire country of Austria. Let me check... Austria PBF is about 1.6GB, and storing all nodes (each as a tuple of two floats + dict overhead) plus ways could be memory-intensive but should be manageable on modern hardware with ~8GB+ RAM.

However, let me think about whether I can optimize. Actually, storing all nodes as a dict of id->(lon,lat) for Austria (roughly 200-300 million nodes?) could be a lot of memory.

Hmm, actually Austria OSM has about 300-400 million nodes. Each node stored as a Python object would be... a dict entry with an integer key and a tuple value. That's about 100+ bytes per entry = 30-40 GB. That's too much!

Let me think of a more memory-efficient approach.

Alternative 1: Process in two stages
- First pass: Just find Gürtel ways and their bounding box
- Extract just the Vienna area using osmium extract
- Process the smaller Vienna PBF

Alternative 2: Use osmium to filter first
- Use `osmium extract -b <bbox>` to get just the Vienna area
- Then process the smaller file

Let me use the extract approach. The bbox for Vienna is roughly:
- South: 48.1
- North: 48.35
- West: 16.1
- East: 16.6

But I also need some buffer since the Gürtel is near the center of Vienna. Let me use a generous bbox.

Actually, the Gürtel is a ring road around the inner city of Vienna. Vienna is roughly at 48.2°N, 16.37°E. The Gürtel is approximately within [48.18, 48.23] N and [16.33, 16.40] E. But to be safe for a 500m buffer plus all surrounding highways, let me use a larger bbox.

Let me use: south=48.05, north=48.40, west=16.10, east=16.70

Let me use osmium extract first, then process the smaller file.

```bash
osmium extract -b 16.1,48.05,16.7,48.4 austria-latest.osm.pbf -o vienna_area.osm.pbf --strategy=smart
```

This would give us a much smaller file to work with.

Let me adjust the script.

Actually, let me also reconsider - for the highway filter, I need all highway ways in a 500m buffer around the Gürtel. The Gürtel is about 12km long, so a 500m buffer around it covers a significant area of Vienna but not all of it. The bbox I chose (16.1,48.05,16.7,48.4) is generous enough.

Let me use this two-step approach:
1. Download PBF
2. Extract bbox around Vienna
3. Process the extracted file

Let me write the complete script now.

Actually, I realize I can do the extraction as a shell command before running the Python script, or incorporate it into the script. Let me just write it as a Python script that does everything.

One more thing: I need to be careful about the `osmium extract` command - the `--strategy=smart` is usually best.

Let me finalize the script.

For the relation handling, I also need to think about the fact that relation members can be ways that are not in highway_ways (e.g., they could be railway tracks for trains/subways). For bus routes, members are typically highway ways. For subway/train, members might be railway ways. Let me also collect railway ways.

Actually, the user wants "every public-transport route relation (bus, tram, subway, train, trolleybus, light_rail) that crosses the same buffer, with each route as a single feature." The route's geometry is the union of its member ways. For bus routes, members are typically highway=* ways. For tram routes, members might be highway=* or railway=* or both. For subway/train, members are typically railway=*. 

So I need to also collect ways that are not highway but could be part of PT routes. Let me cast a wider net and store all ways that could be part of these PT routes.

Actually, the simplest approach: collect ALL ways (not just highway) for the Vienna area. That way when processing relations, I can look up any member way. But this increases memory usage.

Hmm, let me think... For the Vienna area PBF (extracted with generous bbox), the number of ways should be manageable. Let me just collect all ways.

Let me rewrite:

```python
# Store all ways, not just highway, to handle PT route members
all_ways = {}  # way_id -> {tags, nodes}
```

Then when filtering for highways, I only add the highway ones to the output. For PT routes, I use all_ways to look up member geometries.

This is cleaner.

Let me finalize the full script now.

Wait, there's another consideration. For route relations, the members might include not just ways but also nodes (for stops). I should skip non-way members when building geometries.

Also, for PT routes, I want routes that "cross" the buffer, meaning at least part of the route is within the buffer. So I check if any member way intersects the buffer.

Let me write the final script.

Actually, let me just iterate and finalize. Here's the plan:

1. Download PBF
2. Extract Vienna area with osmium
3. Single pass through Vienna PBF:
   a. Collect all nodes
   b. Collect all ways (store tags and node refs)
   c. Identify Gürtel ways (name ends with 'Gürtel')
   d. Collect PT route relations
4. Build Gürtel geometry, buffer by 500m
5. Filter highway ways in buffer
6. Build PT route geometries from member ways in buffer
7. Write to GeoPackage

Memory for Vienna area should be manageable. Let me estimate: Vienna has about 10-20 million nodes? No, the bbox covers a larger area. But still much less than all of Austria.

Let me write the script now.

The `osmium` Python module is installed, but not the `osmium` CLI tool. I need to use a different approach. I can either:

1. Use `osmium` Python bindings to handle extraction
2. Use `osmnx` (which is installed) to download data directly using the OSMnx API
3. Just read the full Austria PBF with osmium and filter in-memory using the bounding box

Let me rethink the approach. Since the PBF is already downloaded (~1.6 GB), I can use osmium Python to read it and filter by bounding box. But reading 1.6 GB might be memory intensive.

Actually, a better approach: use `osmnx` which has great tools for this kind of thing. But OSMnx normally uses the Overpass API, not local PBF files.

Another option: use `pyosmium` to filter. The Python `osmium` package has `osmium.get_osmosis()` or similar tools? Actually, the Python bindings just have SimpleHandler etc.

Let me try a different approach:
1. Use `osmium` Python to directly read the Austria PBF and build the geometries
2. Or use `osmnx` with the Overpass API to query specifically (smaller data transfer)

Actually, the PBF is already downloaded (1.6 GB). Let me just process it in Python directly with osmium, reading everything but only keeping what's in the Vienna area. This should work fine.

Let me rewrite the script to avoid the osmium CLI.

The script is ready. But I need to think about the PBF size - it's ~1.6 GB, which is large. The first pass reads all nodes and ways from the whole Austria PBF. Memory might be an issue. Let me think about a better approach.

Actually, looking at this more carefully - the full Austria PBF is about 1.6 GB compressed, and OSM objects for all of Austria will be quite large in memory. The nodes alone for Austria could be tens or hundreds of millions. This could exceed memory limits.

Let me think of a more memory-efficient approach:
1. First, extract just the Vienna area using osmium CLI... but we don't have that.
2. Alternatively, we can use osmium via Python's library to do filtering.

Wait, the `osmium` Python package has indexing capabilities. But let me check what tools we have.

Actually, `pyosmium` has `osmium.extract` module? Let me check.

Actually, I can use `osmium.index` and create a filtered PBF by using osmium's handler to write filtered output. Or I can use a different approach: use `osmium` filter to extract a bounding box.

Let me check if we can use `osmium` Python library to extract a region without the CLI tool.

Actually, the Python `osmium` package includes an `Extract` class? Let me check.

Let me try a different approach - use `osmium` with a filter to process only objects within a bounding box around Vienna, and use a two-pass approach more efficiently:

1. First, filter the PBF to just the Vienna bounding box using Python osmium
2. Then process the filtered file

For the filtering step, I can use `osmium.SimpleHandler` and write to a new PBF file using `osmium.io.Writer`. But that's complex.

Actually, a simpler approach: let's use osmnx to download Vienna boundary, then use its graph functions. But osmnx works over the Overpass API, not PBF files.

Let me think about this differently. The script I've written reads everything into memory. For Austria's ~60M nodes and ways, this might be too much.

Better approach: 
1. Read the PBF file in two passes using osmium
2. First pass: Get the Gürtel ways' node IDs and PT relations
3. Second pass: Only read nodes and ways that are relevant

Actually, let me use `osmium.index` with a map file to create a filtered version. Or better yet, use `osmium.SimpleHandler` to write a filtered PBF in the first pass.

Let me use the osmium Python library to write a filtered PBF with just Vienna bbox.

Actually, a much simpler approach: use `osmnx` (which is available) to download the OSM data directly via Overpass API. OSMnx can get the graph for a driving network around a specific area. But the user specifically asked for the Geofabrik PBF.

Let me check if we can use `osmium.extract.Extract` from the Python library.

The script timed out after 10 minutes. Processing the entire Austria PBF (765 MB) with a single-pass osmium handler that collects ALL nodes, ALL ways, and ALL relations is too slow especially when we need to scan the entire file. Let me think of a more efficient approach.

The issue is:
1. The PBF is 765 MB and has millions of objects
2. We're reading all nodes, ways, and relations into memory
3. The script timed out after 10 minutes

Better approach:
1. First, do a quick scan to find Gürtel ways and their nodes + the buffer
2. Then use that buffer to filter only relevant objects

Actually, let me use a two-pass approach with osmium:
- Pass 1: Find Gürtel ways (name ends with Gürtel), their node coordinates, build buffer
- Pass 2: Use the buffer to select intersecting highways and PT route relations

But osmium's handler API processes everything anyway. Let me think...

Actually, the issue might be that the script just needs more time. Let me try a more optimized approach:

1. Use `osmium getid` or `osmium tags-filter` to extract just the Gürtel area more efficiently
2. Or use a faster approach

Actually, let's try a different tack. Since osmium-tool (the CLI) isn't installed, let me install it via conda or pip from a different source. Or maybe use `osmium` Python API more efficiently.

Wait, the Python osmium package (pyosmium) includes a binary tool. Let me check:

```
python -m osmium
```

Or maybe I can use `osmium` command from the Python package location.

The full PBF is 765 MB and the script timed out at 600s. I need to speed this up significantly. The issue is that:

1. We're reading the entire Austria PBF which contains lots of nodes/ways
2. We're doing O(n) geometry operations for every way

Let me think about optimizations:

1. First, we can filter to just the Vienna area using a bounding box filter during the osmium pass
2. We can use osmium's index-based filtering, or just pre-filter with a bounding box

Actually, the osmium python module does allow filtering by area using the `osm.filter` module. Let me check if we can use `osmium.index` and `osmium.filter` to do spatial filtering.

Actually, a better approach: use `osmium` with `apply_file` but with a location handler option that can do on-the-fly filtering. The key bottleneck is likely the geometry creation for every single way.

Let me think about this differently. We need:
1. Ways with name ending in "Gürtel" → we can find them from the full file quickly
2. All highways that intersect their 500m buffer
3. PT route relations that cross the same buffer

Optimization ideas:
- First pass: only collect Gürtel ways + PT relations (fast, since few features)
- Build the buffer
- Second pass: filter all ways using the buffer, only keeping highways or ways referenced by PT relations

This way we don't build geometries for the millions of non-interesting ways.

Actually, even faster: we can use osmnx to get the Gürtel data directly from the PBF, or use a two-pass approach where pass 1 only collects what we need.

Wait, the real slowdown is probably the first pass reading every node and way. With 765 MB and probably tens of millions of nodes, storing all nodes in memory is expensive.

Let me use a smarter approach:
1. First pass: read relations (PT routes) and ways (Gürtel identification, and collect way IDs for PT route members)
2. Build buffer from Gürtel ways
3. Second pass: read ways that are either highways or PT route member ways, and test against buffer
4. Third pass: read nodes referenced by those ways

Actually, osmium allows using a NodeLocationStore or creating an index. Let me use a file-based index.

Better plan:
1. Pass 1: Read relations, collect PT route members (way refs). Read ways, identify Gürtel ways (name ends with 'Gürtel'). Store ALL way geometries temporarily? No...

Let me use a different approach: use `osmium.index.create_map` for node locations, then do two passes.

Actually, the simplest optimization: use a bounding box filter first. Vienna is roughly 16.1-16.7E, 48.05-48.4N. We can use `osmium.filter` to create a filtered PBF, then read that.

Wait, let me check if pyosmium supports `osmium.filter` as a Python function...

Actually, let me just use `osmnx` to get the Gürtel data via Overpass API, or use `osmnx.graph_from_bbox`. But the user asked to use the PBF.

Let me try a different approach: use `osmium.SimpleHandler` but with a node cache that doesn't store every node. Use `osmium.index.create_map` with a flexible key.

Actually the simplest fix: let's do it in phases.

Phase 1: Read only relations and find Gürtel ways (don't store all nodes, just store way IDs and their node lists). This is fast because we don't process nodes.

Phase 2: Read the Gürtel ways' node locations from a second pass to build the buffer.

Phase 3: Read all ways again, test intersection, and for highways + PT member ways that intersect, collect their node refs.

Phase 4: Read nodes for those collected refs.

This avoids storing all nodes in memory.

Let me implement this carefully.

Wait, I'm using `create_map` and `node_location_index[n.id] = (n.location.lon, n.location.lat)` but `create_map` returns an index that stores locations (as osmium.osm.Location objects) not tuples. Let me check the osmium API more carefully.

Actually looking at the osmium docs, `create_map` returns an index map. The typical pattern is:

```python
from osmium.index import create_map
node_locations = create_map("flex_mem")
handler = SomeHandler(node_locations)
handler.apply_file(filepath, locations=True)
```

And then in the handler you can access `node_locations[node_id]` to get a Location.

But writing `node_location_index[n.id] = (n.location.lon, n.location.lat)` won't work properly because the sparse_mem_array expects osmium-specific types.

Let me think about a better approach. Actually the standard osmium approach for this is:

1. Create a location handler that populates a node location store
2. Use `apply_file` with `locations=True` flag which automatically populates the store

Or alternatively, the simplest approach: use the handler that tracks locations itself.

Actually, the fastest approach would be to use OSMnx directly. OSMnx can download the data from a bounding polygon. Let me reconsider.

Wait - OSMnx is installed (2.1.0). We can use osmnx to download the Gürtel area. Actually, osmnx uses the Overpass API which would be slow for large areas. But we have the PBF file locally.

Let me think of a cleaner approach. The issue with the previous script was timeout. Let me use a more memory-efficient but faster approach.

Actually, the key insight: the osmium library can process the entire 765MB PBF quite fast if we use `locations=True` in apply_file. The previous script timed out because the first pass handler tried to store all nodes in a Python dict (nodes = {}), which is slow for millions of nodes.

Let me use the osmium index system properly:

```python
from osmium.index import create_map

class Handler(osmium.SimpleHandler):
    def __init__(self, node_locations):
        super().__init__()
        self.locations = node_locations
    
    def node(self, n):
        if n.location.valid():
            self.locations[n.id] = n.location

node_locations = create_map('flex_mem')
handler = Handler(node_locations)
handler.apply_file(PBF_FILE, locations=True)
```

With `locations=True`, the C++ backend handles node location storage efficiently.

But actually, the issue is simpler: I was storing all nodes as Python tuples in a Python dict, which is super slow for 100M+ nodes. Let me use osmium's built-in location index.

Let me rewrite with this approach.

The `sparse_mem_array` creates a `LocationTable` which can store node IDs -> (lon, lat). Let me check how to use it properly. The `get` method returns a tuple, and `set` sets values. But the `osmium.SimpleHandler` with `locations=True` automatically populates a location table. However, the issue is more fundamental - the node location table is used internally.

Let me simplify the approach. Instead of using complex location tables, I can use osmium's built-in node location handler which is very fast.

Actually, the simplest and fastest approach:
1. Use `osmium.SimpleHandler` with `locations=True` which populates a node location cache automatically during the way/relation handler callbacks.
2. But we need to do multiple passes... 

Actually, the cleanest approach: Use osmium's built-in `NodeLocationsForWays` handler. You attach it to a reader and it stores node locations in memory, then calls the way handler with resolved locations.

But that requires single-pass processing. Let me think differently.

The key issue with my first script was that it was trying to store ALL nodes (millions) in a Python dict, which caused memory and speed issues. But the osmium C++ LocationTable is very efficient.

Let me try a different approach: use osmium's NodeLocationsForWays handler with a single pass to do everything at once.

Actually, let me check how NodeLocationsForWays works in the Python bindings.

In osmium, you can do:
```python
import osmium as o

class Handler(o.SimpleHandler):
    def way(self, w):
        # w.nodes now have valid locations if we use NodeLocationsForWays
        pass

h = Handler()
h.apply_file("file.osm.pbf", locations=True)
```

Wait, the `locations=True` parameter... Let me check.

Actually, in osmium 4.x Python bindings, `SimpleHandler.apply_file` has a `locations` parameter. When set to `True`, it uses `osmium.NodeLocationsForWays` to automatically store node locations and resolve way node references.

But this requires the handler to process both nodes and ways. Let me restructure to do it in one pass with locations enabled.

However, we need multiple passes because:
1. Pass 1: Find Gürtel ways and PT relations
2. Pass 2: Build buffer from Gürtel geometries
3. Pass 3: Collect highways and PT way geometries intersecting buffer

But with locations=True, the handler processes nodes internally and provides resolved node locations in way callbacks.

Let me try this approach:
- Pass 1: Read with locations=True. Store Gürtel way geometries + PT relation info + all highway/way geometries that intersect buffer. But we don't know the buffer yet!

So we need at least 2 passes. Let me optimize:

Pass 1: Read ways only (no locations needed). Find Gürtel way IDs and PT relations. 
Pass 2: Read Gürtel ways with locations=True to build their geometry. Compute buffer.
Pass 3: Read all ways with locations=True. For each highway way or PT member way, test intersection with buffer.

Actually, let's try a simpler approach. The first script was storing ALL nodes in a Python dict which is memory-heavy and slow. Let me use osmium's LocationTable instead, which is C++-backed and much faster.

Actually, rethinking. The problem with the first script was that it timed out after 10 minutes. That script stored all 765M PBF data in Python dicts. Let me measure how many nodes and ways there are.

For Austria, there are roughly:
- ~300-400 million nodes
- ~30-40 million ways
- ~200k relations

Storing 400M node tuples in a Python dict is ~40+ GB of memory and very slow. But storing in osmium's LocationTable is very compact (about 16 bytes per node).

So the right approach:
1. Use osmium's LocationTable (sparse_mem_array) which stores node locations in C++ memory.
2. Read ways and test intersection on-the-fly.

Let me rewrite using the LocationTable properly. The `get` method on the LocationTable returns a `Location` object or raises KeyError.

Actually, I realize the indexing approach is complex. Let me just use `locations=True` flag with the SimpleHandler, which automatically resolves node locations for way nodes. But we need to be careful because with `locations=True`, it processes nodes internally.

Let me try a cleaner approach:

```python
class MyHandler(osmium.SimpleHandler):
    def __init__(self):
        super().__init__()
        # stores results
        
    def way(self, w):
        # w.nodes[i].lat, w.nodes[i].lon are available when locations=True
        pass

handler = MyHandler()
handler.apply_file(PBF_FILE, locations=True)
```

The `locations=True` flag uses osmium's internal `NodeLocationsForWays` handler which stores node locations in an internal buffer. This is very efficient.

Let me write a cleaner script using this approach.

The issue is that reading the full PBF with locations=True is very slow. Let me think of a better approach.

Actually, the problem is that osmium Python's `locations=True` stores nodes at C++ speed internally, but even 800M node PBF is huge. Let me try a different strategy:

1. Use `osmium` to filter the PBF to just the Gürtel and a bounding box around Vienna, much tighter
2. Then work on the smaller filtered PBF

Actually wait - the issue is that the full Austria PBF has ~800M nodes and even reading it is expensive. Let me use a different approach:

1. Use `osmium` in C++ mode to filter ways by tag (highway ways near Gürtel) - but we don't have the CLI tool
2. Use GeoPandas/osmnx approach - osmnx can download from OSM API but that might be too slow for large areas
3. Use a bounding-box filter on the PBF to drastically reduce data

Actually let me try a smarter approach:
- First, find the bbox of the Gürtel ways by scanning just the ways (fast without locations)
- Then filter the PBF to that bbox + buffer using osmium's built-in filtering
- Then work with the much smaller extract

Let me try using osmium's `apply_file` with the `locations=True` flag which stores node locations internally - this is in C++ and should be much faster than the Python handler.

Actually wait, the previous test with locations=True timed out after 5 minutes. That's because it's iterating over all 800M nodes in Python.

Let me try a different approach: use the osmium C++ handler with `locations=True` in the apply_file call. When `locations=True` is passed, osmium internally builds a node location index in C++, and when we access `w.nodes` in the way handler, the locations are already available without needing to separately process nodes. This is the fast approach.

Let me rewrite - use a single pass with locations=True but only filter in the way handler, and skip the node handler callbacks entirely.

This version uses `locations=True` which tells osmium to resolve node coordinates in C++ internally. The key optimization:

1. Pass 1: Just scan for Gürtel way IDs and PT relations (no node resolution needed)
2. Pass 2: With `locations=True`, build Gürtel geometries and compute buffer
3. Pass 3: With `locations=True`, only check ways that are highways or PT member ways, with bbox pre-filter
4. Pass 4: Build PT route MultiLineStrings from collected way geometries

The `locations=True` should be much faster than Python callbacks for every node because osmium handles the node cache in C++. Let me run it.

It timed out at 480 seconds (8 minutes). The process is too slow because Pass 3 is iterating through all ~2.1M+ ways in Austria, checking bbox and building geometries. The buffer bbox covers a very large area (15.16 to 16.40 lon) which is basically all of eastern Austria including Vienna but also lots of other areas.

Let me think about a better approach. The buffer bbox is very large because the Gürtel is in Vienna but the buffer extends quite far. Actually wait - the buffer is only 500m around the Gürtel. Let me check the bbox output: "Buffer bbox: 15.1603, 47.0245, 16.4054, 48.2408". That's 15.16 to 16.40 longitude and 47.02 to 48.24 latitude. That seems like a very large area for a 500m buffer around the Gürtel in Vienna. The Gürtel is in the center of Vienna at roughly 48.20N, 16.37E. So the buffer should be quite small.

Wait - looking at the coordinates: 15.16 to 16.40 longitude. That's about 1.24 degrees. Vienna is at about 16.37E. I'm confused... Let me check. Actually those numbers look wrong. The Gürtel is around 48.20N, 16.37E. A 500m buffer in WGS84 would be about 0.005 degrees in lat and about 0.0045 degrees in lon. But the bbox shows 15.16 to 16.40 which is huge!

Oh wait, I think the issue is the buffer is being reprojected incorrectly. Let me look at the code more carefully.

The code does:
1. Build Gürtel geometry in WGS84
2. Convert to TARGET_CRS (EPSG:31287)
3. Buffer 500m in projected CRS
4. Convert back to WGS84

But wait - the Gürtel ways might not all be near each other. "Gürtel" in Vienna is just one circular road, but maybe there are other ways named ending in "Gürtel" across Austria (Gürtelstrasse, etc.)? Let me check - the filter is `name.endswith('Gürtel')`. In Vienna the main road is called "Wiener Gürtel" (or just "Gürtel") but there might be other Gürtel streets across Austria. That's fine for the filter, but they should be geographically clustered near Vienna.

Actually wait, 415 Gürtel ways - there could be Gürtel ways all over Austria. But when we compute `unary_union(guertel_lines)` from all 415 Gürtel ways across Austria, the union might produce a huge geometry spanning the whole country, and then buffering that by 500m gives a massive area.

Hmm, but the user said "Vienna's Gürtel, defined as the highway ways whose name ends in 'Gürtel'". So we need to be smarter - we need to find the Gürtel IN Vienna. The problem is that "Gürtel" ending names exist across Austria, not just in Vienna.

Let me fix this:
1. First find ways ending in "Gürtel" 
2. Filter to those near Vienna (say, within a bbox of Vienna)
3. Then compute the buffer

Let me also significantly speed up Pass 3 by being smarter about it.

Actually, I think the biggest issue is that we're reading ~2M ways in Python, which is slow. Let me use a much more targeted approach. Let me use osmium to extract a small area around Vienna first, then work with that.

Actually, maybe a better approach: extract just a small area around Vienna using osmium's extract functionality, then process that. But osmium CLI isn't available.

Let me try a different approach:
1. Identify the Vienna Gürtel ways by looking at which Gürtel ways are clustered around Vienna
2. Create a small buffer around the Gürtel
3. Use osmium's filtered apply to only process ways in that area

Or even simpler: since we know Vienna Gürtel is around 48.2N, 16.37E, let me use a tight bbox filter first, then compute the buffer.

Let me think about the performance issue. The buffer bbox is too large because there are Gürtel-named ways spread across Austria. For instance, "Grazer Gürtel" in Graz (south of Vienna) and similar streets in other cities. Let me check what the bbox looks like.

Wait, the previous run showed:
```
Buffer bbox: 15.1603, 47.0245, 16.4054, 48.2408
```

That's actually not too bad - min_lon=15.16, min_lat=47.02, max_lon=16.41, max_lat=48.24. This covers from roughly Wiener Neustadt to Vienna, and from ~Styria to north of Vienna. That's still a huge area.

But the issue is that the previous run was checking ALL highways that had any node in that large bbox. Let me think...

Actually, the bbox is large because some Gürtel-named ways in other parts of Austria are being included in the Vienna Gürtel group. The algorithm detected them as being within 0.5 degrees of Vienna center... Let me recalculate: 0.5 degrees is ~55km. 

Vienna center: 16.37, 48.20
The bbox min coords: 15.16, 47.02

15.16 is at the longitude of... St. Pölten area (west of Vienna). And 47.02 is way south, near Graz level. 

So some Gürtel-named ways in the Graz area (~15.4°E, 47.1°N) are being included. The distance from Vienna center to Graz center is about 1.1° in latitude and 0.97° in longitude, so about 1.46 degrees total, which is > 0.5 but with a generous threshold...

Wait, 0.5 degrees radius is ~55km. Graz is ~150km from Vienna. So Graz Gürtel shouldn't be included with a 0.5 threshold.

But the bbox results show 15.16, 47.02 as min. Let me check: Straßwalchen? That's 13.26°E. No. 

Actually 15.16°E is near... maybe there are Gürtel named roads in Wiener Neustadt (16.23°E, 47.81°N), Mödling (16.29°E, 48.08°N), or even further.

Let me just check: the bbox min is 15.16°E, 47.02°N. That's near... Hartberg? Or maybe there's a Brucker Gürtel or something in Bruck an der Mur (15.27°E, 47.41°N).

The issue is some Gürtel way far south got included. Let me fix this by using a tighter spatial filter. I'll compute the centroid of all Vienna Gürtel ways and use a smaller radius.

Actually, let me just check the output to see which ways are Gürtel-named in Austria and then filter more carefully.

Actually, let me try a different approach entirely. Instead of iterating over the whole Austria PBF for pass 2/3, let me:
1. Figure out the tightest possible bounding box for the Vienna Gürtel + 500m buffer
2. Use osmium extract-like filtering to process only that region

But osmium doesn't have a pure-Python extract filter. However, I can do a bbox pre-filter in the handler more efficiently.

The real problem is that `locations=True` builds a full node index of ALL nodes in the PBF, which for 765MB PBF means ~hundreds of millions of nodes. That's the bottleneck - the C++ node location cache needs to process every node in Austria.

Better approach: Instead of using `locations=True` on the full PBF, I should:
1. First extract a tight region around Vienna (using osmium extract, but it failed earlier because the tool wasn't installed)
2. Process just that smaller file

Let me install osmium-tool or use a workaround.

Actually, let me try a completely different strategy using osmium's file-level extraction in Python, or use a simple bounding box approach.

Wait - there's `osmium.extract` but it's a CLI tool. Let me check if there's a Python equivalent.

Actually, I think the fastest approach would be:
1. Find Vienna Gürtel way IDs in Pass 1
2. Determine their bounding box 
3. Use osmium with `locations=True` + the ability to skip nodes outside the bbox? No, osmium doesn't support that.

Alternative: use the `idx` parameter of `apply_file` to use a flexible memory index, which might be faster?

Or better: use a two-file approach:
- Use `osmium` with the `locations=True` but pass a tighter `idx` type
- Or pre-filter with a simple bounding box read

Let me look at osmium's Python API more carefully. There might be a way to use `osmium.io.Reader` with filters, or use pyosmium's `apply_file` with the `locations=True` parameter efficiently.

Actually, the real issue might be memory/performance of the `sparse_mem_array` index on a 765MB PBF. Let me try `flex_mem` (the default) which osmium uses internally.

Actually, looking at the output more carefully - Pass 3 was checking 2.1M+ ways and continuing to count up. That's because the buffer bbox is very large (encompassing most of eastern Austria). The bbox is ~1.2° x 1.2° which at 111km/degree is ~133km x 133km ~ 17,700 km². That's huge - it's basically all of eastern Austria.

The issue is some Gürtel-named ways far from Vienna being included. Let me narrow down to actual Vienna Gürtel ways by looking at their names - the famous Gürtel in Vienna is usually "Gürtel" or has a prefix like "Währinger Gürtel", "Gürtel" etc.

Actually, the name endings are things like:
- "Gürtel" (the main one)
- "Währinger Gürtel"
- "Gürtel" etc.

But there are also streets named "Grazer Gürtel" in Graz, etc.

Let me just be smarter about it. I'll compute the center of mass of all Gürtel ways, filter to those within 0.25 degrees (~28km), which should give us only the Vienna ones.

Wait, but the previous run's output showed bbox: 15.1603, 47.0245, 16.4054, 48.2408

That's a HUGE area for a 500m buffer around the Vienna Gürtel. Let me think...

The Gürtel in Vienna is a ring road around the inner city. Its approximate bounds are:
- West: ~16.32°E (near Westbahnhof)
- East: ~16.38°E 
- North: ~48.23°N
- South: ~48.18°N

So a 500m buffer (about 0.0045°) would extend this to roughly:
- minx: ~16.30, miny: ~48.17, maxx: ~16.39, maxy: ~48.24

But we got 15.16-16.41 x 47.02-48.24. This means some Gürtel way at lon=15.16, lat=47.02 is being included.

15.16°E, 47.02°N is... that's in Styria, near... let me think. Graz is at 15.44°E, 47.07°N. So this is southwest of Graz.

The issue might be that the clustering threshold of 0.5 degrees is too generous. Let me be more precise.

Actually, I know - "Grazer Gürtel" is a street in Graz (47.07°N, 15.44°E). Distance from Vienna (48.20°N, 16.37°E):
dlat = 48.20 - 47.07 = 1.13°
dlon = 16.37 - 15.44 = 0.93°
dist = sqrt(1.13² + 0.93²) = 1.46° ≈ 162 km

But the threshold was 0.5° (~55km). So Grazer Gürtel should NOT be included. Unless there's also a Gürtel in... hmm.

Let me look: 15.16°E, 47.02°N - that's near Oberwart (16.20°E, 47.29°N)? No, 15.16 is west.

Actually 15.16°E, 47.02°N - that's near the border of Styria and Lower Austria, maybe near... Mürzzuschlag (15.68°E, 47.59°N)? No.

Let me just check what Gürtel-named streets exist in that area. But the simpler fix is to use a MUCH tighter threshold. The Vienna Gürtel is very specific. Let me use 0.1 degrees (~11km).

Wait - the previous run also said "Gürtel way IDs: 415" which seems like a lot for the Vienna Gürtel specifically. The actual Vienna Gürtel probably has ~50-100 way segments. So there are ~300+ Gürtel-named ways elsewhere in Austria.

Let me fix the approach: Instead of using a distance threshold, let me just filter using the known approximate bbox of Vienna's Gürtel, or better yet, use the OSM admin boundary. 

Actually, the simplest fix: The Vienna Gürtel is centered around 16.35°E, 48.20°N. Let me use a threshold of 0.1° (roughly 11 km) which should be plenty for the Gürtel + its surrounding area.

Let me also think about a completely different approach - using OSMnx to fetch just the Vienna area. OSMnx can query OSM directly. But we need all the data from a PBF...

Actually, let me try yet another approach. Instead of processing the entire Austria PBF, I can:
1. Find the Vienna Gürtel bbox
2. Use `osmium`'s node cache only on a geoparsed subset
3. Or use `osmium.io.Reader` with filters

Actually, I think the simplest fix is to just tighten the spatial filter and it will work. The bbox will be small enough that the `locations=True` pass will be fast.

Let me also consider: can I pre-filter the PBF file using osmium's extract functionality in Python?

Actually, there's a much simpler approach. Let me check if `osmium` supports file-level operations like extracting by bounding box using its C++ API through Python. Let me look at what's available:

`osmium.io.File`, `osmium.io.Reader`, `osmium.io.Writer`, `osmium.filter.GeoFilter` - these exist!

Let me check if I can use osmium's GeoFilter or if there's a way to do a bbox-based extraction in pure Python.

Actually wait - I can use `osmium.io.Reader` with a `osmium.filter.GeoFilter` to pre-filter nodes. Or I can use `osmium.filter.EntityFilter` with some geometry filtering.

Let me check what's available:

```python
from osmium.filter import GeoFilter, EntityFilter
```

I think the simplest approach is to:
1. Pass 1: Find Vienna Gürtel ways and compute tight bbox
2. Use `osmium` with locations=True but with a much tighter bbox (0.1° search radius for Vienna)
3. This will make the node cache small and the way iteration fast

Let me update the code.

Actually, I just realized - the issue was that in the test run, the buffer bbox was HUGE (15.16 to 16.41). That means the Vienna cluster included Gürtel ways from far away. The distance threshold of 0.5° was too generous. Let me use 0.1° which is ~11km, more than enough to capture the actual Vienna Gürtel.

Let me also consider that the Gürtel in Vienna has many segments. The entire ring is maybe 10-15 km in circumference. So a 0.1° radius from center should definitely capture all of it.

Let me rewrite with this fix.

The key issue is that with `locations=True`, osmium needs to process ALL nodes in the PBF to build the location index, which is slow for a 765MB file. With `locations=True`, EVERY node is loaded into the C++ index. This means each pass takes a while but it should be much faster now that we have a tight bbox.

Let me check how many total nodes Austria has. The previous run was killed at 8 min. Let me see what's happening.

Actually, the timeout was 480 seconds (8 min) and it was still in Pass 3 checking ways. The issue is that even with locations=True, processing ALL nodes (to build location index) takes a lot of time.

But wait - the first run had a huge buffer bbox because it was using all Gürtel ways in Austria (Gürtel in Linz, Graz, etc.). Now with the tight Vienna filter, the buffer bbox should be much smaller.

Let me try running it again. The issue was the bbox was too large before. Now it should be much tighter.

Actually let me think about this more carefully. The bbox in the last timeout was:
```
Buffer bbox: 15.1603, 47.0245, 16.4054, 48.2408
```

This is because there are Gürtel-named ways all over Austria (e.g., "Gürtel" in Linz, etc.). The buffer was computed from ALL Gürtel ways in Austria, not just Vienna's.

My new fix uses tighter 0.15° clustering around Vienna (48.20N, 16.37E) which should give a bbox of roughly 16.2 to 16.5 and 48.05 to 48.35. That's much smaller.

Let me run it.

The output looks great! Everything is working correctly. Let me verify a few more things:

1. The file has two layers: `highways` (LineString) and `pt_routes` (MultiLineString) ✓
2. CRS is EPSG:31287 (MGI / Austria Lambert) ✓
3. Highway attributes: name, highway, maxspeed, lanes, surface, oneway ✓
4. PT route attributes: ref, name, operator, route ✓
5. The Gürtel ways are present (222 with Gürtel names) ✓
6. PT routes include bus (219), train (86), tram (49), subway (10), light_rail (2) ✓
7. No trolleybus routes in the data (that's fine, there may not be any in the area) ✓

The output is perfect. Let me provide a summary.