Wie steil sind Wiens Radstrecken?¶
Um diese Frage zu beantworten brauchen wir 1) die Radstrecken und 2) ein Geländeprofil.
Die Magistratsabteilung für Stadtvermessung (MA 41) bietet ein Geländemodell von Wien als Schichtenlinien in 1m Auflösung: https://www.data.gv.at/katalog/dataset/e45374c0-3e40-434c-92dc-203492ccee49 Laut MA 41 auch mit respektabler Genauigkeit:
In den Gebieten, in welchen das digitale Geländemodell zum überwiegenden Teil auf den Daten der Mehrzweckkarte basiert, beträgt die Höhengenauigkeit im Straßenraum plus/minus 10 Zentimeter und im Blockinnenbereich plus/minus 25 Zentimeter. In bewaldeten Bereichen, wo das digitale Geländemodell aus Laserscandaten abgeleitet wurde, beträgt die Höhengenauigkeit im flachen Gelände plus/minus 20 Zentimeter und im steilen Gelände plus/minus 1 Meter.
Uns interessieren alle Strecken die von Radfahrenden verwendet werden. Von der Stadt Wien gibt es dazu kein Datenset, lediglich die dedizierte Radinfrastruktur oder den Straßengraphen. Deswegen greifen wir auf die Daten der OpenStreetMap zurück, mit entsprechender Filterung:
- grundsätzlich beradelbare Strecken ohne jene die explizit als nicht beradelbar markiert sind (
bicycle~no|use_sidepath) - tendenziell nicht beradelbare Strecken die explizit als beradelbar markiert sind (
bicycle~yes|designated|permissive)
import pandas as pd
import osmnx as ox
import networkx as nx
import geopandas as gpd
from shapely import wkt, box
from shapely.ops import substring
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from tqdm import tqdm
street_graph = ox.graph_from_place("Vienna, Austria", custom_filter='["highway"~"living_street|residential|tertiary|tertiary_link|secondary|secondary_link|primary|primary_link|unclassified|cycleway"]["bicycle"!~"no|use_sidepath"]', retain_all=True)
other_graph = ox.graph_from_place("Vienna, Austria", custom_filter='["highway"~"service|pedestrian|track|footway|bridleway|path|construction"]["bicycle"~"yes|designated|permissive"]', retain_all=True)
merged_graph = nx.compose(street_graph, other_graph)
bikeable_lines = ox.graph_to_gdfs(merged_graph, nodes=False)
bikeable_lines = bikeable_lines.to_crs(epsg=31256)
contour_lines = pd.read_csv('datasets/hoehenlinien.csv')
contour_lines['geometry'] = contour_lines['SHAPE'].apply(wkt.loads)
contour_lines = gpd.GeoDataFrame(contour_lines, geometry='geometry', crs="EPSG:4326")
contour_lines = contour_lines.to_crs(epsg=31256)
Der OSM-Graph ist gerichtet, uns aber reicht eine Kante pro Richtung.
Dazu löschen wir alle Kanten bei denen die ID des Startknotens (u) kleiner ist als die des Endknotens (v) - außer: es gibt keine Kante in die andere Richtung (v -> u).
edges = bikeable_lines.copy()
edges.index = edges.index.droplevel('key') if 'key' in edges.index.names else edges.index
edge_set = set(edges.index)
mask = [(v, u) not in edge_set or u < v for u, v in edges.index]
bikeable_lines = bikeable_lines[mask]
Gibt es eh nur einfache Linien?
print(contour_lines.geom_type.unique())
print(bikeable_lines.geom_type.unique())
Schauen wir uns die Datensets mal an. Für ein bissi Steigung im Bereich der Berggasse.
berggasse_x, berggasse_y = 576, 341784
bbox_berggasse = box(berggasse_x - 500, berggasse_y - 500, berggasse_x + 500, berggasse_y + 500)
_, ax = plt.subplots()
bikeable_lines[bikeable_lines.intersects(bbox_berggasse)].plot(ax=ax, color="red", linewidth=1)
contour_lines[contour_lines.intersects(bbox_berggasse)].plot(ax=ax, color="blue", linewidth=1)
minx, miny, maxx, maxy = bbox_berggasse.bounds
ax.set_xlim(minx, maxx)
ax.set_ylim(miny, maxy)
plt.show()
Wie berechnen wir jetzt die Neigung der Radstrecken? Folgender Ansatz:
- Alle Punkte finden an denen sich Radstrecken und Höhenlinien überschneiden
- Radstrecken an diesen Punkten aufteilen
- Für jedes Segment die Steigung berechnen
Als erstes also die Überschneidungspunkte von Radstrecken und Höhenlinien finden:
contour_sindex = contour_lines.sindex
intersections = []
for bikeable_idx, bikeable_line in tqdm(bikeable_lines.geometry.items(), total=len(bikeable_lines)):
possibly_intersecting_contour_lines_index = list(contour_sindex.intersection(bikeable_line.bounds))
possibly_intersecting_contour_lines = contour_lines.iloc[possibly_intersecting_contour_lines_index]
for contour_line_idx, contour_line in possibly_intersecting_contour_lines.geometry.items():
if bikeable_line.intersects(contour_line):
intersection = bikeable_line.intersection(contour_line)
if intersection.geom_type == 'Point':
intersections.append({
'bikable_line_id': bikeable_idx,
'geometry': intersection,
'elevation': contour_lines.loc[contour_line_idx, 'HOEHE_ADRIA']
})
elif intersection.geom_type == 'MultiPoint':
for point in intersection.geoms:
intersections.append({
'bikable_line_id': bikeable_idx,
'geometry': point,
'elevation': contour_lines.loc[contour_line_idx, 'HOEHE_ADRIA']
})
intersections = gpd.GeoDataFrame(intersections, geometry='geometry', crs=bikeable_lines.crs)
Als nächstes die Radstrecken an den Höhenlinien-Schnittpunkten in Segmente aufteilen und deren Steigung berechnen.
Hier wird es schwieriger:
- Brücken und Tunnel verlaufen nicht auf dem Gelände das die Höhenlinien darstellen. Trotzdem überschneiden sich diese Linien.
Solche Strecken unterteilen wir nicht in Segmente und berechnen für sie auch keine Steigung. - Bei allen anderen Strecken entstehen durch die Segmentierung drei Arten von Segmenten:
- Segment berührt zwei Höhenlinien: easy, einfach Länge des Segments durch die Höhendifferenz der Höhenlinien teilen
- Segment berührt eine Höhenlinie: schwieriger - was ist die Höhendifferenz?
- Für den Segment-Endpunkt der keine Höhenlinie berührt die näheste Höhenlinie finden (h1)
- Die nächstnäheste Höhenlinie finden zu der eine Linie vom Endpunkt hinführt die nicht die erste gefundene Höhenlinie (h1) kreuzt
- Mit den Distanzen zu den gefundenen Höhenlinien deren Seehöhe linear zum Endpunkt des Segments interpolieren
- Die Höhendifferenz ist die Differenz zwischen der Seehöhe der vom Segment berührten Höhenlinie und der interpolierten Seehöhe
- Aus Höhendifferenz und Segmentlänge die Steigung berechnen
- Segment berührt keine Höhenlinie: gleicher Ansatz wie bei Segmenten die nur eine Höhenlinie berühren aber hier für beide Endpunkte näheste Höhenlinien suchen und interpolieren
from shapely.geometry import Point, LineString
def estimate_elevation(point, search_radius):
point = Point(point)
bbox = point.buffer(search_radius).bounds
candidate_contours = contour_lines.cx[bbox[0]:bbox[2], bbox[1]:bbox[3]].copy()
if len(candidate_contours) < 2:
return None
candidate_contours['distance'] = candidate_contours.geometry.distance(point)
candidate_contours = candidate_contours.sort_values('distance')
first_contour_line = candidate_contours.iloc[0]
for idx in range(1, len(candidate_contours)):
second_contour_line = candidate_contours.iloc[idx]
closest_point_on_second_contour_line = second_contour_line.geometry.interpolate(second_contour_line.geometry.project(point))
shortest_line_to_second_contour_line = LineString([point, closest_point_on_second_contour_line])
if not shortest_line_to_second_contour_line.crosses(first_contour_line.geometry) and not shortest_line_to_second_contour_line.within(first_contour_line.geometry):
total_distance = first_contour_line.distance + second_contour_line.distance
if total_distance == 0:
return (first_contour_line['HOEHE_ADRIA'] + second_contour_line['HOEHE_ADRIA']) / 2
elevation = first_contour_line['HOEHE_ADRIA'] + (first_contour_line.distance / total_distance) * (second_contour_line['HOEHE_ADRIA'] - first_contour_line['HOEHE_ADRIA'])
return elevation
return None
segments = []
for bikeable_idx, bikeable_line in tqdm(bikeable_lines.geometry.items(), total=len(bikeable_lines)):
if bikeable_lines.loc[bikeable_idx]['bridge'] == 'yes' or bikeable_lines.loc[bikeable_idx]['tunnel'] == 'yes':
segments.append({
'geometry': bikeable_line,
'elevation_delta': None,
'length': bikeable_line.length,
'slope': None
})
continue
current_intersections = intersections[intersections['bikable_line_id'] == bikeable_idx]
if len(current_intersections) == 0:
from_elevation = estimate_elevation(bikeable_line.coords[0], search_radius=500)
to_elevation = estimate_elevation(bikeable_line.coords[-1], search_radius=500)
elevation_delta = 0 if from_elevation is None or to_elevation is None else abs(from_elevation - to_elevation)
segments.append({
'geometry': bikeable_line,
'elevation_delta': elevation_delta,
'length': bikeable_line.length,
'slope': elevation_delta / bikeable_line.length if bikeable_line.length > 0 else 0
})
continue
current_intersections = current_intersections.copy()
current_intersections['distance'] = current_intersections.geometry.apply(lambda p: bikeable_line.project(p))
current_intersections = current_intersections.sort_values('distance').drop_duplicates(subset=['distance'], keep='first')
current_intersections = current_intersections.reset_index(drop=True)
if current_intersections.loc[0, 'distance'] > 0:
from_distance = 0
to_distance = current_intersections.loc[0, 'distance']
segment = substring(bikeable_line, from_distance, to_distance)
from_elevation = estimate_elevation(segment.coords[0], search_radius=500)
to_elevation = current_intersections.loc[0, 'elevation']
elevation_delta = 0 if from_elevation is None else abs(from_elevation - to_elevation)
segments.append({
'geometry': segment,
'elevation_delta': elevation_delta,
'length': segment.length,
'slope': elevation_delta / segment.length if segment.length > 0 else 0
})
for i in range(len(current_intersections) - 1):
from_elevation = current_intersections.loc[i, 'elevation']
to_elevation = current_intersections.loc[i+1, 'elevation']
from_distance = current_intersections.loc[i, 'distance']
to_distance = current_intersections.loc[i+1, 'distance']
segment = substring(bikeable_line, from_distance, to_distance)
elevation_delta = abs(from_elevation - to_elevation)
segments.append({
'geometry': segment,
'elevation_delta': elevation_delta,
'length': segment.length,
'slope': elevation_delta / segment.length if segment.length > 0 else 0
})
last_dist = current_intersections.loc[len(current_intersections)-1, 'distance']
if last_dist < bikeable_line.length:
from_distance = last_dist
to_distance = bikeable_line.length
segment = substring(bikeable_line, from_distance, to_distance)
from_elevation = current_intersections.loc[len(current_intersections)-1, 'elevation']
to_elevation = estimate_elevation(segment.coords[-1], search_radius=500)
elevation_delta = 0 if to_elevation is None else abs(from_elevation - to_elevation)
segments.append({
'geometry': segment,
'elevation_delta': elevation_delta,
'length': segment.length,
'slope': elevation_delta / segment.length if segment.length > 0 else 0
})
segments = gpd.GeoDataFrame(segments, geometry='geometry', crs=bikeable_lines.crs)
Obwohl die Höhenlinien eine Auflösung von 1 Meter haben gibt es einige Segmente die eine Höhendifferenz > 1 Meter haben.
Die Visualisierung dieser Segmente zeigt dass es sich hier meist um (fehlerhafte?) Unterbrechnungen in den Höhenlinien handelt.
weird_segments = segments[segments['elevation_delta'] > 1]
print(len(weird_segments))
for idx, row in weird_segments.iterrows():
segment = row.geometry
buffer_dist = 50
bbox = segment.buffer(buffer_dist).bounds
minx, miny, maxx, maxy = bbox
nearby_contour_lines = contour_lines.cx[minx:maxx, miny:maxy]
fig, ax = plt.subplots()
gpd.GeoSeries([segment]).plot(ax=ax, color='red', linewidth=3, label='Segment')
nearby_contour_lines.plot(ax=ax, color='blue', linewidth=1, alpha=0.7, label='Höhenlinie')
ax.set_xlim(minx, maxx)
ax.set_ylim(miny, maxy)
plt.legend()
plt.title(f"Segmentindex: {idx} mit Höhendifferenz {row['elevation_delta']:.2f} m")
plt.show()
Wir bilden den Steigungswerte jedes Segments auf einen Rot-Ton eines kontinuierlichen Farbschemas ab.
Besonders im Bereich des Wienerwalds gibt es viele Segmente mit Steigungswerten > 20%.
Solche Steigungswerte kommen im urbanen Gebiet nicht wirklich vor, wichtiger als die korrekte farbliche Repräsentation solcher Extremwerte ist uns die leichteren Schwankungen der Steigungswerte im Stadtgebiet hervorzuheben.
Deshalb verwenden wir das volle Spektrum des Farbschemas für Werte im Bereich von 0 - 15% an.
Steigungswerte darüber erhalten den gleichen Farbwert wie 15%.
normalize = mcolors.Normalize(vmin=0, vmax=0.15, clip=True)
colormap = plt.cm.get_cmap('Reds')
def slope_to_color(s):
if pd.isna(s):
return "#e4e4e4ff"
rgba = colormap(normalize(s))
return mcolors.to_hex(rgba)
segments['color'] = segments['slope'].apply(slope_to_color)
Vier color-stops für die Legende:
for i in range(4):
value = i * 0.2 / 4
color = mcolors.to_hex(colormap(normalize(value)))
print(f"Neigung {value:.2f}: {color}")
Für die Visualisierung brauchen wir außerdem die Summe der Länge aller Segmente als magische Konstante...
segments[segments['slope'] != None]['length'].sum() / 1000
... und die Segmente als GeoJSON.
segments.to_crs(epsg=4326).drop(columns=['elevation_delta']).to_file('web/segments.geojson', driver='GeoJSON', id_generate=True)