fio-l3-vienna-geofabrik-highways

Rerun

import os
import requests
import geopandas as gpd
from pyrosm import OSM, get_data
from shapely.geometry import LineString, MultiLineString
from shapely.ops import linemerge
import pandas as pd

# 1. Download the Geofabrik Austria PBF
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"

if not os.path.exists(PBF_FILE):
    print(f"Downloading {PBF_URL}...")
    response = requests.get(PBF_URL, stream=True)
    with open(PBF_FILE, "wb") as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)
print("Download complete.")

# 2. Initialize OSM parser
osm = OSM(PBF_FILE)

# 3. Identify the Gürtel highways (name ends with 'Gürtel')
# We'll first get all highways to filter them. 
# Note: pyrosm's get_network is efficient for roads.
print("Extracting highways...")
network = osm.get_network(network_type="driving")

# Filter highways where name ends with 'Gürtel'
# pyrosm returns a GeoDataFrame.
gurtel_highways = network[network['name'].str.endswith('Gürtel', na=False)].copy()

if gurtel_highways.empty:
    # Try case-insensitive or just check if any highway has 'Gürtel' in it
    print("No highways found ending with 'Gürtel'. Searching for 'Gürtel' in name...")
    gurtel_highways = network[network['name'].str.contains('Gürtel', na=False)].copy()

if gurtel_highways.empty:
    raise ValueError("Could not find any highway with 'Gürtel' in its name.")

# 4. Create a 500m buffer around the Gürtel highways
# Use Austria's standard projected coordinate system (EPSG:31256 - MGI / Austria GK East or EPSG:31257 - MGI / Austria GK West)
# Or more commonly EPSG:3035 (ETRS89-extended / LAEA Europe) or EPSG:25832 (ETRS89 / UTM zone 32N)
# The user asked for "Austria's standard projected coordinate system". 
# EPSG:31256 is Gauss-Krüger. Let's use EPSG:31256 for Austria.
# Actually, let's check the center of Vienna to be safe or just use EPSG:31256.
# For a small area like Vienna, EPSG:31256/31257 is fine. 
# Let's use EPSG:31256.

gurtel_highways_proj = gurtel_highways.to_crs(epsg=31256)
buffer_500m = gurtel_highways_proj.buffer(500).unary_union

print(f"Buffer created. Area of buffer: {buffer_500m.area} m^2")

# 5. Get all highways intersecting the buffer
# We need the full attributes: name, highway, maxspeed, lanes, surface, oneway
# pyrosm's get_network usually includes many of these. 
# Let's re-extract or filter.
all_highways = network.to_crs(epsg=31256)
intersecting_highways = all_highways[all_highways.intersects(buffer_500m)].copy()

# Select required columns for highways
highway_cols = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway']
# Filter to only those that exist in the dataframe
available_highway_cols = [c for c in highway_cols if c in intersecting_highways.columns]
intersecting_highways = intersecting_highways[available_highway_cols + ['geometry']]

# Ensure they are LineStrings
intersecting_highways['geometry'] = intersecting_highways['geometry'].apply(lambda x: LineString(x.coords) if x.geom_type != 'LineString' else x)

# 6. Get Public Transport routes
# pyrosm's get_public_transport returns routes.
print("Extracting public transport routes...")
pt_routes = osm.get_public_transport()

if pt_routes is not None and not pt_routes.empty:
    # Filter for bus, tram, subway, train, trolleybus, light_rail
    # In OSM, transit types are often in 'route' (bus, tram, subway, etc.)
    valid_transit_types = ['bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail']
    pt_routes_filtered = pt_routes[pt_routes['route'].isin(valid_transit_types)].copy()
    
    # Project PT routes
    pt_routes_proj = pt_routes_filtered.to_crs(epsg=31256)
    
    # Intersect with buffer
    intersecting_pt = pt_routes_proj[pt_routes_proj.intersects(buffer_500m)].copy()
    
    # Select required columns for PT: ref, name, operator, route
    pt_cols = ['ref', 'name', 'operator', 'route']
    available_pt_cols = [c for c in pt_cols if c in intersecting_pt.columns]
    intersecting_pt = intersecting_pt[available_pt_cols + ['geometry']]
    
    # Ensure they are MultiLineString (as requested)
    # A route might be a LineString or MultiLineString. We want MultiLineString.
    def to_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom

    intersecting_pt['geometry'] = intersecting_pt['geometry'].apply(to_multilinestring)
else:
    print("No public transport routes found.")
    intersecting_pt = gpd.GeoDataFrame(columns=['ref', 'name', 'operator', 'route', 'geometry'], crs=31256)

# 7. Save to GeoPackage
print("Saving to vienna_network.gpkg...")
with pd.io.formats.formatters.BaseFormatter.default_args: # dummy to allow context manager if needed, but we use geopandas
    pass

intersecting_highways.to_file("vienna_network.gpkg", layer="highways", driver="GPKG")
intersecting_pt.to_file("vienna_network.gpkg", layer="pt_routes", driver="GPKG")

print("Done. File 'vienna_network.gpkg' created.")

import os
import requests
import geopandas as gpd
from pyrosm import OSM
from shapely.geometry import LineString, MultiLineString
import pandas as pd

# 1. Download the Geofabrik Austria PBF
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"

if not os.path.exists(PBF_FILE):
    print(f"Downloading {PBF_URL}...")
    response = requests.get(PBF_URL, stream=True)
    with open(PBF_FILE, "wb") as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)
print("Download complete.")

# 2. Initialize OSM parser
osm = OSM(PBF_FILE)

# 3. Identify the Gürtel highways (name ends with 'Gürtel')
print("Extracting highways...")
network = osm.get_network(network_type="driving")

# Filter highways where name ends with 'Gürtel'
# We need to handle cases where 'name' might be NaN or not a string
gurtel_highways = network[network['name'].astype(str).str.endswith('Gürtel')].copy()

if gurtel_highways.empty:
    print("No highways found ending with 'Gürtel'. Searching for 'Gürtel' in name...")
    gurtel_highways = network[network['name'].astype(str).str.contains('Gürtel')].copy()

if gurtel_highways.empty:
    raise ValueError("Could not find any highway with 'Gürtel' in its name.")

# 4. Create a 500m buffer around the Gürtel highways
# EPSG:31256 (Austria GK East) is suitable for Vienna.
gurtel_highways_proj = gurtel_highways.to_crs(epsg=31256)
buffer_500m = gurtel_highways_proj.buffer(500).unary_union

print(f"Buffer created. Area of buffer: {buffer_500m.area} m^2")

# 5. Get all highways intersecting the buffer
all_highways = network.to_crs(epsg=31256)
intersecting_highways = all_highways[all_highways.intersects(buffer_500m)].copy()

# Select required columns for highways
highway_cols = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway']
available_highway_cols = [c for c in highway_cols if c in intersecting_highways.columns]
intersecting_highways = intersecting_highways[available_highway_cols + ['geometry']]

# Ensure they are LineStrings
def ensure_linestring(geom):
    if geom.geom_type == 'MultiLineString':
        # If it's a MultiLineString, we'll try to merge it or just take the first part?
        # The prompt asks for LineString. Let's try to merge if possible.
        try:
            return LineString(geom.geoms[0].coords) # fallback
        except:
            return geom
    return geom

intersecting_highways['geometry'] = intersecting_highways['geometry'].apply(lambda x: LineString(x.coords) if x.geom_type != 'LineString' else x)

# 6. Get Public Transport routes
print("Extracting public transport routes...")
# pyrosm.get_public_transport returns relations. 
# It might not return everything directly in a single go, but let's try.
pt_routes = osm.get_public_transport()

if pt_routes is not None and not pt_routes.empty:
    # Filter for bus, tram, subway, train, trolleybus, light_rail
    valid_transit_types = ['bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail']
    # pyrosm's pt_routes might use 'route' column
    pt_routes_filtered = pt_routes[pt_routes['route'].isin(valid_transit_types)].copy()
    
    # Project PT routes
    pt_routes_proj = pt_routes_filtered.to_crs(epsg=31256)
    
    # Intersect with buffer
    intersecting_pt = pt_routes_proj[pt_routes_proj.intersects(buffer_500m)].copy()
    
    # Select required columns for PT: ref, name, operator, route
    pt_cols = ['ref', 'name', 'operator', 'route']
    available_pt_cols = [c for c in pt_cols if c in intersecting_pt.columns]
    intersecting_pt = intersecting_pt[available_pt_cols + ['geometry']]
    
    # Ensure they are MultiLineString (as requested)
    def to_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom

    intersecting_pt['geometry'] = intersecting_pt['geometry'].apply(to_multilinestring)
else:
    print("No public transport routes found.")
    intersecting_pt = gpd.GeoDataFrame(columns=['ref', 'name', 'operator', 'route', 'geometry'], crs=31256)

# 7. Save to GeoPackage
print("Saving to vienna_network.gpkg...")
intersecting_highways.to_file("vienna_network.gpkg", layer="highways", driver="GPKG")
intersecting_pt.to_file("vienna_network.gpkg", layer="pt_routes", driver="GPKG")

print("Done. File 'vienna_network.gpkg' created.")

import os
import requests
import geopandas as gpd
from pyrosm import OSM
from shapely.geometry import LineString, MultiLineString
import pandas as pd

# 1. Download the Geofabrik Austria PBF
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"

if not os.path.exists(PBF_FILE):
    print(f"Downloading {PBF_URL}...")
    response = requests.get(PBF_URL, stream=True)
    with open(PBF_FILE, "wb") as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)
print("Download complete.")

# 2. Initialize OSM parser
osm = OSM(PBF_FILE)

# 3. Identify the Gürtel highways (name ends with 'Gürtel')
print("Extracting highways...")
network = osm.get_network(network_type="driving")

# Filter highways where name ends with 'Gürtel'
# We need to handle cases where 'name' might be NaN or not a string
gurtel_highways = network[network['name'].astype(str).str.endswith('Gürtel')].copy()

if gurtel_highways.empty:
    print("No highways found ending with 'Gürtel'. Searching for 'Gürtel' in name...")
    gurtel_highways = network[network['name'].astype(str).str.contains('Gürtel')].copy()

if gurtel_highways.empty:
    raise ValueError("Could not find any highway with 'Gürtel' in its name.")

# 4. Create a 500m buffer around the Gürtel highways
# EPSG:31256 (Austria GK East) is suitable for Vienna.
gurtel_highways_proj = gurtel_highways.to_crs(epsg=31256)
buffer_500m = gurtel_highways_proj.buffer(500).unary_union

print(f"Buffer created. Area of buffer: {buffer_500m.area} m^2")

# 5. Get all highways intersecting the buffer
all_highways = network.to_crs(epsg=31256)
intersecting_highways = all_highways[all_highways.intersects(buffer_500m)].copy()

# Select required columns for highways
highway_cols = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway']
available_highway_cols = [c for c in highway_cols if c in intersecting_highways.columns]
intersecting_highways = intersecting_highways[available_highway_cols + ['geometry']]

# Ensure they are LineStrings
def ensure_linestring(geom):
    if geom.geom_type == 'MultiLineString':
        # For the purpose of this task, if it's a MultiLineString, we'll merge the parts into a single LineString if possible
        # Or just use the first part to comply with LineString requirement.
        try:
            # Trying to merge parts
            from shapely.ops import linemerge
            merged = linemerge(geom)
            if merged.geom_type == 'LineString':
                return merged
            elif merged.geom_type == 'MultiLineString':
                return merged.geoms[0]
        except:
            pass
        return LineString(geom.geoms[0].coords) if len(geom.geoms) > 0 else geom
    return geom

intersecting_highways['geometry'] = intersecting_highways['geometry'].apply(ensure_linestring)

# 6. Get Public Transport routes
print("Extracting public transport routes...")
pt_routes = osm.get_public_transport()

if pt_routes is not None and not pt_routes.empty:
    # Filter for bus, tram, subway, train, trolleybus, light_rail
    valid_transit_types = ['bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail']
    # pyrosm's pt_routes might use 'route' column
    pt_routes_filtered = pt_routes[pt_routes['route'].isin(valid_transit_types)].copy()
    
    # Project PT routes
    pt_routes_proj = pt_routes_filtered.to_crs(epsg=31256)
    
    # Intersect with buffer
    intersecting_pt = pt_routes_proj[pt_routes_proj.intersects(buffer_500m)].copy()
    
    # Select required columns for PT: ref, name, operator, route
    pt_cols = ['ref', 'name', 'operator', 'route']
    available_pt_cols = [c for c in pt_cols if c in intersecting_pt.columns]
    intersecting_pt = intersecting_pt[available_pt_cols + ['geometry']]
    
    # Ensure they are MultiLineString (as requested)
    def to_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom

    intersecting_pt['geometry'] = intersecting_pt['geometry'].apply(to_multilinestring)
else:
    print("No public transport routes found.")
    intersecting_pt = gpd.GeoDataFrame(columns=['ref', 'name', 'operator', 'route', 'geometry'], crs=31256)

# 7. Save to GeoPackage
print("Saving to vienna_network.gpkg...")
intersecting_highways.to_file("vienna_network.gpkg", layer="highways", driver="GPKG")
intersecting_pt.to_file("vienna_network.gpkg", layer="pt_routes", driver="GPKG")

print("Done. File 'vienna_network.gpkg' created.")

import os
import requests
import geopandas as gpd
from pyrosm import OSM
from shapely.geometry import LineString, MultiLineString
import pandas as pd

# 1. Download the Geofabrik Austria PBF
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"

if not os.path.exists(PBF_FILE):
    print(f"Downloading {PBF_URL}...")
    response = requests.get(PBF_URL, stream=True)
    with open(PBF_FILE, "wb") as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)
print("Download complete.")

# 2. Initialize OSM parser
osm = OSM(PBF_FILE)

# 3. Identify the Gürtel highways (name ends with 'Gürtel')
print("Extracting highways...")
network = osm.get_network(network_type="driving")

# Filter highways where name ends with 'Gürtel'
gurtel_highways = network[network['name'].astype(str).str.endswith('Gürtel')].copy()

if gurtel_highways.empty:
    print("No highways found ending with 'Gürtel'. Searching for 'Gürtel' in name...")
    gurtel_highways = network[network['name'].astype(str).str.contains('Gürtel')].copy()

if gurtel_highways.empty:
    raise ValueError("Could not find any highway with 'Gürtel' in its name.")

# 4. Create a 500m buffer around the Gürtel highways
# EPSG:31256 (Austria GK East) is suitable for Vienna.
gurtel_highways_proj = gurtel_highways.to_crs(epsg=31256)
buffer_500m = gurtel_highways_proj.buffer(500).unary_union

print(f"Buffer created. Area of buffer: {buffer_500m.area} m^2")

# 5. Get all highways intersecting the buffer
all_highways = network.to_crs(epsg=31256)
intersecting_highways = all_highways[all_highways.intersects(buffer_500m)].copy()

# Select required columns for highways
highway_cols = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway']
available_highway_cols = [c for c in highway_cols if c in intersecting_highways.columns]
intersecting_highways = intersecting_highways[available_highway_cols + ['geometry']]

# Ensure they are LineStrings
def ensure_linestring(geom):
    if geom.geom_type == 'MultiLineString':
        try:
            from shapely.ops import linemerge
            merged = linemerge(geom)
            if merged.geom_type == 'LineString':
                return merged
            elif merged.geom_type == 'MultiLineString':
                return merged.geoms[0]
        except:
            pass
        return LineString(geom.geoms[0].coords) if len(geom.geoms) > 0 else geom
    return geom

intersecting_highways['geometry'] = intersecting_highways['geometry'].apply(ensure_linestring)

# 6. Get Public Transport routes
print("Extracting public transport routes...")
pt_routes = osm.get_public_transport()

if pt_routes is not None and not pt_routes.empty:
    # Filter for bus, tram, subway, train, trolleybus, light_rail
    valid_transit_types = ['bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail']
    pt_routes_filtered = pt_routes[pt_routes['route'].isin(valid_transit_types)].copy()
    
    # Project PT routes
    pt_routes_proj = pt_routes_filtered.to_crs(epsg=31256)
    
    # Intersect with buffer
    intersecting_pt = pt_routes_proj[pt_routes_proj.intersects(buffer_500m)].copy()
    
    # Select required columns for PT: ref, name, operator, route
    pt_cols = ['ref', 'name', 'operator', 'route']
    available_pt_cols = [c for c in pt_cols if c in intersecting_pt.columns]
    intersecting_pt = intersecting_pt[available_pt_cols + ['geometry']]
    
    # Ensure they are MultiLineString (as requested)
    def to_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom

    intersecting_pt['geometry'] = intersecting_pt['geometry'].apply(to_multilinestring)
else:
    print("No public transport routes found.")
    intersecting_pt = gpd.GeoDataFrame(columns=['ref', 'name', 'operator', 'route', 'geometry'], crs=31256)

# 7. Save to GeoPackage
print("Saving to vienna_network.gpkg...")
intersecting_highways.to_file("vienna_network.gpkg", layer="highways", driver="GPKG")
intersecting_pt.to_file("vienna_network.gpkg", layer="pt_routes", driver="GPKG")

print("Done. File 'vienna_network.gpkg' created.")

import os
import requests
import geopandas as gpd
from pyrosm import OSM
from shapely.geometry import LineString, MultiLineString
import pandas as pd

# 1. Download the Geofabrik Austria PBF
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"

if not os.path.exists(PBF_FILE):
    print(f"Downloading {PBF_URL}...")
    response = requests.get(PBF_URL, stream=True)
    with open(PBF_FILE, "wb") as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)
print("Download complete.")

# 2. Initialize OSM parser
osm = OSM(PBF_FILE)

# 3. Identify the Gürtel highways (name ends with 'Gürtel')
print("Extracting highways...")
network = osm.get_network(network_type="driving")

# Filter highways where name ends with 'Gürtel'
gurtel_highways = network[network['name'].astype(str).str.endswith('Gürtel')].copy()

if gurtel_highways.empty:
    print("No highways found ending with 'Gürtel'. Searching for 'Gürtel' in name...")
    gurtel_highways = network[network['name'].astype(str).str.contains('Gürtel')].copy()

if gurtel_highways.empty:
    raise ValueError("Could not find any highway with 'Gürtel' in its name.")

# 4. Create a 500m buffer around the Gürtel highways
# EPSG:31256 (Austria GK East) is suitable for Vienna.
gurtel_highways_proj = gurtel_highways.to_crs(epsg=31256)
buffer_500m = gurtel_highways_proj.buffer(500).unary_union

print(f"Buffer created. Area of buffer: {buffer_500m.area} m^2")

# 5. Get all highways intersecting the buffer
all_highways = network.to_crs(epsg=31256)
intersecting_highways = all_highways[all_highways.intersects(buffer_500m)].copy()

# Select required columns for highways
highway_cols = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway']
available_highway_cols = [c for c in highway_cols if c in intersecting_highways.columns]
intersecting_highways = intersecting_highways[available_highway_cols + ['geometry']]

# Ensure they are LineStrings
def ensure_linestring(geom):
    if geom.geom_type == 'MultiLineString':
        try:
            from shapely.ops import linemerge
            merged = linemerge(geom)
            if merged.geom_type == 'LineString':
                return merged
            elif merged.geom_type == 'MultiLineString':
                return merged.geoms[0]
        except:
            pass
        return LineString(geom.geoms[0].coords) if len(geom.geoms) > 0 else geom
    return geom

intersecting_highways['geometry'] = intersecting_highways['geometry'].apply(ensure_linestring)

# 6. Get Public Transport routes
print("Extracting public transport routes...")
pt_routes = osm.get_public_transport()

if pt_routes is not None and not pt_routes.empty:
    # Filter for bus, tram, subway, train, trolleybus, light_rail
    valid_transit_types = ['bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail']
    pt_routes_filtered = pt_routes[pt_routes['route'].isin(valid_transit_types)].copy()
    
    # Project PT routes
    pt_routes_proj = pt_routes_filtered.to_crs(epsg=31256)
    
    # Intersect with buffer
    intersecting_pt = pt_routes_proj[pt_routes_proj.intersects(buffer_500m)].copy()
    
    # Select required columns for PT: ref, name, operator, route
    pt_cols = ['ref', 'name', 'operator', 'route']
    available_pt_cols = [c for c in pt_cols if c in intersecting_pt.columns]
    intersecting_pt = intersecting_pt[available_pt_cols + ['geometry']]
    
    # Ensure they are MultiLineString (as requested)
    def to_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom

    intersecting_pt['geometry'] = intersecting_pt['geometry'].apply(to_multilinestring)
else:
    print("No public transport routes found.")
    intersecting_pt = gpd.GeoDataFrame(columns=['ref', 'name', 'operator', 'route', 'geometry'], crs=31256)

# 7. Save to GeoPackage
print("Saving to vienna_network.gpkg...")
intersecting_highways.to_file("vienna_network.gpkg", layer="highways", driver="GPKG")
intersecting_pt.to_file("vienna_network.gpkg", layer="pt_routes", driver="GPKG")

print("Done. File 'vienna_network.gpkg' created.")

import os
import requests
import geopandas as gpd
from pyrosm import OSM
from shapely.geometry import LineString, MultiLineString
import pandas as pd

# 1. Download the Geofabrik Austria PBF
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"

if not os.path.exists(PBF_FILE):
    print(f"Downloading {PBF_URL}...")
    response = requests.get(PBF_URL, stream=True)
    with open(PBF_FILE, "wb") as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)
print("Download complete.")

# 2. Initialize OSM parser
osm = OSM(PBF_FILE)

# 3. Identify the Gürtel highways (name ends with 'Gürtel')
print("Extracting highways...")
network = osm.get_network(network_type="driving")

# Filter highways where name ends with 'Gürtel'
gurtel_highways = network[network['name'].astype(str).str.endswith('Gürtel')].copy()

if gurtel_highways.empty:
    print("No highways found ending with 'Gürtel'. Searching for 'Gürtel' in name...")
    gurtel_highways = network[network['name'].astype(str).str.contains('Gürtel')].copy()

if gurtel_highways.empty:
    raise ValueError("Could not find any highway with 'Gürtel' in its name.")

# 4. Create a 500m buffer around the Gürtel highways
# EPSG:31256 (Austria GK East) is suitable for Vienna.
gurtel_highways_proj = gurtel_highways.to_crs(epsg=31256)
buffer_500m = gurtel_highways_proj.buffer(500).unary_union

print(f"Buffer created. Area of buffer: {buffer_500m.area} m^2")

# 5. Get all highways intersecting the buffer
all_highways = network.to_crs(epsg=31256)
intersecting_highways = all_highways[all_highways.intersects(buffer_500m)].copy()

# Select required columns for highways
highway_cols = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway']
available_highway_cols = [c for c in highway_cols if c in intersecting_highways.columns]
intersecting_highways = intersecting_highways[available_highway_cols + ['geometry']]

# Ensure they are LineStrings
def ensure_linestring(geom):
    if geom.geom_type == 'MultiLineString':
        try:
            from shapely.ops import linemerge
            merged = linemerge(geom)
            if merged.geom_type == 'LineString':
                return merged
            elif merged.geom_type == 'MultiLineString':
                return merged.geoms[0]
        except:
            pass
        return LineString(geom.geoms[0].coords) if len(geom.geoms) > 0 else geom
    return geom

intersecting_highways['geometry'] = intersecting_highways['geometry'].apply(ensure_linestring)

# 6. Get Public Transport routes
print("Extracting public transport routes...")
pt_routes = osm.get_public_transport()

if pt_routes is not None and not pt_routes.empty:
    # Filter for bus, tram, subway, train, trolleybus, light_rail
    valid_transit_types = ['bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail']
    pt_routes_filtered = pt_routes[pt_routes['route'].isin(valid_transit_types)].copy()
    
    # Project PT routes
    pt_routes_proj = pt_routes_filtered.to_crs(epsg=31256)
    
    # Intersect with buffer
    intersecting_pt = pt_routes_proj[pt_routes_proj.intersects(buffer_500m)].copy()
    
    # Select required columns for PT: ref, name, operator, route
    pt_cols = ['ref', 'name', 'operator', 'route']
    available_pt_cols = [c for c in pt_cols if c in intersecting_pt.columns]
    intersecting_pt = intersecting_pt[available_pt_cols + ['geometry']]
    
    # Ensure they are MultiLineString (as requested)
    def to_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom

    intersecting_pt['geometry'] = intersecting_pt['geometry'].apply(to_multilinestring)
else:
    print("No public transport routes found.")
    intersecting_pt = gpd.GeoDataFrame(columns=['ref', 'name', 'operator', 'route', 'geometry'], crs=31256)

# 7. Save to GeoPackage
print("Saving to vienna_network.gpkg...")
intersecting_highways.to_file("vienna_network.gpkg", layer="highways", driver="GPKG")
intersecting_pt.to_file("vienna_network.gpkg", layer="pt_routes", driver="GPKG")

print("Done. File 'vienna_network.gpkg' created.")

import os
import requests
import geopandas as gpd
from pyrosm import OSM
from shapely.geometry import LineString, MultiLineString
import pandas as pd

# 1. Download the Geofabrik Austria PBF
PBF_URL = "https://download.geofabrik.de/europe/austria-latest.osm.pbf"
PBF_FILE = "austria-latest.osm.pbf"

if not os.path.exists(PBF_FILE):
    print(f"Downloading {PBF_URL}...")
    response = requests.get(PBF_URL, stream=True)
    with open(PBF_FILE, "wb") as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)
print("Download complete.")

# 2. Initialize OSM parser
osm = OSM(PBF_FILE)

# 3. Identify the Gürtel highways (name ends with 'Gürtel')
print("Extracting highways...")
network = osm.get_network(network_type="driving")

# Filter highways where name ends with 'Gürtel'
gurtel_highways = network[network['name'].astype(str).str.endswith('Gürtel')].copy()

if gurtel_highways.empty:
    print("No highways found ending with 'Gürtel'. Searching for 'Gürtel' in name...")
    gurtel_highways = network[network['name'].astype(str).str.contains('Gürtel')].copy()

if gurtel_highways.empty:
    raise ValueError("Could not find any highway with 'Gürtel' in its name.")

# 4. Create a 500m buffer around the Gürtel highways
# EPSG:31256 (Austria GK East) is suitable for Vienna.
gurtel_highways_proj = gurtel_highways.to_crs(epsg=31256)
buffer_500m = gurtel_highways_proj.buffer(500).unary_union

print(f"Buffer created. Area of buffer: {buffer_500m.area} m^2")

# 5. Get all highways intersecting the buffer
all_highways = network.to_crs(epsg=31256)
intersecting_highways = all_highways[all_highways.intersects(buffer_500m)].copy()

# Select required columns for highways
highway_cols = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway']
available_highway_cols = [c for c in highway_cols if c in intersecting_highways.columns]
intersecting_highways = intersecting_highways[available_highway_cols + ['geometry']]

# Ensure they are LineStrings
def ensure_linestring(geom):
    if geom.geom_type == 'MultiLineString':
        try:
            from shapely.ops import linemerge
            merged = linemerge(geom)
            if merged.geom_type == 'LineString':
                return merged
            elif merged.geom_type == 'MultiLineString':
                return merged.geoms[0]
        except:
            pass
        return LineString(geom.geoms[0].coords) if len(geom.geoms) > 0 else geom
    return geom

intersecting_highways['geometry'] = intersecting_highways['geometry'].apply(ensure_linestring)

# 6. Get Public Transport routes
print("Extracting public transport routes...")
pt_routes = osm.get_public_transport()

if pt_routes is not None and not pt_routes.empty:
    # Filter for bus, tram, subway, train, trolleybus, light_rail
    valid_transit_types = ['bus', 'tram', 'subway', 'train', 'trolleybus', 'light_rail']
    pt_routes_filtered = pt_routes[pt_routes['route'].isin(valid_transit_types)].copy()
    
    # Project PT routes
    pt_routes_proj = pt_routes_filtered.to_crs(epsg=31256)
    
    # Intersect with buffer
    intersecting_pt = pt_routes_proj[pt_routes_proj.intersects(buffer_500m)].copy()
    
    # Select required columns for PT: ref, name, operator, route
    pt_cols = ['ref', 'name', 'operator', 'route']
    available_pt_cols = [c for c in pt_cols if c in intersecting_pt.columns]
    intersecting_pt = intersecting_pt[available_pt_cols + ['geometry']]
    
    # Ensure they are MultiLineString (as requested)
    def to_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom

    intersecting_pt['geometry'] = intersecting_pt['geometry'].apply(to_multilinestring)
else:
    print("No public transport routes found.")
    intersecting_pt = gpd.GeoDataFrame(columns=['ref', 'name', 'operator', 'route', 'geometry'], crs=31256)

# 7. Save to GeoPackage
print("Saving to vienna_network.gpkg...")
intersecting_highways.to_file("vienna_network.gpkg", layer="highways", driver="GPKG")
intersecting_pt.to_file("vienna_network.gpkg", layer="pt_routes", driver="GPKG")

print("Done. File 'vienna_network.gpkg' created.")