User:Maxbe/Satteldrehung nach Höhenlinien
Bemerkung
Das hier beschriebene Verfahren ist irgendwie ganz nett, hat aber ausserhalb einer kleinen Testumgebung nie wirklich funktioniert. Zum einen ergaben sich massive Performanceprobleme wegen der vielen geometrischen Abfragen (benachbarte Höhenlinen suchen, nach Abstand sortieren und nächstgelegenen Punkt suchen). Zum anderen ist es auch vom theoretischen Ansatz her etwas schräg: Man nimmt Höhenpunkte, abstrahiert daraus Höhenlinien und unternimmt dann grosse Anstrengungen, aus dieser Abstraktion wieder ein Modell der Höhen zu formen. Die aussichtsreichere Methode ist die nebenan beschriebene Berechnung aus den originalen Höhenwerten.
Worum gehts?
Sättel werden in vielen Karten mit einem Symbol dargestellt, das die Richtung des Sattels angibt. Mit "Richtung" ist hier die gedachte Linie von Tal zu Tal gemeint. Es gibt in OSM auch ein Tag dafür, leider wird direction=* nur selten eingetragen, noch seltener wird es auf unseren Karten dargestellt.
Die beste Lösung wäre natürlich, wir hätten bei jedem natural=saddle auch gleich ein direction-Tag, das dem Renderer die Ausrichtung verrät. Bis es soweit ist, kann der Renderer Sättel nur als ungerichtetes Symbol darstellen oder anderweitig versuchen, die Richtung abzuschätzen:
- Man könnte aus der Richtung der nächstgelegenen Straße die Richtung des Sattels abschätzen. Das funktioniert bestimmt auch bei grossen Pässen wie Brenner- oder Berninapass. Es scheitert aber bei kleinen Sätteln, die bestenfalls von Wanderwegen erschlossen werden. Letztere verbinden oft nicht Tal mit Tal, sondern sind lediglich ein markanter Punkt auf dem Weg den Grat entlang, oder ein Knick, wo der Weg vom Tal auf einen der benachbarten Gipfel schwenkt. Manche Sättel haben überhaupt keinen Weg, sondern sind nur als benamter Orientierungspunkt eingetragen.
- Man könnte aus den benachbarten Höhenwerten die Richtungen mit Steigung und Gefälle berechnen und daraus die Richtung abschätzen.
- Man könnte die Richtung aus den Höhenlinien abschätzen. Höhenlinien sind bei Renderern häufig anzutreffen, gelegentlich liegen sie sogar in der Datenbank neben den OSM-Objekten. Man sucht sich einfach die nächst niedrigere Höhenlinie und die Verbindungsline von Sattelpunkt zu dieser Linie gibt die Richtung an.
Sättel bei OSM
Der ideale Sattel...
... hat natürlich seine Richtung als direction=* eingetragen, aber um die ganz idealen Sättel gehts hier nicht.
Der fast ideale Sattel ohne direction-Tag liegt irgendwo in der Mitte eines sanduhrförmigen Musters der Höhenlinen. Er hat die Höhe korrekt und maschinenlesbar eingetragen. Hier könnte man tatsächlich die nächstniedrigere Höhenline suchen und die Richtung dorthin gibt die Richtung des Sattels an.
Der reale Sattel...
... hat kein direction-Tag, liegt irgendwo neben den Höhenlinien und hat keine oder eine zu den Höhenlinien unpassende Höhenangabe eingetragen. Gründe dafür gibt es viele:
- Die Höhenlinien sind in der Ebene oder in der Höhe verschoben, weil der Satellit nicht gut in enge Schluchten blicken konnte oder durch Reflexionen gestört war.
- Das Gelände ist so fein strukturiert, dass wesentliche Geländedetails vom groben Raster der Höhenmessung nicht erfasst wurden, z.B. falls der Grat sehr schmal ist oder einer der beiden Gipfel neben dem Sattel nur eine dünne Nadel ist.
- Der Mapper hat sich verschätzt, weil sein GPS im Gebirge nicht sehr exakt arbeitet, weil er von einem verschobenen Luftbild abgemalt hat oder ungenaue Daten importiert hat.
- Es wurde nicht der tiefste Punkt des Sattels eingetragen, sondern eine Stelle ein Stück den Hang rauf. Das kommt vor, weil auch die Wegebauer den Weg manchmal nicht an die tiefste Stelle legen, sondern ein bisschen oberhalb den Kamm überqueren. Dort wird dann das Hinweisschild mit Passname und Höhe aufgestellt oder der Aussichtsparkplatz mit Kiosk eingerichtet und diese Stelle wird auch in OSM als Sattel eingetragen.
Lösungsversuch
Abschätzung mit einer Höhenlinie
Um Probleme mit falschen (oder jedenfalls zu den Höhenlinien unpassenden) Höhenangaben zu umgehen, wird das ele-Tag des Sattels ignoriert. Lediglich die Lage wird ausgewertet.
- Für die Bestimmung der Höhe wird die nächstgelegene Höhenlinie herangezogen, in den beiden nächsten Bildern die 100m-Linie.
- Dann wird die übernächste Linie bestimmt und falls diese unter der eben gefundenen "Sattelhöhe" liegt, wird die Richtung zu dieser Höhenlinie als "Richtung des Sattels" festgelegt.
- Falls die übernächste Linie über der gefundenen Höhe liegt, geht es dort aufwärts und der Sattel zeigt senkrecht zur Verbindungslinie von Sattel und Höhenlinie.
Abschätzung mit mehreren Höhenlinien
Da der reale Sattel etwas neben der richtigen Position liegt, führ die Berücksichtigung nur einer Höhenlinie häufig zu falschen Richtungen. Berücksichtigt man mehrere Linien kann man einen Kompromiss finden.
- Man nimmt die 1. Linie zur Höhenbestimmung
- Sucht die 2. Höhenlinie und bestimmt die erste Näherung der Richtung
- Sucht die 3. Linie und korrigiert die Richtung. Der Einfluss dieser Linie ergibt sich aus (Distanz der 2. Linie zum Sattel)/(Distanz der 3. Linie zum Sattel)
- Sucht die 4. Linie, korrigiert wieder mit der Gewichtung (Distanz der 2. Linie zum Sattel)/(Distanz der 4. Linie zum Sattel)
- Sucht die 5. Linie ...
- ... macht das so lange, bis man eine bestimmte Anzahl von Linien durch hat, oder keine weitere Linie in einem bestimmten Umkreis findet oder die Gewichtung (Distanz der 2. Linie)/(Distanz der N. Linie) einen Grenzwert unterschreitet.
Probleme
(... viele ...)
- Das Verfahren ist nur in einem kleinen Umkreis vernünftig einsetzbar, es bevorzugt die nächstgelegenen Höhenlinien. Für Karten in kleinem Maßstab möchte man vielleicht eine eher "abstrakte Passrichtung" haben, die quer zu einem Gebirgszug liegt, statt tangential zur obersten Serpentinenkehre, die man erst im Maßstab 1:25000 sieht.
- Das Verfahren ist langsam. Abfragen wie "welche Höhenlinien liegen in 200m Umkreis, wo liegt der nächste Punkt an diesen Linien, in welcher Richtung und wie weit ist der entfernt?" kosten Rechenzeit und Plattenzugriffe.
- In flachem Gelände funktioniert es nicht, man muss mindestens 2, besser 4 oder 5 Höhenlinien finden.
- Bei stark verrauschten Höhenlinien funktioniert es nicht.
- Das Verfahren richtet sich bevorzugt nach den ansteigenden Hängen und ignoriert die Hälfte der unter dem Sattel liegenden Höhenlinien. Üblicherweise gibt es nämlich zwei Punkte, an denen der Sattelpunkt einer Höhenlinie nahe kommt, im Beispiel rechts kommen die 80 und 90m-Linie nach der Umrundung des linken Berges wieder am Sattel vorbei. Da aber nur der nächstgelegene Punkt berücksichtigt wird, werden die blau gekringelten Punkte ignoriert. Bei den höher gelegenen Linien passiert dieser Fehler nicht, die Berge links und rechts des Sattels haben je eigene Höhenlinien (bei Sätteln an Kraterrändern wäre der Effekt anders rum, aber Krater sind selten)
Ignorierte Näherungspunkte - Das ganze Verfahren hat keinerlei mathematische Grundlage, sondern ist das Ergebnis von Rumprobieren. Jeder einzelne Parameter (Suchradius, Anzahl der berücksichtigten Höhenlinien, Gewichtungsfaktoren ...) kann diskutiert und beliebig geändert werden, ohne "falscher" zu werden.
- Falls man die so ermittelte Richtung zusammen mit den OSM-Daten in eine Datenbank schreibt, sollte man sich die Lizenzen aller beteiligten Daten ansehen und seinen Anwalt fragen, ob man das darf und welche Daten unter welcher Lizenz man einem möglichen Interessierten zur Verfügung stellen muss. Man sollte überhaupt immer seinen Anwalt fragen, wenn man OSM-Daten zusammen mit Daten aus anderen Quellen nutzt und das Ergebnis darstellt. Bei der Verwendung von SRTM-Daten dürfte die Lösung aber einfach sein, die darf man ja beliebig weitergeben.
Beispiel
Rechts im Bild ein realer Sattel. Die blauen Punkte liegen auf den niedrigeren Höhenlinien, die roten auf den höheren. Die Striche oben geben die Ausrichtung des Sattels an, wenn man nur den einen Punkt berücksichtigen würde. Die schwarzen Linien unten zeigen die Richtung nach Berücksichtigung der Punkte 2-N.
Die Tabelle unten gibt die Richtung zur Höhenlinie und das Zwischenergebnis für jeden Schritt an. Man muss beachten, dass die Richtung um maximal 90° falsch sein kann. Statt im Schritt 3 um 122° nach links zu drehen wird als nur um 58° nach rechts gedreht. Zwischenergebnisse werden immer auf 0°-179° normiert, ein Sattel der nach Westen (=-90°) zeigt, zeigt genau so nach Osten (=+90°).
| Schritt | Höhe | Richtung | Entfernung | Differenz | Gewichtung | Ergebnis | |
|---|---|---|---|---|---|---|---|
| 1 | 1340 | Erster Schritt bestimmt die Höhe | |||||
| 2 | 1320 | 130° | 48m | - | 1 | 130° | erstes Setzen der Richtung |
| die nördliche 1320m-Linie wird ignoriert wegen Punkt Nr. 2 | |||||||
| die östliche 1340m-Linie wird ignoriert, weil sie gleich der Sattelhöhe ist | |||||||
| 3 | 1360 | 8° | 98m | 58° | 48/98 | 158° | Korrektur um 0.48*(188°-130°)=28° nach rechts |
| 4 | 1300 | 21° | 137m | 43° | 48/137 | 173° | 0.35*(201°-158°)=15° nach rechts |
| die südliche 1300m-Linie wird ignoriert wegen Punkt Nr. 4 | |||||||
| 5 | 1280 | 167° | 217m | 6° | 48/215 | 172° | eine kleine Korrektur nach links |
| 6 | 1380 | 4° | 221m | 12° | 48/221 | 175° | eine kleine Korrektur nach rechts |
| 7 | 1360 | 175° | 244m | 0° | 48/244 | 175° | dass sich hier nichts ändert, ist Zufall |
Programmierung
Die Funktion getsaddledirection(GEOMETRY,TEXT)
Die Funktion getsaddledirection(GEOMETRY,TEXT) nimmt einen Punkt (in EPSG:3857) als erstes Argument umd berechnet die Richtung des Sattels. Das zweite Argument ist die bereits gemappte Richtung als Zahl oder Text (z.b. '12' oder 'north' oder 'ssw').
Rückgabewert ist eine ganze Zahl zwischen 0 und 179 (0=Nord, 90=Ost, 179=fast Süd). Ein negativer Rückgabewert zeigt Fehler an, z.B. wenn keine Höhenlinien in der Gegend gefunden wurden. Bei Fehlern wird -135 zurückgegeben, was ein Renderer entweder auswerten könnte, oder einfach rendern, dann zeigt der Sattel nach Nordost.
Wer den Quelltext ohne Beschreibung haben will, findet ihn hier auch als eigene Datei
CREATE OR REPLACE FUNCTION getsaddledirection(point IN GEOMETRY,osmdirection IN TEXT) RETURNS INTEGER AS $$ BEGIN
Zuerst wird versucht, die gemappte Richtung zu interpretieren
IF (osmdirection ~ '^[0-9]+$') THEN direction=((osmdirection::INTEGER)+360)%180;
ELSEIF (osmdirection ~ '^[0-9]+\.[0-9]+$') THEN direction=((ROUND(osmdirection::FLOAT)::INTEGER)+360)%180;
ELSEIF (osmdirection='s' OR osmdirection='south' OR osmdirection='n' OR osmdirection='north')
THEN direction:=0;
ELSEIF (osmdirection='ssw' OR osmdirection='south-southwest' OR osmdirection='nne' OR osmdirection='north-northeast')
THEN direction:=22;
....
ELSE
Falls das nichts hilft, werden die Höhenlinien zur Auswertung herangezogen. Da diese Höhenlinien in einer anderen Datenbank liegen, wird die Verbindung mit dblink_connect_u() hergestellt und dann die Daten mit dblink() geholt.
IF ((dblink_get_connections() IS NULL) OR ('saddle_contours_connection' != ANY(dblink_get_connections()))) THEN
PERFORM dblink_connect_u('saddle_contours_connection', 'dbname=contours');
END IF;
Der String mit der Abfrage an die Höhenliniendatenbank wird zusammengebaut und soll ungefähr aussehen wie "SELECT geom,height,distance() FROM contours WHERE st_intersects(searcharea,geom) ORDER BY distance() LIMIT searchlimit"
querystring='SELECT wkb_geometry,height,ST_Distance(st_setsrid(wkb_geometry,3857),' ' ' || saddlepoint || ' ' '::geometry) as dist
FROM contours WHERE ST_Intersects(ST_Expand( || saddlepoint || ::geometry,' || SearchAreaMerc || '),st_setsrid(wkb_geometry,3857))
ORDER BY dist ASC LIMIT ' || SearchLimit;
Grosse Schleife um die Ergebnisse dieser Abfrage
<<getcontourloop>>
FOR result IN ( SELECT height::FLOAT,ST_ClosestPoint(st_setsrid(way,3857),st_setsrid(saddlepoint::geometry,3857)) AS cp,dist::FLOAT
FROM dblink('saddle_contours_connection',querystring)
AS t1(way geometry,height integer,dist float)
) LOOP
i:=i+1;
Erster Punkt? Dann ist das die Sattelhöhe
IF (i=1) THEN
saddleheight:=result.height;
saddleheight:=result.height;
ELSE
Alle anderen Punkte bestimmen die Richtung
IF (ABS(saddleheight-result.height)>MinDescent) THEN
SELECT CAST(d AS INTEGER) INTO thisdirection FROM (SELECT degrees(ST_Azimuth(saddlepoint,result.cp)) AS d) AS foo;
Abhängig ob es rauf oder runter geht wird die Richtung um 90° gedreht
IF (result.height<saddleheight) THEN thisdirection:=(360+thisdirection)%180; END IF;
IF (result.height>saddleheight) THEN thisdirection:=(360+thisdirection+90)%180; END IF;
Zweiter Punkt? Dann wird die aktuelle Richtung zum ersten Zwischenergebnis
IF(firstdistance<0) THEN
direction:=thisdirection;
firstdistance:=result.dist;
ELSE
Dritter bis N. Punkt? Dann wird der Differenzwinkel bestimmt (-90° bis +90°)
diffdirection=thisdirection-direction;
IF(diffdirection> 90) THEN diffdirection=180-diffdirection; END IF;
IF(diffdirection<-90) THEN diffdirection=diffdirection+180; END IF;
Und das neue Zwischenergebnis wird das bisherig + der gewichteten Differenz
direction:=(round(direction+diffdirection*(firstdistance/result.dist))::INTEGER)%180;
IF (direction<0) THEN direction:=direction+180; END IF;
END IF;
END IF;
END IF;
END LOOP getcontourloop;
Ende der Schleife
RETURN direction; END; $$ LANGUAGE plpgsql;
Installation
Einfach einmal obige Funktion laden:
psql gis < saddledirection.sql
Die Funktion macht allerdings Gebrauch von dblink_connect_u() und dblink(). Diese Funktionen stehen aus Sicherheitsgründen nur dem Datenbankverwalter postgis zur Verfügung, jedenfalls in der hier genutzten Version ohne Passwort in der Query. Nach Lektüre der Sicherheitsbedenken kann man entweder auch dem Konto die Rechte für dblink erteilen mit dem man die Datenbankabfragen ausführt (also z.B. "gis" oder "tirex"):
psql -d gis -c "GRANT EXECUTE ON FUNCTION dblink_connect_u(text) to username;" psql -d gis -c "GRANT EXECUTE ON FUNCTION dblink_connect_u(text,text) to username;"
oder man verzichtet darauf. Dann kann man allerdings während des Renderns diese Funktion nicht verwenden. Man könnte sie aber trotzdem nutzen, um z.B. alle paar Tage allen Sätteln die Richtung in die Datenbank einzutragen:
UPDATE planet_osm_point SET direction=getsaddledirection(way,direction) WHERE "natural" in ('saddle','col','notch') AND direction IS NULL;
Das muss man dann allerdings als "postgis" machen. Und man müsste sich überlegen, wie man mit korrekt getaggten Sätteln umgeht, die aber keine Zahl eingetragen haben, sondern z.B. "SSW". Die könnte man auch überschreiben (dann halt die Bedingung "direction IS NULL" weglassen), und durch Zahlen ersetzen, aber da könnte man auch Daten zerstören, die man später noch braucht.
Testen
Man sucht sich einfach einen Sattel raus, und testet die Funktion mit dessen Geometrie und verschiedenen Werten als 2. Parameter
SELECT osm_id,name,ele,getsaddledirection(way,NULL) FROM planet_osm_point WHERE "natural"='saddle' AND name='Glaswandscharte';
sollte irgendwas in grober Nord-Süd-Richtung liefern (2 oder 175 oder so)
SELECT osm_id,name,ele,getsaddledirection(way,'177') FROM planet_osm_point WHERE "natural"='saddle' AND name='Glaswandscharte';
liefert "177" zurück
SELECT osm_id,name,ele,getsaddledirection(way,'SSE') FROM planet_osm_point WHERE "natural"='saddle' AND name='Glaswandscharte';
liefert "157" für südsüdost
SELECT osm_id,name,ele,getsaddledirection(way,'links oben') FROM planet_osm_point WHERE "natural"='saddle' AND name='Glaswandscharte';
liefert das gleiche Ergebnis wie die erste Abfrage, weil "links oben" nicht interpretiert werden konnte und stattdessen die Höhenlinien befragt wurden.