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.
Leider ist der Download der gesamten Daten im OGD-Portal etwas umständlich gestaltet. Der abgebotene Downloadlink hat einen restriktiven Filter drinnen: https://data.wien.gv.at/daten/geo?service=WFS&version=1.1.0&request=GetFeature&typeName=ogdwien:HOEHENLINIEOGD&bbox=16.2%2C48.2%2C16.201%2C48.201%2CEPSG:4326&srsName=EPSG:4326&outputFormat=json&maxfeatures=5 Also maximal fünf Höhenlinien und nur jene innerhalb der definierten Box. Diese Filter einfach rauslöschen, absurd lange laden lassen und voilá - der Downlaod für das gesamte Datenset beginnt (560MB): https://data.wien.gv.at/daten/geo?service=WFS&version=1.1.0&request=GetFeature&typeName=ogdwien:HOEHENLINIEOGD&srsName=EPSG:4326&outputFormat=json
Als Radstrecken verwenden wir das Hauptradverkehrsnetz der Stadt Wien: https://www.data.gv.at/datasets/1ea3d3e8-fa07-4c37-af68-eb588d439de2 Da sind zwar nicht alle Wege enthalten die viele von uns auch mit dem Rad befahren aber die wichtigsten Strecken sollten wohl drinnen sein!
Ein drittes Datenset brauchen wir auch noch: die Flächenmehrzweckkarte (FMZK) der Stadt Wien: https://www.data.gv.at/datasets/7cf0da04-1f77-4321-929e-78172c74aa0b Da sind nämlich Brückenflächen drinnen und die brauchen wir um für Wege die über Brücken verlaufen keine falschen Neigungswerte zu berechnen.
import geopandas as gpd
cycling_lines = gpd.read_file('datasets/hauptradverkehrsnetz.json').to_crs(epsg=31256)
contour_lines = gpd.read_file('datasets/hoehenlinien.json').to_crs(epsg=31256)
fmzk = gpd.read_file('datasets/fmzk.gpkg').to_crs(epsg=31256)
bridge_areas = fmzk[fmzk['BW_BRK_ID'].notnull()]
Gibt es eh nur einfache Linien bzw. Polygone bei den Brücken?
print(bridge_areas.geom_type.unique())
print(contour_lines.geom_type.unique())
print(cycling_lines.geom_type.unique())
cycling_lines = cycling_lines.explode(index_parts=False).reset_index(drop=True)
Jetzt aber:
print(cycling_lines.geom_type.unique())
Schauen wir uns die Datensets mal an. Für ein bissi Steigung im Bereich der Berggasse.
import matplotlib.pyplot as plt
from shapely import box
berggasse_x, berggasse_y = 576, 341784
bbox_berggasse = box(berggasse_x - 500, berggasse_y - 500, berggasse_x + 500, berggasse_y + 500)
_, ax = plt.subplots()
cycling_lines[cycling_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()
Vor wir mit der Neigungsberechnung anfangen filtern wir mal alle Streckenabschnitte raus die über Brücken verlaufen - für die können wir nämlich keine Neigung berechnen (die Höhenlinien beziehen sich auf den Boden unter der Brücke).
Dazu die Radstrecken an den Überschneidungspunkten mit Brückenflächen aufsplitten und Strecken die über Brückenflächen verlaufen getrennt ablegen.
from shapely.ops import split
from tqdm import tqdm
import pandas as pd
bridge_tagged_cycling_lines = []
for idx, line in tqdm(cycling_lines.geometry.items(), total=len(cycling_lines)):
possible_idx = list(bridge_areas.sindex.intersection(line.bounds))
possible_bridges = bridge_areas.iloc[possible_idx]
overlapping_bridges = possible_bridges[possible_bridges.intersects(line)]
if overlapping_bridges.empty:
new_row = cycling_lines.loc[idx].copy()
new_row['bridge'] = False
bridge_tagged_cycling_lines.append(new_row)
else:
merged_bridge_area = overlapping_bridges.union_all()
splitted = split(line, merged_bridge_area)
for segment in splitted.geoms:
new_row = cycling_lines.loc[idx].copy()
new_row['geometry'] = segment
intersection = segment.intersection(merged_bridge_area)
fraction_inside = 0
if segment.length > 0:
fraction_inside = intersection.length / segment.length
new_row['bridge'] = fraction_inside > 0.8
bridge_tagged_cycling_lines.append(new_row)
cycling_lines = gpd.GeoDataFrame(bridge_tagged_cycling_lines, geometry='geometry', crs=cycling_lines.crs).reset_index(drop=True)
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:
intersections = []
for bikeable_idx, bikeable_line in tqdm(cycling_lines.geometry.items(), total=len(cycling_lines)):
possibly_intersecting_contour_lines_index = list(contour_lines.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=cycling_lines.crs)
Als nächstes die Radstrecken an den Höhenlinien-Schnittpunkten in Segmente aufteilen und deren Steigung berechnen.
Hier wird es schwieriger:
- Radstrecken-Segmente die über Brücken verlaufen ignorieren wir. Für sie wird kein Neigungswert berechnet.
- Bei allen anderen Strecken entstehen durch die Segmentierung an den Schnittpunkten mit den Höhenlinien 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 nur 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 gar 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
from shapely.ops import substring
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(cycling_lines.geometry.items(), total=len(cycling_lines)):
if cycling_lines.loc[bikeable_idx]['bridge']:
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=cycling_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%.
import matplotlib.colors as mcolors
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
... außerdem die Straßennamen ...
import json
unique_strnam = cycling_lines['STRNAM'].dropna().unique().tolist()
with open('web/streetnames.json', 'w', encoding='utf-8') as f:
json.dump(unique_strnam, f, ensure_ascii=False, indent=2)
... und die Segmente als GeoJSON.
segments.to_crs(epsg=4326).drop(columns=['elevation_delta']).to_file('web/segments.geojson', driver='GeoJSON', id_generate=True)