User:Maxbe/Kartenversuch
Vorbemerkung und Warnung
Ich habe die Karte nur zum Spielen gemacht, um halt auch mal eine Karte zu rendern. Und nur für eigene Zwecke, deshalb reicht mir auch eine kleine Gegend (9°Ost bis 15°Ost und 46°Nord bis 49.1°Nord, erzeugt aus dem europe-Extrakt der Geofabrik) südlich und östlich von München sowie die Alpen südlich von München.
Ich kann kaum SQL und ich mag eigentlich keine Datenbanken. Viele Abfragen im preprocessing sind echte quick'n dirty-Programmierung, wo ganze Tabellen durchsucht werden statt eines Indexes. Insbesondere wird häufig der spatial index von PostGIS ignoriert. Sowas ist pfui, das funktioniert mit winzigen Gebieten, skaliert aber nicht mit größeren Datenmengen. Trotzdem erzähle ich aber gern, wie ich die Karte mache. Vielleicht kann ja ein besserer Datenbanker was damit anfangen...
Worum geht es?
Ziel war eine Karte, die als Hintergrund für die typische "Das ist mein Laden"- oder "das ist der Wanderweg auf meinen Lieblingsberg"-Karte im Web dienen kann. Dazu werden relativ wenige POIs gebraucht (z.B. nicht der Konkurrenzladen in der Nachbarschaft), aber alle wichtigen Straße, Sehenswürdigkeiten, Bahnhöfe und Ortsnamen sollten drauf sein. Für die Berge sollten Höhenlinien und Schummerung drauf sein. Zu besichtigen ist die Karte hier: http://geo.dianacht.de/topo/
Das ganze sollte mit dem mapserver (Version 6) von mapserver.org laufen. Weil damit kenne ich mich aus und kann recht einfach Grafikexporte in beliebiger Größe
erzeugen. Betriebssystem ist Debian Lenny 32 Bit, die Datenbank sollte nicht mehr als 3GB auf der Platte brauchen und auf einen vServer der 10-Euro-Klasse 1-2 gleichzeitige User bedienen können.
Vorbereitung
OSM-Datenbank
Man braucht eine laufende PostGIS-DB mit OSM-Daten. Da habe ich mich an http://trac.osgeo.org/mapserver/wiki/RenderingOsmData gehalten und fürs die gelegentlichen Updates an http://wiki.openstreetmap.org/wiki/DE:HowTo_minutely_hstore
Das Einspielen und Updaten der Datenbank ist dort auch gut beschrieben. Letztendlich hatte ich eine Datenbank der gewünschten Gegend namens "osm" mit Koordinaten in EPSG 900913.
Höhendaten
Für den Import der Höhendaten braucht man die Tools von http://www.gdal.org/ und für die Schummerung die perrygeo-Tools von http://www.perrygeo.net/wordpress/ (zumindest das Programm "hillshade").
Die Gebrauchsanweisung für das Einspielen der Höhendaten steht bei http://wiki.openstreetmap.org/wiki/Shaded_relief_maps_using_mapnik
Dummerweise ist dort das Einspielen für OSM-Datenbanken mit EPSG 4326-Daten beschreiben, weswegen ich die Daten noch nachträglich konvertiere. Die eigentlichen Daten nehme ich in Form eines geoTIFFs von http://www.opendem.info die sind schön geglättet und geflickt an den Löchern der originalen NASA-Daten.
Das Einspielen sieht dann bei mir so aus:
Runterladen und auspacken
wget "http://koenigstuhl.geog.uni-heidelberg.de/~mover/srtm_germany_dtm.zip" unzip srtm_germany_dtm.zip mv srtm_germany_dtm.tif srtm.tif rm srtm_germany_dtm.zip
Höhenlinien im 25m-Abstand erzeugen und in die DB namens "srtm" einspielen:
psql srtm -c "drop table contours;" gdal_translate -projwin 10 48.655 13 47 -of GTiff -co "TILED=YES" -a_srs "+proj=latlong" srtm.tif input_adapted.tif gdal_contour -i 25 -snodata 32767 -a height input_adapted.tif input_adapted.shp shp2pgsql -p -I -g way input_adapted.shp contours | psql -q srtm shp2pgsql -a -g way input_adapted.shp contours | psql -q srtm
Konvertieren in spherical mercator
psql srtm -c "UPDATE geometry_columns SET SRID=900913 WHERE f_table_name='contours';" psql srtm -c "alter table contours drop constraint enforce_srid_way;" psql srtm -c "UPDATE contours set way=Transform(ST_SetSRID(way,4326),900913);"
In ein Zusätzliches Feld den "Typ" der Höhenlinie eintragen. "25", "50" oder "100", das erleichtert später das rendern von z.B. nur den 100er-Linien
psql srtm -c "alter table contours add column typ text;" psql srtm -c "update contours set typ='25';" psql srtm -c "update contours set typ='50' where CAST(height as int)%50=0;" psql srtm -c "update contours set typ='100' where CAST(height as int)%100=0;"
Die Schummerung wird ebenfalls aus der oben erzeugten Datei input_adapted.tif erzeugt, nach
http://wiki.openstreetmap.org/wiki/Hillshading_with_Mapnik
gdalwarp -of GTiff -co "TILED=YES" -srcnodata 32767 -t_srs "+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m" \ -rcs -order 3 -tr 24 24 -wt Float32 -ot Float32 -wo SAMPLE_STEPS=100 -multi input_adapted.tif warped.tif hillshade warped.tif hillshade.tif -z 3 cp hillshade.tif /home/www/wms/schummerung/hillshade.tif
Jetzt hab ich also eine DB "osm mit den Openstreetmapdaten, eine DB "srtm" mit den Höhen und ein TIFF mit den Graustufen der Schummerung.
Wasser ausschneiden
Die Höhenlinien überschneiden sich gelegentlich mit Wasserflächen. Das kommt in der Natur nicht oft vor, weil Wasser oft sehr waagerecht liegt. Um diese Auffälligkeit der schlechten Höhendaten zu verbergen, werden alle Höhenlinien an größeren Wasserflächen abgeschnitten:
$db=pg_pconnect('host=localhost port=5432 dbname=osm user=osm password=osm') or die("Keine DB, sonst alles in Ordnung");
$query=" select st_astext(way) as w,name as n from osm_polygon where \"natural\"='water' and way_area>10000;";
$result=pg_query($db,$query);
while($wasser=pg_fetch_array($result)) {
echo "-- ".$wasser['n']."\n";
echo "delete from contours where st_contains(ST_GeomFromText('".$wasser['w']."',900913),way);\n";
echo "update contours set way=ST_Difference(way,ST_GeomFromText('".$wasser['w']."',900913))
where st_intersects(way,ST_GeomFromText('".$wasser['w']."',900913));\n";
}
Das gleiche kann man natürlich noch mit landuse=basin oder reservoir oder riverbank machen. Bei mir kam noch eine Fehlermeldung "... violates check constraint enforce_geotype_way ...". Die kommt daher, dass irgendwo beim Import ein constraint eingebaut wurde, dass Höhenlinien zwingend MULTILINESTRING sein müssen, ST_Difference aber auch simple LINESTRING erzeugen kann. Zum Rendern ist das egal, also hab ich das constraint gelöscht ("ALTER TABLE contours DROP CONSTRAINT enforce_geotype_way;"). Vorsichtige Menschen können ja nach dem Programmlauf noch "select geometrytype(way) from contours" ausführen und nach ungewöhnlichen Typen suchen.
Bessere Höhendaten für Österreich und Bayern
Österreich hat 10-Meter-Höhendaten unter cc-by-sa-Lizenz veröffentlicht. Der Einbau lohnt sich und wurde freundlicherweise von Andreas Binder im Forum beschrieben: http://forum.openstreetmap.org/viewtopic.php?id=30626
Die Bayerische Vermessungsverwaltung hat 50-Meter-Daten unter cc-by-Lizenz veröfentlicht. Das ist zwar nicht ganz so toll wie bei den Österreichern, aber gut genug, um es einzubauen. Gegenüber den SRTM-Daten sind die Messpunkte nur unwesentlich dichter, dafür gibt es weniger Stellen, wo der Satellit nicht hinsah, geblendet wurde oder wild phantasiert hat.
Die Daten gibt es bei http://www.vermessung.bayern.de/opendata zum runterladen und haben ein Format, das eher unhandlich ist. Jeweils 10x10km in einer Datei der Form (rechtswert hochwert hoehe)
~/ head 437_556_t50dgm.txt 4370000 5569950 363.6 4370050 5569950 360.2 4370100 5569950 356.1 ....
Die steckt man alle in eine Datei mit genau einem Leerzeichen als Trenner:
cat *dgm.txt | awk '{print $1" "$2" "$3;}' > baydgm50.unsorted_xyz
oder fuer kleinere Bereiche
cat *_5[1234]*dgm.txt | awk '{print $1" "$2" "$3;}' > baydgm50.unsorted_xyz
Nachsehen, wie gross der Bereich ist
cat baydgm50.unsorted_xyz | awk 'BEGIN{FS=" ";xma=yma=-1;xmi=ymi=100000000}NR>1{if($1>xma){xma=$1;}if($1<xmi){xmi=$1;}if($2>yma){yma=$2;}if($2<ymi){ymi=$2;}}END{print xmi" "xma" "ymi"}
Fuer ganz Bayern kommen dann diese Werte raus: 4279050 4636750 5236050 5605500
Mit diesen Werten schreibt man sich dann ein C-Programm und wandelt baydgm50.unsorted_xyz in baydgm50.xyz um. In der xyz sind die Werte nach y sortiert und haben dort, wo es keine Werte gibt den Wert -999999 stehen. Die bbox muss nämlich für die Weiterverarbeitung mit "nodata"-Werten aufgefuellt sein.
Die xyz-Datei kann man jetzt in das handhabbare tif überführen:
gdal_translate -a_nodata -999999.0 baydgm50.xyz baydgm50.tif
Das Ergebnis ist ein tif mit den Ecken 4279025,5236025 4636775, 5605525 und 50x50 Pixelgroesse und Koordinaten in EPSG:31468 (wobei letzteres nicht in der Datei steht).
Dieses tif kann man dann in Mercator oder EGSG:4326 umprojizieren und weiter verarbeiten
gdalwarp -of GTiff -co "TILED=YES" -s_srs "EPSG:31468" -t_srs "EPSG:900913" -srcnodata -999999 baydgm50.tif baydgm50_merc.tif
Die SRTM-Daten sollte man aber trotzdem behalten, oder seinen Anwalt fragen, was lizenzrechtlich passiert, wenn man z.B. die Dominanz von Gipfeln oder die Ausrichtung von Sätteln (s.u.) mit cc-by-sa-Daten rechnet und die so errechneten Werte in seine odbl-Datenbank zurückschreibt...
Preprocessing
Um nicht bei jedem Rendervorgang mühsam die POIs aus der Datenbank zu extrahieren, die mich interessieren, wird eine eigene Tabelle mit den POIs angelegt ("osm_poi"). Dabei wird auch gleich das lästige Problem beseitigt, dass viele Objekte doppelt in der DB vorhanden sind, einmal als Fläche und einmal als Node innerhalb derselben. Einer der beiden wir weggeworfen.
Das gesamte Preprocessing gibt es auf der oben verlinkten Seite zum runterladen. Zum nachbauen dürfte es nicht gut geeignet sein, aber vielleicht zum nachsehen, wie ein bestimmtes Icon auf die Karte kommt. Im wesentlichen besteht es aus einer langen Liste von Befehlen der Art
select way,osm_id,leisure,sport,name,amenity from osm_point where ( leisure = 'pitch' or sport is not null )
union (select centroid(way) as way,osm_id,leisure,sport,name,amenity from osm_polygon
where (leisure = 'pitch' or sport is not null)
and not exists (select way from osm_point where
(osm_point.leisure=osm_polygon.leisure or osm_point.sport=osm_polygon.sport)
and ST_Contains(osm_polygon.way,osm_point.way)))
und dem anschliessenden kopieren dieses Nodes in osm_poi (in diesem Fall werden Sportplätze gesucht, falls es einen sowohl als Node als auch als Polygon gibt, gewinnt der Node).
Es gibt noch ein paar weitere Aktionen im preprocessing, das 1x täglich durchgeführt wird. Die meisten hab ich aus Diskussionen im Forum entnommen.
S- und U-Bahnhöfe
Idee aus dem Forum: http://forum.openstreetmap.org/viewtopic.php?id=13907
Zum einen wollte ich alle U-Bahnstationen mit einem "U" und alle S-Bahnen mit "S" beschriftet haben. Keine ganz leichte Aufgabe, weil sowas nur unvollständig in den Daten steht. Aber in kleinen Schritten kommt man hin:
Verwendet wird das Feld "service". Das braucht bei Bahnhöfen keiner und ist gut geeignet, die Begriffe "station", "ubahn" und "sbahn aufzunehmen".
S-Bahnhöfe sind Stationen, die zu Bahnstrecken gehören, die der MVV betreibt:
update osm_point set service='station' where railway='station'; update osm_point set service='sbahn' where railway='station' and exists \ (select tags from osm_rels where 'n'||osm_id=ANY(members) and \ ARRAY['route','train']<@tags and ARRAY['network','MVV']<@tags);
U-Bahnhöfe sind Stationen, die zu U-Bahnstrecken gehören, die der MVV betreibt:
update osm_point set service='ubahn' where railway='station' and exists \ (select tags from osm_rels where 'n'||osm_id=ANY(members) and \ ARRAY['route','subway']<@tags and ARRAY['network','MVV']<@tags);
Manche sind aber nicht Teil der Relation, sondern Teil der Schiene:
for i in `psql osm -c "select nodes from osm_ways where exists (select osm_id \
from osm_line where id=osm_id and railway='subway');" | \
sed 's/{//g; s/}/,/g; s/,/ /g; s/[^0-9 ]//g'` ; do
echo $i
psql osm -c "update osm_point set service='ubahn' where osm_id='"$i"' and railway='station' ;"
done
Und jetzt der Holzhammer für alles was übrigbleibt. Ich habe eine Tabelle "bahnhof" mit allen S- und U-Bahnhaltestellen:
update osm_point as o1 set service='sbahn' where railway='station' and service='station' \ and exists(select name from bahnhof where bahn='s-bahn' and name=o1.name) and not exists \ (select * from osm_point as o2 where railway='station' and service='sbahn' and o2.name=o1.name);
Sowas funktioniert natürlich nur, wenn man in seinem Gebiet nur einen Betreiber für S-Bahnen hat und nur eine Stadt mit U- und S-Bahn (an dieser Stelle nochmal der Hinweis auf die Vorbemerkung...).
Dominanz von Gipfeln
Idee aus dem Forum: http://forum.openstreetmap.org/viewtopic.php?id=14210
zur besseren Übersichtlichkeit von Bergregionen in kleineren Maßstäben wird bei Gipfeln noch die Entfernung zum nächstgelegenen höheren Gipfel bestimmt und im unbenutzten Feld "population" abgelegt. Das ganze dient dazu, z.B. im Zoomlevel 10 nur die "wichtigsten" Berge einer Region darzustellen. Die Höhe alleine als Mass dafür ist nicht sonderlich hilfreich, weil hohe Berge oft nur unbedeutende Nachbarn noch höherer sind. Niedrige freistehende Berge sind oft wichtiger.
Im Beispiel rechts gewinnt die Leilachspitze im Zoom 10. Die ist auch der höchste Berg dieser Gegend. Bei Zoomlevel 12 wird die Sulzspitze angezeigt, die höhere Lachenspitze dagegen noch nicht. Das liegt daran, dass die Lachenspitze dank der geringen Entfernung zur höheren Roten Spitze eine geringere Dominanz hat. Im Zoomlevel 13 gibt es keine Generalisierung nach Dominanz mehr, da entscheidet dann einfach die Höhe.
Möglicherweise würde den Fall "Lachenspitze gegen Sulzspitze" ein Betrachter vor Ort anders sehen. Häufig passt die Sortierung aber meines Erachtens schon gut. Ausnahmen sind Berge, die aus besonderen Gründen als prominent angesehen werden (wegen ihrer Form, oder weil dort eine Bahn rauf geht z.B.)
Dominanz mit Abstand zur Höhenline
In Wirklichkeit wird die Dominanz nicht mit dem Abstand vom nächsten Gipfel definiert, sondern mit dem Abstand zum nächsthöheren Punkt. Wer also Zeit hat, mit mehr Punkten zu rechnen, sollte auch die Höhendaten ausserhalb der Gipfel berücksichtigen.
Den Unterschied sieht man am nebenstehenden Beispiel. Im oberen Teil ist der Waxenstein der dominanteste Berg in der Gegend und stellt seinen kleineren Nachbarn, die Demelspitze in den Schatten. Sein grösserer Nachbar, das Brauneck, hat ebenfalls eine kleinere Dominanz als der Waxenstein, weil es nur 1200m vom höheren Stangeneck entfernt ist.
Berechnet man allerdings die Entfernung nach den Höhendaten, werden Waxenstein und Demelspitze sehr undominant (sind sie tatsächlich, die Demelspitze ist ein unscheinbarer Felszacken mitten im steilen Wald, der Waxenstein ein wenig bedeutender Vorgipfel). Das Brauneck verliert zwar auch etwas, weil der Schnittpunkt mit dem Stangeneck ein Stück "vor" dessen Gipfel liegt, allerdings ist es nun wesentlich dominanter als der Waxenstein und wird früher dargestellt.
Angesichts der eher unbefriedigenden Qualität der SRTM- oder ASTER-Höhendaten, schlage ich noch ein paar Meter zur Gipfelhöhe dazu. Ausserdem werden Linien im Umkreis von 100m um den Gipfel nicht berücksichtigt. Sonst passiert es nämlich, dass man bei der Berechnung der Dominanz auf die eigene Höhenlinie um den Gipfel stösst, nur weil die Höhe um ein paar Meter falsch ist oder der Gipfel ein bisschen neben der Stelle steht, wo laut Höhenlinien der höchste Punkt sein soll. Zur Optimierung begrenze ich die Suche nach der höheren Höhenlinie auf die Entfernung zum nächsthöheren Gipfel (die Berechnung oben). Sollte in diesem Umkreis keine Höhenlinie gefunden werden, ist irgend etwas kaputt an den Höhendaten...
Die Berechnung wird mit einem Programm aus der OpenTopoMap berechnet, dessen Beschreibung man hier nebenan findet.
Dominanz von Kneipen
Das Problem bei Gastronomiebetrieben ist, dass ich frei stehende Biergärten und abgelegene Bergasthöfe (auch wenn diese nicht als alpine_hut o.ä. eingetragen sind) gerne bei Zoom 15 oder 16 rendern würde. Restaurants in der Münchner Innenstadt aber bestenfalls erst bei Zoom 18, und dann eher nicht so deutlich. Deshalb wird der Gastronomie ebenfalls eine "population" zugewiesen. Dieser Wert gibt die Anzahl der Gastronomiebetriebe in der 1000-"Meter"-bbox an und schwankt zwischen 1 (einzelner Biergarten), 2 (Gasthaus mit einzeln getaggtem Biergarten) und 300 (München, Glockenbachviertel).
Beschriftung von Gebirgen
- Available in English
Idee aus dem Forum: http://forum.openstreetmap.org/viewtopic.php?id=16094
Das ist noch nicht so toll gelöst und auch nicht im preprocessing dauerhaft untergebracht. Der Herr Netzwolf hat eine ganze Menge Relationen der Form (name=Wetterstein;place=region;region:type=mountain_area;type=multipolygon) erstellt. Die fassen im wesentlichen Bachläufe oder andere Grenzlinien zusammen und umreißen jeweils ein Gebirge (z.B. http://www.openstreetmap.org/browse/relation/2127978).
- Umriss der Relation holen
- In allen vier Himmelsrichtungen nach dem Punkt suchen, der am weitesten vom Centroid der Relation entfernt ist (wobei beispielsweise ein Punkt zugleich der weitest entfernte im Norden als auch im Westen sein kann)
- Einen Kreis durch die beiden gefundenen Punkte und das Centroid legen
- Das Kreissegment durch diese drei Punkte bestimmen
- An den Enden dieses Segments 20% frei lassen und in den verbleibenden Teil in der Mitte die Beschriftung setzen
- Die Größe der Schrift hängt von der Textlänge und den Länge des Segments ab. Die Ausrichtung der Schrift hängt von der Lage des Segments ab
Weitere Bastelei wäre nötig, um das auch für Untergruppen (ab Zoom 8 könnte man z.B. "Wetterstein" und "Allgäuer Alpen" durch "Alpen" ersetzen) oder Gebirge mit ungewöhnlichen Formen hinzubekommen.
Ausrichtung von Pässen und Sätteln
Pässe sollten nach Möglichkeit so dargestellt werden, dass die Richtung des Übergangs angedeutet wird. Versuche, diese Richtung anhand der Strasse festzustellen, die darüber führt, sind gescheitert: Oft ist die Zuordnung zur Strasse nicht klar. Kleinere Sättel sind oft nicht nur Übergänge, sondern auch Kreuzungen von "Passüberquerung" und "Kammstrasse". Oft gibt es einen Abstecher, der am Pass zum nahegelegenen Gipfel führt...
Das vielleicht ausichtsreichere Verfahren wurde beim Bild rechts angewendet:
Man
- holt sich die Höhe des Sattels aus den Höhendaten
- geht von dort 100m in alle Richtungen und bestimmt die Höhe dort
- sucht sich die Richtung in der es vorne und hinten am weitesten runter und links/rechts am weitesten rauf geht
- dreht das Icon in diese Richtung
Das Bild wurde mit ASTER-Daten im 90m-Raster berechnet, die einzelnen Messpunkte in 15°-Schritten durch Interpolation dieses Rasters. Gewonnen hat die Richtung mit minimalem Wert für (Vorne-Passhöhe) + (Hinten-Passhöhe) - (Links-Passhöhe) - (Rechts-Passhöhe).
Es hat sich als günstig herausgestellt, auch die Nachbarwerte von Links/Rechts... mit einfliessen zu lassen, hier wurden die direkten Nachbarn mit 1/3 berücksichtigt. Das hilft, auch Sättel darzustellen, die z.B. in Krümmungen eines Gebirgszuges liegen, oder wo der Sattel ein bisschen neben dem Sattelpunkt in der echten Welt (oder wie diese echte Welt laut Höhendaten aussieht), oder wo einfach die verrauschten Höhendaten geglättet werden sollten...
Bis Januar 2017 wurde ein PHP-Progrämmchen zum Rotieren verwendet, inzwischen werden die Daten wie hier nebenan beschrieben aus Aster-Höhendaten im 90m-Raster erzeugt.
Gitternetze
Auf der Seite lassen sich noch ein paar Gitternetze einblenden, um ein Kilometer-Raster in UTM oder Gauß-Krüger-Koordinaten auszugeben. Das ist eigentlich Spielerei, die entsprechenden Passagen im Mapfile lassen sich auch ohne Verlust an Nährwert weglassen. Ich wollte einfach nur mal die Möglichkeiten des Mapservers ausnützen, praktisch beliebige Projektionen auszugeben. Außerdem funktioniert die Maßstabsanzeige des Mapservers nur richtig, wenn man eine Projektion mit echten Metern verwendet. Die Maßeinheit bei Spherical Mercator wird zwar auch mit "m" bezeichnet, ist aber in unseren Breiten deutlich mehr. Das ist auch der Grund, warum die Karten-Exporte in den anderen Projektionen deutlich detailierter erscheinen, die Werte für minscale und maxscale beziehen sich auch auf diese "Meter".
Rendern
Die Icons für Sportplätze, Kirchen... sind von http://www.sjjb.co.uk/mapicons/ Die sind CC-0, was praktisch ist, wenn man die Karte beliebig verwenden möchte. Ein paar sind selbstgemalt. Der Dank für die 2,5D-Darstellung der Häser gebührt User:SunCobalt.
Das Mapfile gibt es hier zum runterladen: User:Maxbe/Kartenversuch/mapfile.
Wers brauche kann, kann sich bedienen. Der Stil steht unter AGPL, wer ihn findet, darf ihn behalten, verändern, verbessern, nutzen und verbreiten. Es sollte aber irgendwo erwähnen, wo er den Stil her hat und seine eigenen Änderungen wieder unter dieser Lizenz zur Verfügung stellen.
Viel Spaß damit!