fio-l3-vienna-geofabrik-highways

Rerun

import osmium

class HighwayFilter(osmium.SimpleReader):
    def __init__(self, handler):
        super(HighwayFilter, self).__init__(handler)

class WaysHandler(osmium.SimpleHandler):
    def __init__(self):
        super(WaysHandler, self).__init__()
        self.ways = []

    def way(self, w):
        if 'highway' in w.tags:
            # We need the name to filter by 'Gürtel'
            # But we also need to know if it's a highway.
            # We store it for later processing.
            self.ways.append(w)

# Since we need to filter by name 'Gürtel' and intersection with a buffer,
# we should first find the geometries of the Gürtel highways.

import osmium
import geopandas as gpd
import pandas as pd
from shapely.geometry import LineString, MultiLineString
from pyproj import Transformer
import requests

# 1. Download Geofabrik Austria
# 2. Parse PBF to find Gürtel highways and their geometries
# 3. Create 500m buffer around Gürtel
# 4. Filter highways in buffer
# 5. Filter PT routes in buffer
# 6. Export to GeoPackage

def download_data(url, dest):
    print(f"Downloading {url}...")
    response = requests.get(url, stream=True)
    with open(dest, 'wb') as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)
    print("Download complete.")

# For the sake of this environment, I'll assume the file might be large.
# I will use osmium to extract ways and relations.
# I'll first identify the Gürtel ways.

import osmium
import os
import requests
import geopandas as gpd
import pandas as pd
from shapely.geometry import LineString, MultiLineString
from shapely.ops import linemerge
import pyproj

# Configuration
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"
OUTPUT_FILE = "vienna_network.gpkg"
BUFFER_DIST = 500 # meters

# EPSG for Austria: EPSG:31257 (MGI / Austria GK East) or EPSG:3035 (ETRS89 / LAEA Europe)
# The prompt asks for Austria's standard projected coordinate system. 
# EPSG:31257 is commonly used for Austria.
TARGET_CRS = "EPSG:31257"

def download_data():
    if not os.path.exists(PBF_FILE):
        print(f"Downloading {PBF_URL}...")
        r = requests.get(PBF_URL, stream=True)
        with open(PBF_FILE, 'wb') as f:
            for chunk in r.iter_content(chunk_size=1024*1024):
                if chunk:
                    f.write(chunk)
        print("Download finished.")
    else:
        print("PBF file already exists.")

class GürtelFinder(osmium.SimpleHandler):
    def __init__(self):
        super().__init__()
        self.guertel_ways = []

    def way(self, w):
        name = w.tags.get('name', '')
        if name.endswith('Gürtel') and 'highway' in w.tags:
            # Store the way for geometry extraction
            self.guertel_ways.append(w)

# Note: osmium SimpleHandler doesn't provide geometry directly in way() without a location handler.
# We need a LocationContext.

class GeometryHandler(osmium.SimpleHandler):
    def __init__(self):
        super().__init__()
        self.guertel_geometries = []
        self.all_highways = []
        self.pt_relations = []

    def way(self, w):
        # Check if it's a highway
        if 'highway' in w.tags:
            # We will collect all highways to filter later
            # But we need geometries.
            pass

    def relation(self, r):
        # Handle PT relations
        pass

# Actually, osmium's way to get geometry is by using a LocationContext.
# Let's rewrite the approach.

def extract_data():
    import osmium
    
    class Collector(osmium.SimpleHandler):
        def __init__(self):
            super().__init__()
            self.guertel_ways = []
            self.highways = []
            self.pt_relations = []

        def way(self, w):
            if 'highway' in w.tags:
                self.highways.append(w)
                if w.tags.get('name', '').endswith('Gürtel'):
                    self.guertel_ways.append(w)

        def relation(self, r):
            # Filter for PT routes
            pt_types = {'bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail'}
            # We check tags of the relation if available, or look for route type in members
            # Standard OSM: relation has 'route' tag
            if r.tags.get('route') in pt_types:
                self.pt_relations.append(r)

    # However, collecting ALL highways from Austria in memory might crash.
    # Austria PBF is ~700MB. The ways list will be huge.
    # Better approach: 
    # 1. Find Gürtel way IDs and geometries first.
    # 2. Use those to create a buffer (in EPSG:31257).
    # 3. Run a second pass to collect highways and PT relations that intersect the buffer.

    print("Pass 1: Finding Gürtel geometries...")
    guertel_geoms = []

    class GürtelGeometryFinder(osmium.SimpleHandler):
        def __init__(self):
            super().__init__()
            self.geoms = []

        def way(self, w):
            if w.tags.get('name', '').endswith('Gürtel') and 'highway' in w.tags:
                try:
                    # Get geometry
                    # osmium way objects don't have geometry unless you use a location handler
                    # but we can use osmium.geom.WKBFactory or similar if we have a location provider
                    pass
                except:
                    pass

    # Let's use a more robust approach using osmium.SimpleHandler with a location handler.
    pass

if __name__ == "__main__":
    pass

import os
import osmium
import osmium.geom
import geopandas as gpd
import pandas as pd
from shapely.geometry import LineString, MultiLineString
import requests

# Configuration
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"
OUTPUT_FILE = "vienna_network.gpkg"
BUFFER_DIST = 500  # meters
TARGET_CRS = "EPSG:31257"

def download_data():
    if not os.path.exists(PBF_FILE):
        print(f"Downloading {PBF_URL}...")
        r = requests.get(PBF_URL, stream=True)
        with open(PBF_FILE, 'wb') as f:
            for chunk in r.iter_content(chunk_size=1024*1024):
                if chunk:
                    f.write(chunk)
        print("Download finished.")
    else:
        print("PBF file already exists.")

class GuertelFinder(osmium.SimpleHandler):
    def __init__(self, geom_factory):
        super().__init__()
        self.geom_factory = geom_factory
        self.guertel_geoms = []

    def way(self, w):
        if w.tags.get('name', '').endswith('Gürtel') and 'highway' in w.tags:
            try:
                geom = self.geom_factory.create_linestring(w)
                if geom:
                    self.guertel_geoms.append(geom)
            except Exception:
                pass

class NetworkCollector(osmium.SimpleHandler):
    def __init__(self, geom_factory, buffer_geom):
        super().__init__()
        self.geom_factory = geom_factory
        self.buffer_geom = buffer_geom
        self.highways = []
        self.pt_routes = []

    def way(self, w):
        if 'highway' in w.tags:
            try:
                geom = self.geom_factory.create_linestring(w)
                if geom and geom.intersects(self.buffer_geom):
                    # Extract attributes
                    attrs = {
                        'name': w.tags.get('name', ''),
                        'highway': w.tags.get('highway', ''),
                        'maxspeed': w.tags.get('maxspeed', ''),
                        'lanes': w.tags.get('lanes', ''),
                        'surface': w.tags.get('surface', ''),
                        'oneway': w.tags.get('oneway', '')
                    }
                    self.highways.append({'geometry': geom, 'properties': attrs})
            except Exception:
                pass

    def relation(self, r):
        pt_types = {'bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail'}
        if r.tags.get('route') in pt_types:
            try:
                # For relations, we need to reconstruct geometry from members
                # This is tricky with osmium without a full location handler/cache.
                # However, we can use the 'osm_id' to get geometry if we have a way to resolve it,
                # or we can use the fact that osmium handles relations if we use a proper handler.
                # Let's use a different approach for PT: we'll collect relation members and build geometry.
                pass
            except Exception:
                pass

# Re-evaluating: The relation geometry is complex. 
# Let's use a simpler approach:
# 1. Find Gürtel geometries.
# 2. Buffer them.
# 3. Use a tool or a more efficient way to get intersecting geometries.

# Given the constraints and the complexity of OSM relations in osmium:
# I'll use a simpler two-pass approach.
# Pass 1: Find Gürtel way geometries.
# Pass 2: Collect all highways and PT routes.

# Since I cannot easily reconstruct MultiLineStrings for PT relations from a single pass 
# in osmium without a heavy cache, I will try to use the 'osmium.SimpleHandler' 
# with 'apply' and 'location' if possible, or just collect all way geometries 
# that are part of a relation.

# Actually, a very common way to handle this is to:
# 1. Find Gürtel IDs.
# 2. Find all ways/relations that intersect the buffer.

# Let's try a simpler script that uses osmium to extract ONLY the necessary data.

def solve():
    download_data()
    
    geom_factory = osmium.geom.WKBFactory()
    
    print("Pass 1: Finding Gürtel geometries...")
    guertel_finder = GuertelFinder(osmium.geom.FeatureGeometryFactory())
    # Use a simpler reader for Pass 1
    reader = osmium.SimpleReader(PBF_FILE)
    guertel_finder.apply(reader)
    
    if not guertel_finder.guertel_geoms:
        print("No Gürtel highways found!")
        return

    # Create buffer
    # We need to convert to projected CRS to buffer in meters
    # Let's use a GeoSeries for this
    from shapely.geometry import shape
    import pyproj
    from shapely.ops import transform

    # Create a dummy GeoDataFrame to handle CRS and buffering
    gdf_guertel = gpd.GeoDataFrame({'geometry': guertel_finder.guertel_geoms}, crs="EPSG:4326")
    gdf_guertel_proj = gdf_guertel.to_crs(TARGET_CRS)
    buffer_geom = gdf_guertel_proj.buffer(BUFFER_DIST).union_all()
    
    # Convert buffer back to 4326 for the second pass (to work with osmium/osm coords)
    buffer_geom_4326 = transform(
        pyproj.Transformer.from_crs(TARGET_CRS, "EPSG:4326", always_xy=True).transform,
        buffer_geom
    )
    # Actually, it's easier to keep everything in 4326 for the pass and then project.
    # But buffer is in meters.
    
    # Let's refine: buffer_geom is in EPSG:4326.
    # We'll use it to filter in Pass 2.
    
    print("Pass 2: Collecting highways and PT routes...")
    # We will collect ALL ways and ALL relations, then filter.
    # To handle PT routes as single features, we need to collect all their member ways.
    
    # To avoid memory issues, we'll use a more targeted approach.
    # We'll collect all ways that intersect the buffer.
    # And all relations that intersect the buffer.
    
    # Note: osmium's 'apply' is efficient.
    
    class NetworkCollectorV2(osmium.SimpleHandler):
        def __init__(self, buffer_geom_4326):
            super().__init__()
            self.buffer_geom_4326 = buffer_geom_4326
            self.highways = []
            self.pt_relations_data = {} # relation_id -> {attrs, member_ways}

        def way(self, w):
            # For every way, check if it's a highway and intersects buffer
            if 'highway' in w.tags:
                try:
                    # Use a simpler way to get geometry
                    # This might be slow, but let's see
                    pass
                except:
                    pass
    
    # Actually, let's use the fact that we can iterate over ways and relations.
    # Given the time and memory limits, I'll write a robust script.
    pass

if __name__ == "__main__":
    solve()

import os
import osmium
import osmium.geom
import geopandas as gpd
import pandas as pd
from shapely.geometry import LineString, MultiLineString
import pyproj
from shapely.ops import transform
import requests

# Configuration
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"
OUTPUT_FILE = "vienna_network.gpkg"
BUFFER_DIST = 500  # meters
TARGET_CRS = "EPSG:31257"

def download_data():
    if not os.path.exists(PBF_FILE):
        print(f"Downloading {PBF_URL}...")
        r = requests.get(PBF_URL, stream=True)
        with open(PBF_FILE, 'wb') as f:
            for chunk in r.iter_content(chunk_size=1024*1024):
                if chunk:
                    f.write(chunk)
        print("Download finished.")
    else:
        print("PBF file already exists.")

class GuertelFinder(osmium.SimpleHandler):
    def __init__(self, geom_factory):
        super().__init__()
        self.geom_factory = geom_factory
        self.guertel_geoms = []

    def way(self, w):
        if w.tags.get('name', '').endswith('Gürtel') and 'highway' in w.tags:
            try:
                geom = self.geom_factory.create_linestring(w)
                if geom:
                    self.guertel_geoms.append(geom)
            except Exception:
                pass

class NetworkCollector(osmium.SimpleHandler):
    def __init__(self, buffer_geom_4326, geom_factory):
        super().__init__()
        self.buffer_geom_4326 = buffer_geom_4326
        self.geom_factory = geom_factory
        self.highways = []
        self.pt_relations = []
        # For PT routes, we need to collect all ways in a relation.
        # We'll store them in a dict keyed by relation ID.
        self.relation_ways = {} # relation_id -> list of ways

    def way(self, w):
        # Check if it's a highway
        if 'highway' in w.tags:
            try:
                geom = self.geom_factory.create_linestring(w)
                if geom and geom.intersects(self.buffer_geom_4326):
                    attrs = {
                        'name': w.tags.get('name', ''),
                        'highway': w.tags.get('highway', ''),
                        'maxspeed': w.tags.get('maxspeed', ''),
                        'lanes': w.tags.get('lanes', ''),
                        'surface': w.tags.get('surface', ''),
                        'oneway': w.tags.get('oneway', '')
                    }
                    self.highways.append({'geometry': geom, 'properties': attrs})
            except Exception:
                pass
        
        # Check if this way is part of a PT relation we are tracking
        # We'll need to handle this in the relation() or by checking relation members.
        # A better way: 
        # 1. Find all PT relations that intersect the buffer.
        # 2. Collect all ways belonging to those relations.

    def relation(self, r):
        # This is called for every relation.
        # If it's a PT route, we check if it intersects.
        # But 'r' doesn't have geometry directly. 
        # We need to check its members.
        pass

# Since the above is getting complex for a single script, 
# I will use a simpler, more direct approach.
# 1. Use osmium to get all ways and relations.
# 2. Store only those that intersect the buffer.

def solve():
    download_data()
    
    geom_factory = osmium.geom.WKBFactory()
    
    print("Pass 1: Finding Gürtel geometries...")
    guertel_finder = GuertelFinder(osmium.geom.FeatureGeometryFactory())
    guertel_finder.apply(osmium.SimpleReader(PBF_FILE))
    
    if not guertel_finder.guertel_geoms:
        print("No Gürtel highways found!")
        return

    # Create buffer in 4326
    gdf_guertel = gpd.GeoDataFrame({'geometry': guertel_finder.guertel_geoms}, crs="EPSG:4326")
    gdf_guertel_proj = gdf_guertel.to_crs(TARGET_CRS)
    buffer_geom_proj = gdf_guertel_proj.buffer(BUFFER_DIST).union_all()
    
    # Transform buffer back to 4326
    transformer = pyproj.Transformer.from_crs(TARGET_CRS, "EPSG:4326", always_xy=True)
    buffer_geom_4326 = transform(transformer.transform, buffer_geom_proj)

    print("Pass 2: Collecting all highways and PT routes...")
    
    # To handle PT relations, we need to:
    # 1. Identify PT relations.
    # 2. For each PT relation, get its member ways.
    # 3. Merge the ways to form the route geometry.
    
    # We will do this in one pass by:
    # - Collecting all highways that intersect buffer.
    # - Collecting all PT relations and their members.
    # - Then filtering PT relations by buffer.

    class FullCollector(osmium.SimpleHandler):
        def __init__(self, buffer_geom_4326, geom_factory):
            super().__init__()
            self.buffer_geom_4326 = buffer_geom_4326
            self.geom_factory = geom_factory
            self.highways = []
            self.pt_relations_info = {} # relation_id -> {attrs, member_ways_geoms}

        def way(self, w):
            # 1. Highways
            if 'highway' in w.tags:
                try:
                    geom = self.geom_factory.create_linestring(w)
                    if geom and geom.intersects(self.buffer_geom_4326):
                        self.highways.append({
                            'geometry': geom,
                            'name': w.tags.get('name', ''),
                            'highway': w.tags.get('highway', ''),
                            'maxspeed': w.tags.get('maxspeed', ''),
                            'lanes': w.tags.get('lanes', ''),
                            'surface': w.tags.get('surface', ''),
                            'oneway': w.tags.get('oneway', '')
                        })
                except Exception:
                    pass
            
            # 2. PT Relation members (we need to know which relation they belong to)
            # This is hard because 'way' doesn't know its parent relation.
            # We'll use a different approach for PT routes.

    # Let's try another way for PT routes: 
    # Use osmium to find all PT relations first.
    
    class PTRelationFinder(osmium.SimpleHandler):
        def __init__(self):
            super().__init__()
            self.pt_relation_ids = []
            self.pt_relation_attrs = {}

        def relation(self, r):
            pt_types = {'bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail'}
            if r.tags.get('route') in pt_types:
                self.pt_relation_ids.append(r.id)
                self.pt_relation_attrs[r.id] = {
                    'ref': r.tags.get('ref', ''),
                    'name': r.tags.get('name', ''),
                    'operator': r.tags.get('operator', ''),
                    'route': r.tags.get('route', '')
                }

    print("Finding PT relations...")
    pt_finder = PTRelationFinder()
    pt_finder.apply(osmium.SimpleReader(PBF_FILE))
    
    # Now we need to get the geometry for these PT relations.
    # This is the hardest part with osmium. 
    # A relation's geometry is the union of its member ways.
    # We'll do a third pass to collect all ways that are members of these PT relations.

    class PTWayCollector(osmium.SimpleHandler):
        def __init__(self, pt_relation_ids, geom_factory):
            super().__init__()
            self.pt_relation_ids = set(pt_relation_ids)
            self.geom_factory = geom_factory
            self.relation_to_ways = {rid: [] for rid in pt_relation_ids}

        def relation(self, r):
            # We don't need to do anything in relation()
            pass

        def way(self, w):
            # This is still not quite right. osmium's way() doesn't give parent relations.
            # BUT, osmium.SimpleHandler.relation(r) gives us r.members.
            pass

    # Let's use the 'osm_id' approach.
    # Pass 1: Get Gürtel buffer.
    # Pass 2: Get all highways in buffer.
    # Pass 3: Get all PT relations in buffer.
    # To do Pass 3:
    #   - Find all PT relations.
    #   - For each relation, get its member way IDs.
    #   - Collect geometries for all those way IDs.
    #   - Check if the resulting MultiLineString intersects the buffer.

    # I'll use a simpler strategy for the PT routes to ensure it works.
    # I'll collect all ways that are part of *any* PT relation.
    # Then I'll group them by relation.

    class WayRelationCollector(osmium.SimpleHandler):
        def __init__(self, geom_factory):
            super().__init__()
            self.geom_factory = geom_factory
            self.ways_in_relations = [] # List of (relation_id, way_geom, way_attrs)
            self.pt_relations_attrs = {}

        def relation(self, r):
            pt_types = {'bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail'}
            if r.tags.get('route') in pt_types:
                self.pt_relations_attrs[r.id] = {
                    'ref': r.tags.get('ref', ''),
                    'name': r.tags.get('name', ''),
                    'operator': r.tags.get('operator', ''),
                    'route': r.tags.get('route', '')
                }
                # We'll use this to collect members in the next step? 
                # No, osmium provides members in the relation object.
                
        def way(self, w):
            # We can't easily link ways to relations in 'way()'
            pass
            
    # Actually, the easiest way is to use 'apply' on the reader and 
    # in the 'relation' handler, iterate over members.
    
    class PTCollector(osmium.SimpleHandler):
        def __init__(self, buffer_geom_4326, geom_factory):
            super().__init__()
            self.buffer_geom_4326 = buffer_geom_4326
            self.geom_factory = geom_factory
            self.highways = []
            self.pt_routes = []

        def way(self, w):
            if 'highway' in w.tags:
                try:
                    geom = self.geom_factory.create_linestring(w)
                    if geom and geom.intersects(self.buffer_geom_4326):
                        self.highways.append({
                            'geometry': geom,
                            'name': w.tags.get('name', ''),
                            'highway': w.tags.get('highway', ''),
                            'maxspeed': w.tags.get('maxspeed', ''),
                            'lanes': w.tags.get('lanes', ''),
                            'surface': w.tags.get('surface', ''),
                            'oneway': w.tags.get('oneway', '')
                        })
                except Exception:
                    pass

        def relation(self, r):
            pt_types = {'bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail'}
            if r.tags.get('route') in pt_types:
                # Collect all geometries of member ways
                route_geoms = []
                for member in r.members:
                    if member.type == 'way':
                        try:
                            # This is the tricky part: we need the geometry of the member way.
                            # osmium's relation.members only gives IDs and types.
                            # To get the geometry, we'd need to have the ways in a cache.
                            pass
                        except Exception:
                            pass
                # This approach is hitting a wall because of how osmium works.
                # I will use a different approach.

    # FINAL ATTEMPT AT STRATEGY:
    # 1. Find Gürtel geometries and create buffer.
    # 2. Use osmium to collect ALL ways and ALL relations.
    # 3. For highways: if 'highway' in tags and intersects buffer, save.
    # 4. For relations: if 'route' in pt_types, collect all its member way geometries.
    #    To do this efficiently:
    #    - First, collect all ways that are part of ANY PT relation into a dictionary {way_id: geometry}.
    #    - Second, collect all PT relations and their member way IDs.
    #    - Third, reconstruct PT route geometries and filter by buffer.

    # Let's do this. It's two passes.

    print("Pass 1: Finding Gürtel geometries...")
    guertel_finder = GuertelFinder(osmium.geom.FeatureGeometryFactory())
    guertel_finder.apply(osmium.SimpleReader(PBF_FILE))
    
    if not guertel_finder.guertel_geoms:
        print("No Gürtel highways found!")
        return

    gdf_guertel = gpd.GeoDataFrame({'geometry': guertel_finder.guertel_geoms}, crs="EPSG:4326")
    gdf_guertel_proj = gdf_guertel.to_crs(TARGET_CRS)
    buffer_geom_proj = gdf_guertel_proj.buffer(BUFFER_DIST).union_all()
    transformer = pyproj.Transformer.from_crs(TARGET_CRS, "EPSG:4326", always_xy=True)
    buffer_geom_4326 = transform(transformer.transform, buffer_geom_proj)

    print("Pass 2: Collecting highways and PT way-to-relation mappings...")
    
    class Pass2Collector(osmium.SimpleHandler):
        def __init__(self, buffer_geom_4326, geom_factory):
            super().__init__()
            self.buffer_geom_4326 = buffer_geom_4326
            self.geom_factory = geom_factory
            self.highways = []
            self.pt_way_geoms = {} # way_id -> geometry
            self.pt_relations = [] # List of (relation_id, {attrs}, [member_way_ids])

        def way(self, w):
            # 1. Highways
            if 'highway' in w.tags:
                try:
                    geom = self.geom_factory.create_linestring(w)
                    if geom and geom.intersects(self.buffer_geom_4326):
                        self.highways.append({
                            'geometry': geom,
                            'name': w.tags.get('name', ''),
                            'highway': w.tags.get('highway', ''),
                            'maxspeed': w.tags.get('maxspeed', ''),
                            'lanes': w.tags.get('lanes', ''),
                            'surface': w.tags.get('surface', ''),
                            'oneway': w.tags.get('oneway', '')
                        })
                except Exception:
                    pass
            
            # We'll also store ALL way geometries that *might* be part of a PT relation
            # But that's too many. Let's just collect all ways and their geometries 
            # but only for the ways that are actually members of PT relations.
            # Wait, I'll just collect all way geometries for ALL ways in Pass 2.
            # If the PBF is large, this might be too much memory.
            # Let's only collect ways that are members of PT relations.
            # But I don't know which ones are in PT relations yet.
            
            # OK, let's collect ALL ways in Pass 2 and their geometries.
            # To save memory, I'll only store geometry if it's a highway OR 
            # if it's a way (we don't know if it's PT yet).
            # Actually, let's just store all way geometries in a dict.
            try:
                # We'll use a more memory-efficient way: only store ways in a dict 
                # if we see them in a relation. But we need to see the relation first.
                pass
            except:
                pass

    # This is getting very circular. Let's use the most reliable method for osmium:
    # 1. Find Gürtel buffer.
    # 2. Find all PT relations (IDs and attributes).
    # 3. Find all highways that intersect the buffer.
    # 4. Find all ways that are members of the PT relations found in step 2.
    # 5. Reconstruct PT routes and check if they intersect the buffer.

    print("Pass 1: Finding Gürtel buffer...")
    guertel_finder = GuertelFinder(osmium.geom.FeatureGeometryFactory())
    guertel_finder.apply(osmium.SimpleReader(PBF_FILE))
    gdf_guertel = gpd.GeoDataFrame({'geometry': guertel_finder.guertel_geoms}, crs="EPSG:4326")
    buffer_geom_proj = gdf_guertel.to_crs(TARGET_CRS).buffer(BUFFER_DIST).union_all()
    buffer_geom_4326 = transform(pyproj.Transformer.from_crs(TARGET_CRS, "EPSG:4326", always_xy=True).transform, buffer_geom_proj)

    print("Pass 2: Finding PT relations...")
    class PTFinder(osmium.SimpleHandler):
        def __init__(self):
            super().__init__()
            self.pt_relations = {} # id -> {attrs, member_way_ids}
        def relation(self, r):
            pt_types = {'bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail'}
            if r.tags.get('route') in pt_types:
                self.pt_relations[r.id] = {
                    'attrs': {
                        'ref': r.tags.get('ref', ''),
                        'name': r.tags.get('name', ''),
                        'operator': r.tags.get('operator', ''),
                        'route': r.tags.get('route', '')
                    },
                    'member_way_ids': [m.ref for m in r.members if m.type == 'way']
                }
    pt_finder = PTFinder()
    pt_finder.apply(osmium.SimpleReader(PBF_FILE))
    
    print(f"Found {len(pt_finder.pt_relations)} PT relations.")

    print("Pass 3: Collecting highways and PT way geometries...")
    class Collector(osmium.SimpleHandler):
        def __init__(self, buffer_4326, pt_way_ids, geom_factory):
            super().__init__()
            self.buffer_4326 = buffer_4326
            self.pt_way_ids = set(pt_way_ids)
            self.geom_factory = geom_factory
            self.highways = []
            self.pt_way_geoms = {} # way_id -> geometry

        def way(self, w):
            # Highway
            if 'highway' in w.tags:
                try:
                    geom = self.geom_factory.create_linestring(w)
                    if geom and geom.intersects(self.buffer_4326):
                        self.highways.append({
                            'geometry': geom,
                            'name': w.tags.get('name', ''),
                            'highway': w.tags.get('highway', ''),
                            'maxspeed': w.tags.get('maxspeed', ''),
                            'lanes': w.tags.get('lanes', ''),
                            'surface': w.tags.get('surface', ''),
                            'oneway': w.tags.get('oneway', '')
                        })
                except Exception:
                    pass
            
            # PT way
            if w.id in self.pt_way_ids:
                try:
                    geom = self.geom_factory.create_linestring(w)
                    if geom:
                        self.pt_way_geoms[w.id] = geom
                except Exception:
                    pass

    target_way_ids = set()
    for rel in pt_finder.pt_relations.values():
        target_way_ids.update(rel['member_way_ids'])
    
    collector = Collector(buffer_geom_4326, target_way_ids, osmium.geom.FeatureGeometryFactory())
    collector.apply(osmium.SimpleReader(PBF_FILE))

    # Reconstruct PT routes
    pt_routes_list = []
    for rid, info in pt_finder.pt_relations.items():
        route_geoms = []
        for wid in info['member_way_ids']:
            if wid in collector.pt_way_geoms:
                route_geoms.append(collector.pt_way_geoms[wid])
        
        if route_geoms:
            # Merge the ways into a single MultiLineString
            try:
                merged_geom = MultiLineString(route_geoms)
                # Check if the route intersects the buffer
                if merged_geom.intersects(buffer_geom_4326):
                    pt_routes_list.append({
                        'geometry': merged_geom,
                        **info['attrs']
                    })
            except Exception:
                pass

    # Create GeoDataFrames
    print("Creating GeoDataFrames...")
    highways_gdf = gpd.GeoDataFrame(collector.highways, crs="EPSG:4326")
    if not highways_gdf.empty:
        highways_gdf = highways_gdf.set_geometry('geometry')
    else:
        # Create empty GeoDataFrame with correct columns
        highways_gdf = gpd.GeoDataFrame(columns=['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway'], geometry='geometry', crs="EPSG:4326")

    pt_routes_gdf = gpd.GeoDataFrame(pt_routes_list, crs="EPSG:4326")
    if not pt_routes_gdf.empty:
        pt_routes_gdf = pt_routes_gdf.set_geometry('geometry')
    else:
        pt_routes_gdf = gpd.GeoDataFrame(columns=['ref', 'name', 'operator', 'route'], geometry='geometry', crs="EPSG:4326")

    # Project both to TARGET_CRS
    highways_gdf = highways_gdf.to_crs(TARGET_CRS)
    pt_routes_gdf = pt_routes_gdf.to_crs(TARGET_CRS)

    # Write to GeoPackage
    print(f"Writing to {OUTPUT_FILE}...")
    highways_gdf.to_file(OUTPUT_FILE, layer='highways', driver="GPKG")
    pt_routes_gdf.to_file(OUTPUT_FILE, layer='pt_routes', driver="GPKG")
    print("Done!")

if __name__ == "__main__":
    solve()