Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Berechnung der längsten Gerade im Schienennetz der ÖBB mittels OpenStreetMap Daten
-- alle Zuglinien als Relation abfragen
with recursive get_lines as (
select osm_id, way, row_number() over() as tmp_id from at_osm.planet_osm_line
where route in ('train', 'railway')
),
-- die einzenen Punkte der Relationen ermitteln
dump_points as (
select tmp_id, osm_id, path[1] as path_id, geom from get_lines as a inner join lateral st_dumppoints(way) as b on true
),
-- den Winkel zwischen den einzelnen Punkten berechnen
calculate_azimuth as (
select a.tmp_id, a.osm_id, a.path_id,
round(st_azimuth(a.geom, b.geom)::numeric, 1) as azimuth,a.geom
from dump_points as a inner join dump_points as b on (a.osm_id=b.osm_id and a.path_id+1=b.path_id and a.tmp_id=b.tmp_id)
),
-- herausinden, ob der Winkel davor auch gleich ist, wie der aktuelle. Falls nicht, dann handelt es sich um einen Startpunkt
calculate_lag as (
select tmp_id, osm_id, azimuth, geom, path_id,
case when (azimuth IS DISTINCT FROM lag(azimuth) over (partition by tmp_id, osm_id order by path_id)) then true end as
calculate_recursion
from calculate_azimuth
),
-- rekursiv die Ausdehnung der einzelnen Startpunkt berechnen
calculate_longest_line as (
select tmp_id, osm_id, path_id as start_path_id, path_id, geom, geom as end_point, azimuth from calculate_lag where calculate_recursion=true
union all
select a.tmp_id, a.osm_id, b.start_path_id, a.path_id, b.geom, a.geom, a.azimuth from calculate_lag as a, calculate_longest_line as b
where a.osm_id=b.osm_id and a.tmp_id=b.tmp_id and a.path_id=b.path_id+1 and a.calculate_recursion is null
),
-- maximale Ausdehnung berechnen
calculate_path_range as (
select tmp_id, osm_id, start_path_id, max(path_id) as end_path_id from calculate_longest_line group by 1,2,3
),
-- die Länge ermitteln. Einerseits als korrekte Länge im geography-Typ und andrerseits in der Kartenprojektion
calculate_length as (
select a.tmp_id, a.osm_id, a.start_path_id, end_path_id, azimuth,
ST_Distance(ST_Transform(geom, 4326)::geography, ST_Transform(end_point, 4326)::geography) as length_geography,
ST_Distance(geom, end_point) as length_geom
from calculate_longest_line as a,
calculate_path_range as b where a.osm_id=b.osm_id and a.tmp_id=b.tmp_id and a.start_path_id=b.start_path_id
and b.end_path_id=a.path_id
),
-- die maximale Länge pro Relation berechnen
get_max_length as (
select tmp_id, osm_id, max(length_geography) as max_length_geography,
max(length_geom) as max_length_geom from calculate_length
group by 1,2
)
-- und sich die Resultate geordnet nach Länge ausgeben
select * from get_max_length order by max_length_geography desc
;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment