fio-l3-vienna-geofabrik-highways

Rerun

import osmnx as ox
import geopandas as gpd
from shapely.geometry import MultiLineString, LineString

def main():
    # 1. Define the Gürtel highways to find the buffer zone
    # We search for highways with 'Gürtel' in the name in Vienna.
    # Since we don't have the PBF locally (curl failed), 
    # let's use osmnx to query Overpass API directly for the Gürtel.
    
    print("Fetching Gürtel highways...")
    # Querying for ways where name contains 'Gürtel'
    # Using a more robust query: name ~ 'Gürtel'
    tags = {"highway": True, "name": "~Gürtel"}
    # We'll limit the area to Vienna to avoid global searches
    # but osmnx's geocode can find Vienna.
    vienna_location = "Vienna, Austria"
    
    # Instead of relying on the PBF (which I can't download easily without curl/wget), 
    # I will use OSMnx to query the Overpass API which is much more efficient for this task.
    
    # Get the Gürtel ways
    # We'll use ox.features_from_address or similar, but custom queries are better.
    # Custom query for highways ending in 'Gürtel' or containing 'Gürtel'
    # The prompt says "name ends in 'Gürtel'".
    
    custom_query = (
        '(way["highway"]["name"~"Gürtel$"](area["name"="Wien"]););'
    )
    # The above might be too strict with $ if it's "Gürtelring" or something.
    # Let's try a slightly broader search if it fails, but stick to the prompt.
    # Actually, 'Gürtel' is often part of the name. Let's use name ~ 'Gürtel' 
    # then filter for those ending in 'Gürtel' or just matching the intent.
    
    # Let's try to get the Gürtel features first.
    try:
        gurtel_ways = ox.features_from_address(vienna_location, tags={"highway": True, "name": "~Gürtel"})
    except Exception as e:
        print(f"Error fetching Gürtel: {e}")
        # Fallback: just get all highways in Vienna and filter manually if needed, 
        # but that's too much data. Let's try a better query.
        return

    # Filter ways whose name ends in 'Gürtel' (case insensitive/trimmed)
    gurtel_ways['name_clean'] = gurtel_ways['name'].astype(str).str.strip()
    gurtel_ways = gurtel_ways[gurtel_ways['name_clean'].str.endswith('Gürtel')]
    
    if gurtel_ways.empty:
        print("No highways found matching 'Gürtel' suffix. Trying broader search...")
        gurtel_ways = ox.features_from_address(vienna_location, tags={"highway": True, "name": "~Gürtel"})
        gurtel_ways['name_clean'] = gurtel_ways['name'].astype(str).str.strip()
        gurtel_ways = gurtel_ways[gurtel_ways['name_clean'].str.endswith('Gürtel')]

    if gurtel_ways.empty:
        print("Could not find Gürtel highways.")
        return

    print(f"Found {len(gurtel_ways)} Gürtel segments.")

    # 2. Create 500m buffer
    # Use EPSG:31256 (Austria Austria GK94 / Gauss-Krüger zone 3) or EPSG:3035 (ETRS89-extended / LAEA Europe)
    # But for Vienna, EPSG:31256 or a UTM zone is better for meters. 
    # Let's use EPSG:32633 (WGS 84 / UTM zone 33N) which covers Vienna.
    
    gurtel_ways_proj = gurtel_ways.to_crs(epsg=32633)
    buffer_zone = gurtel_ways_proj.buffer(500).unary_union
    
    # 3. Get all highways intersecting the buffer
    print("Fetching surrounding highways...")
    # Get all highways in Vienna area (roughly) and filter by buffer
    # Or more efficiently: get highways within the bounding box of the buffer
    bbox = buffer_zone.bounds # (minx, miny, maxx, maxy)
    
    # Convert bbox back to lat/lon for osmnx
    from shapely.ops import transform
    import pyproj
    
    project_to_wgs84 = pyproj.Transformer.from_crs(epsg=32633, to_latlon=True).transform
    buffer_wgs84 = transform(project_to_wgs84, buffer_zone)
    minx, miny, maxx, maxy = buffer_wgs84.bounds
    
    # Query highways in the bbox
    # To be safe and thorough, we'll use the bounding box.
    # osmnx.features_from_bbox(north, south, east, west, tags)
    # bbox is (minx, miny, maxx, maxy) in lon/lat.
    # ox.features_from_bbox(north=maxy, south=miny, east=maxx, west=minx, ...)
    
    highway_tags = {"highway": True}
    # This might be a large query. Let's hope it works.
    highways_in_area = ox.features_from_bbox(maxy, miny, maxx, minx, tags=highway_tags)
    
    # Filter highways that intersect the buffer
    highways_in_area_proj = highways_in_area.to_crs(epsg=32633)
    highways_intersecting = highways_in_area_proj[highways_in_area_proj.intersects(buffer_zone)].copy()
    
    # Filter to only LineStrings and keep requested attributes
    # Note: osmnx features can be Points, Lines, Polygons. We want LineStrings.
    highways_intersecting = highways_intersecting[highways_intersecting.geometry.type.isin(['LineString', 'MultiLineString'])]
    
    # Clean up attributes
    cols_to_keep = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway', 'geometry']
    # Only keep columns that actually exist
    existing_cols = [c for c in cols_to_keep if c in highways_intersecting.columns]
    highways_final = highways_intersecting[existing_cols].copy()
    
    # Ensure LineString type (if MultiLineString, we'll explode or keep as is? 
    # The prompt says "highways layer as LineString". If it's a MultiLineString, 
    # we should probably explode it to individual LineStrings to satisfy the requirement.)
    if not highways_final.empty:
        highways_final = highways_final.explode(index_parts=False)
        # After explode, check if we have MultiLineStrings left (shouldn't if we only take LineString)
        # Actually, explode converts MultiLineString to LineString.
    
    # 4. Get PT routes
    print("Fetching PT routes...")
    pt_tags = {
        "route": ["bus", "tram", "subway", "train", "trolleybus", "light_rail"]
    }
    # Note: the prompt says "public-transport route relation". 
    # In OSM, these are usually 'route' tags on relations.
    # osmnx features_from_bbox will fetch these.
    
    pt_features = ox.features_from_bbox(maxy, miny, maxx, minx, tags=pt_tags)
    
    # Filter PT features that intersect the buffer
    pt_features_proj = pt_features.to_crs(epsg=32633)
    pt_intersecting = pt_features_proj[pt_features_proj.intersects(buffer_zone)].copy()
    
    # Keep requested attributes: ref, name, operator, route
    pt_cols_to_keep = ['ref', 'name', 'operator', 'route', 'geometry']
    existing_pt_cols = [c for c in pt_cols_to_keep if c in pt_intersecting.columns]
    pt_final = pt_intersecting[existing_pt_cols].copy()
    
    # Each route as a single feature, MultiLineString
    # OSM relations are already returned as single features (as MultiLineStrings/GeometryCollections)
    # We should ensure they are MultiLineStrings.
    if not pt_final.empty:
        # Some features might be GeometryCollections if they contain points and lines. 
        # We only want the linestring parts.
        def extract_linestrings(geom):
            if geom.geom_type == 'LineString':
                return geom
            elif geom.geom_type == 'MultiLineString':
                return geom
            elif geom.geom_type == 'GeometryCollection':
                lines = [g for g in geom.geoms if g.geom_type in ['LineString', 'MultiLineString']]
                if not lines:
                    return None
                from shapely.ops import unary_union
                return unary_union(lines)
            else:
                return None

        pt_final['geometry'] = pt_final['geometry'].apply(extract_linestrings)
        pt_final = pt_final[pt_final['geometry'].notnull()]
        
        # Ensure MultiLineString
        def ensure_multilinestring(geom):
            if geom.geom_type == 'LineString':
                return MultiLineString([geom])
            return geom
        
        pt_final['geometry'] = pt_final['geometry'].apply(ensure_multilinestring)

    # 5. Write to GeoPackage
    print("Writing to GeoPackage...")
    # Convert back to Austria standard projected system. 
    # prompt says "Austria's standard projected coordinate system". 
    # Usually this is EPSG:31256 or EPSG:3035. Let's use EPSG:31256.
    target_crs = "EPSG:31256"
    
    highways_final = highways_final.to_crs(target_crs)
    pt_final = pt_final.to_crs(target_crs)
    
    output_file = "vienna_network.gpkg"
    if os.path.exists(output_file):
        os.remove(output_file)
        
    highways_final.to_file(output_file, layer='highways', driver="GPKG")
    pt_final.to_file(output_file, layer='pt_routes', driver="GPKG")
    
    print(f"Success! Saved to {output_file}")

if __name__ == "__main__":
    main()

import osmnx as ox
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiLineString, LineString
import os
import pyproj
from shapely.ops import transform

def main():
    # 1. Fetch Gürtel highways
    # We'll search for "Gürtel, Vienna, Austria" to find the segments.
    print("Locating Gürtel...")
    try:
        # Use a larger distance to ensure we capture the full Gürtel ring
        # The Gürtel is a large ring around the center of Vienna.
        # We'll get all highways in the vicinity of 'Gürtel, Vienna' 
        # and then filter for those whose name ends in 'Gürtel'.
        gurtel_features = ox.features_from_address("Gürtel, Vienna, Austria", tags={"highway": True}, dist=5000)
    except Exception as e:
        print(f"Error fetching Gürtel: {e}")
        return

    if gurtel_features.empty:
        print("No Gürtel segments found.")
        return

    # Filter: ways whose name ends in 'Gürtel'
    # Note: name might be a list or a single string.
    def name_ends_with_gurtel(name):
        if isinstance(name, list):
            return any(str(n).strip().endswith('Gürtel') for n in name)
        return str(name).strip().endswith('Gürtel')

    gurtel_ways = gurtel_features[gurtel_features['name'].apply(name_ends_with_gurtel)].copy()

    if gurtel_ways.empty:
        print("No highways found with name ending in 'Gürtel'. Trying broader search (contains 'Gürtel')...")
        gurtel_ways = gurtel_features[gurtel_features['name'].astype(str).str.contains('Gürtel')].copy()

    if gurtel_ways.empty:
        print("Could not identify Gürtel segments.")
        return

    print(f"Found {len(gurtel_ways)} Gürtel segments.")

    # 2. Create 500m buffer
    # Use EPSG:32633 (UTM 33N) for Vienna for accurate buffering in meters
    gurtel_proj = gurtel_ways.to_crs(epsg=32633)
    buffer_zone = gurtel_proj.buffer(500).unary_union
    
    # Get bbox of buffer for querying
    bbox = buffer_zone.bounds # (minx, miny, maxx, maxy)
    
    # Convert bbox to lat/lon for osmnx
    # transformer from 32633 to 4326
    transformer = pyproj.Transformer.from_crs(epsg=32633, crs="EPSG:4326", always_xy=True)
    minx_lon, miny_lat, maxx_lon, maxy_lat = [transformer.transform(x, y) for x, y in [(bbox[0], bbox[1]), (bbox[2], bbox[3])]]
    # Wait, the transform above returns (lon, lat) for each point. 
    # bbox[0] is minx (lon), bbox[1] is miny (lat)... 
    # Actually bbox is (minx, miny, maxx, maxy) in the projected system.
    
    # Correct conversion:
    def project_point(x, y):
        return transformer.transform(x, y)
    
    # We need the lat/lon bbox
    p1 = project_point(bbox[0], bbox[1]) # (lon, lat)
    p2 = project_point(bbox[2], bbox[3]) # (lon, lat)
    
    # osmnx.features_from_bbox(north, south, east, west)
    # south = min(lat1, lat2), north = max(lat1, lat2), west = min(lon1, lon2), east = max(lon1, lon2)
    north = max(p1[1], p2[1])
    south = min(p1[1], p2[1])
    west = min(p1[0], p2[0])
    east = max(p1[0], p2[0])

    # 3. Fetch surrounding highways
    print("Fetching surrounding highways...")
    highway_tags = {"highway": True}
    highways_in_area = ox.features_from_bbox(north, south, east, west, tags=highway_tags)
    
    highways_in_area_proj = highways_in_area.to_crs(epsg=32633)
    # Filter by intersection with buffer
    highways_intersecting = highways_in_area_proj[highways_in_area_proj.intersects(buffer_zone)].copy()
    
    # Ensure we only have LineStrings/MultiLineStrings for the highway layer
    highways_intersecting = highways_intersecting[highways_intersecting.geometry.type.isin(['LineString', 'MultiLineString'])]
    
    # Select and keep requested attributes
    cols_to_keep = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway']
    existing_cols = [c for c in cols_to_keep if c in highways_intersecting.columns]
    # We must include 'geometry'
    highways_final = highways_intersecting[existing_cols + ['geometry']].copy()
    
    # Ensure they are LineStrings (explode MultiLineString)
    if not highways_final.empty:
        highways_final = highways_final.explode(index_parts=False)
        # After explode, ensure everything is LineString (some might still be MultiLineString if the geometry was weird)
        # But for the sake of the requirement "highways layer as LineString"
        highways_final = highways_final[highways_final.geometry.type == 'LineString']

    # 4. Fetch PT routes
    print("Fetching PT routes...")
    pt_tags = {
        "route": ["bus", "tram", "subway", "train", "trolleybus", "light_rail"]
    }
    pt_features = ox.features_from_bbox(north, south, east, west, tags=pt_tags)
    
    pt_features_proj = pt_features.to_crs(epsg=32633)
    pt_intersecting = pt_features_proj[pt_features_proj.intersects(buffer_zone)].copy()
    
    # Select and keep requested attributes: ref, name, operator, route
    pt_cols_to_keep = ['ref', 'name', 'operator', 'route']
    existing_pt_cols = [c for c in pt_cols_to_keep if c in pt_intersecting.columns]
    pt_final = pt_intersecting[existing_pt_cols + ['geometry']].copy()
    
    # Requirement: "each route as a single feature" and "MultiLineString"
    if not pt_final.empty:
        # Extract linestrings from potentially complex geometries (like GeometryCollections)
        def extract_linestrings(geom):
            if geom.is_empty:
                return None
            if geom.geom_type in ['LineString', 'MultiLineString']:
                return geom
            if geom.geom_type == 'GeometryCollection':
                lines = [g for g in geom.geoms if g.geom_type in ['LineString', 'MultiLineString']]
                if not lines:
                    return None
                from shapely.ops import unary_union
                return unary_union(lines)
            return None

        pt_final['geometry'] = pt_final['geometry'].apply(extract_linestrings)
        pt_final = pt_final[pt_final['geometry'].notnull()]
        
        # Ensure MultiLineString
        def ensure_multilinestring(geom):
            if geom.geom_type == 'LineString':
                return MultiLineString([geom])
            return geom
        
        pt_final['geometry'] = pt_final['geometry'].apply(ensure_multilinestring)

    # 5. Write to GeoPackage
    # Austria's standard projected system: EPSG:31256 (MGI / Austria GK Hayne)
    # or EPSG:3035 (ETRS89 / LAEA Europe). Let's use 31256 as requested.
    target_crs = "EPSG:31256"
    
    print(f"Writing to {target_crs}...")
    if not highways_final.empty:
        highways_final = highways_final.to_crs(target_crs)
    if not pt_final.empty:
        pt_final = pt_final.to_crs(target_crs)
    
    output_file = "vienna_network.gpkg"
    if os.path.exists(output_file):
        os.remove(output_file)
        
    if not highways_final.empty:
        highways_final.to_file(output_file, layer='highways', driver="GPKG")
        print(f"Layer 'highways' written.")
    
    if not pt_final.empty:
        pt_final.to_file(output_file, layer='pt_routes', driver="GPKG")
        print(f"Layer 'pt_routes' written.")
    
    print(f"Finished! File saved: {output_file}")

if __name__ == "__main__":
    main()

import osmnx as ox
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiLineString, LineString
import os
import pyproj
from shapely.ops import transform

def main():
    # Since we cannot use external APIs or download PBFs in this environment,
    # and the prompt implies I should be able to, I will provide the script 
    # that would work in a standard GIS environment.
    # However, I will attempt to simulate a small part of the process 
    # if I can find any way to get data, but the error "Connection refused" 
    # suggests the environment is air-gapped or highly restricted.
    
    # In a real scenario, the script below is the solution.
    
    # 1. Define the query parameters
    vienna_location = "Vienna, Austria"
    gurtel_search_term = "Gürtel"
    buffer_dist = 500 # meters
    target_crs = "EPSG:31256" # Austria standard

    print("Attempting to fetch data via OSMnx...")
    
    try:
        # 1. Get Gürtel highways
        # We look for highways in Vienna that contain 'Gürtel' in the name
        # We'll use a large area to be safe or a specific point
        gurtel_features = ox.features_from_address("Gürtel, Vienna, Austria", tags={"highway": True}, dist=5000)
        
        def name_ends_with_gurtel(name):
            if isinstance(name, list):
                return any(str(n).strip().endswith('Gürtel') for n in name)
            return str(name).strip().endswith('Gürtel')

        gurtel_ways = gurtel_features[gurtel_features['name'].apply(name_ends_with_gurtel)].copy()
        
        if gurtel_ways.empty:
             gurtel_ways = gurtel_features[gurtel_features['name'].astype(str).str.contains('Gürtel')].copy()

        if gurtel_ways.empty:
            print("No Gürtel segments found.")
            return

        # 2. Buffer
        gurtel_proj = gurtel_ways.to_crs(epsg=32633)
        buffer_zone = gurtel_proj.buffer(buffer_dist).unary_union
        
        # 3. Highway extraction
        bbox = buffer_zone.bounds
        # Conversion logic (simplified)
        transformer = pyproj.Transformer.from_crs(epsg=32633, crs="EPSG:4326", always_xy=True)
        p1 = transformer.transform(bbox[0], bbox[1])
        p2 = transformer.transform(bbox[2], bbox[3])
        north, south, west, east = max(p1[1], p2[1]), min(p1[1], p2[1]), min(p1[0], p2[0]), max(p1[0], p2[0])
        
        highways_in_area = ox.features_from_bbox(north, south, east, west, tags={"highway": True})
        highways_in_area_proj = highways_in_area.to_crs(epsg=32633)
        highways_intersecting = highways_in_area_proj[highways_in_area_proj.intersects(buffer_zone)].copy()
        highways_intersecting = highways_intersecting[highways_intersecting.geometry.type.isin(['LineString', 'MultiLineString'])]
        
        cols_h = ['name', 'highway', 'maxspeed', 'lanes', 'surface', 'oneway']
        existing_cols_h = [c for c in cols_h if c in highways_intersecting.columns]
        highways_final = highways_intersecting[existing_cols_h + ['geometry']].copy()
        if not highways_final.empty:
            highways_final = highways_final.explode(index_parts=False)
            highways_final = highways_final[highways_final.geometry.type == 'LineString']

        # 4. PT extraction
        pt_tags = {"route": ["bus", "tram", "subway", "train", "trolleybus", "light_rail"]}
        pt_features = ox.features_from_bbox(north, south, east, west, tags=pt_tags)
        pt_features_proj = pt_features.to_crs(epsg=32633)
        pt_intersecting = pt_features_proj[pt_features_proj.intersects(buffer_zone)].copy()
        
        cols_p = ['ref', 'name', 'operator', 'route']
        existing_cols_p = [c for c in cols_p if c in pt_intersecting.columns]
        pt_final = pt_intersecting[existing_cols_p + ['geometry']].copy()
        
        if not pt_final.empty:
            def extract_linestrings(geom):
                if geom.is_empty: return None
                if geom.geom_type in ['LineString', 'MultiLineString']: return geom
                if geom.geom_type == 'GeometryCollection':
                    lines = [g for g in geom.geoms if g.geom_type in ['LineString', 'MultiLineString']]
                    return MultiLineString(lines) if lines else None
                return None
            pt_final['geometry'] = pt_final['geometry'].apply(extract_linestrings)
            pt_final = pt_final[pt_final['geometry'].notnull()]
            pt_final['geometry'] = pt_final['geometry'].apply(lambda g: MultiLineString([g]) if g.geom_type == 'LineString' else g)

        # 5. Save
        highways_final = highways_final.to_crs(target_crs)
        pt_final = pt_final.to_crs(target_crs)
        
        output = "vienna_network.gpkg"
        if os.path.exists(output): os.remove(output)
        highways_final.to_file(output, layer='highways', driver="GPKG")
        pt_final.to_file(output, layer='pt_routes', driver="GPKG")
        print(f"Successfully created {output}")

    except Exception as e:
        print(f"An error occurred: {e}")
        print("Note: This script requires internet access to query OSM data via Overpass API.")

if __name__ == "__main__":
    main()