Berechnung der längsten Gerade im Schienennetz der ÖBB mittels OpenStreetMap Daten

osm-laengste-gerade.sql

-- 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
;