Die bekannte Mercator-Projektion, allgegenwärtig bei Karten im Internet, hat eine wichtige Eigenschaft – sie weist über das Gebiet, in dem sie verwendet wird, einen stark variablen Maßstab auf. Das ist der Grund, weshalb zum Beispiel Grönland auf diesen Karten größer erscheint als Südamerika. Der Maßstab variiert über den Kartenbereich des üblichen quadratischen Ausschnitts der Internet-Karten um mehr als den Faktor 11 und selbst wenn man den äußersten Rand im Norden und Süden in der zentralen Antarktis und im Arktischen Ozean außen vor lässt ergibt sich immer noch ein Maßstabsverhältnis von etwa 1:8 zwischen Äquator und hochpolaren Breiten.
Im Prinzip ist dies kein großes Problem für Internet-Karten, denn man kann ja frei rein- und rauszoomen. Aus diesem Grund hat man bei der ursprünglichen Wahl der Projektion für solche Karten diese Eigenschaft der ansonsten enormen Vorteile wegen bereitwillig in Kauf genommen. Das Problem ist jedoch, dass traditionell Kartengestaltung und Geodaten-Verarbeitung nicht gewohnt sind, mit so großen Maßstabsunterschieden umzugehen. In anderen Worten ausgedrückt: Während im Prinzip fast jeder über diese Eigenschaft der Projektion Bescheid weiß, ignorieren fast alle diese Tatsache in der Praxis und jeder tut so, als ob die Mercator-Karte sich nicht von Darstellungen in anderen Projektionen unterscheidet.
Eines der offensichtlichsten Beispiele hierfür sind die diversen Karten zur Daten-Statistik, insbesondere die bekannten Darstellungen der Datendichte wie Martin Raifer’s Karte der OpenStreetMap Knoten-Dichte – Diese Karten sind schlicht und einfach falsch, wenn man sie als quantitative Darstellungen der Datendichte interpretiert.
Solche Karten zeigen eigentlich die Anzahl der Elemente pro Quadratmeter in projizierten Koordinaten, was etwas sehr viel anderes ist als die Dichte bezogen auf die tatsächliche Fläche auf der Planeten-Oberfläche. Und dieser Unterschied ist enorm – wenn der lineare Maßstab um den Faktor 11.5 zwischen Äquator und Kartenrand variiert, unterscheidet sich die Dichte um das Quadrat dieses Faktors, also mehr als 130.
Ich habe ein kleines Werkzeug auf github zur Verfügung gestellt, mit dem man diesen Unterschied in der Nachbearbeitung korrigieren kann. Das Werkzeug ist recht generisch, es sollte mit jedem GDAL-kompatiblen Raster-Format funktionieren für welches GDAL Lese- und Schreibunterstützung bietet und wo Koordinatensystem-Informationen enthalten sind sowie für jede Projektion – solange proj4 die inverse Transformation hierfür kennt, was auch für gdalwarp erforderlich ist.
Ich demonstriere das hier für die Punktdichte eines OSM-Küstenlinien-Extraktes – bin nicht geduldig genug, das über den ganzen planet laufen zu lassen. Um Punktdichten zu zählen, kann man das neue node_density aus libosmium verwenden, welches Jochen zur Verfügung gestellt hat oder mein eigenes gdal_nodedensity, welches auch die Darstellung von Liniendichten unabhängig von der Punktdichte der Linien selbst erlaubt. Durch den Aufruf von gdal_nodedensity
gdal_nodedensity -a_srs EPSG:3857 -ot Float32 -ts 1024 1024 -te -20037508.342789244 -20037508.342789244 20037508.342789244 20037508.342789244 osm_coastlines_tmp.osm coast_nodedensity.tif
erzeugt man die folgende falsche Darstellung der Punktdichte, wo die Dichte in Richtung der Pole zu niedrig angezeigt wird.
Führt man nun folgendes aus:
gdal_valscale coast_nodedensity.tif
erhält man eine korrekte Darstellung der Dichten:
Das ist alles – keine weiteren Parameter sind nötig, alle Informationen werden der Bilddatei entnommen.
Im gdal-tools repository findet sich noch ein weiteres Programm, welches eine Raster-Maske um eine wählbare Distanz puffert – wiederum in tatsächlichen Längeneinheiten und nicht in projizierten Koordinaten. Mit der Küstenlinien-Maske als Beispiel (erzeugt aus den Küstenlinien-Shapefiles mittels gdal_rasterize --config GDAL_CACHEMAX 2048 -a_srs EPSG:3857 -ot Byte -ts 1024 1024 -te -20037508.342789244 -20037508.342789244 20037508.342789244 20037508.342789244 -burn 255 -at osm_coastlines_land_3857.shp coast_mask.tif – das -at stellt sicher, dass keine kleinen Inseln verloren gehen):
kann man aufrufen:
gdal_maskbuffer coast_mask.tif 200000
und erhält eine um 200 km expandierte Küstenlinie:
Keines der üblichen Pakete zur Vektordaten-Verarbeitung bietet Puffer-Funktionen mit variablem Abstand. Mittels ST_Buffer() in PostGIS erhält man zum Beispiel einen konstanten Abstand im gewählten Koordinatensystem, auch bei Verwendung des Geography-Datentyps, und wenn sich das Element über den halben Globus erstreckt ist das Ergebnis falsch (siehe auch die Warnung in der verlinkten Dokumentation). Der einzige Weg, in solchen Fällen bei einem globalen Datensatz eine korrekte Pufferung zu erhalten, wäre die Aufteilung aller Geometrien in genügend keine Stücke, so dass die Maßstabsunterschiede pro Stück nicht mehr ins Gewicht fallen, und dann die einzelnen Elemente mit einem jeweils individuell bestimmten Radius zu puffern. Das wäre eine Menge Arbeit und würde auch wieder schief gehen, sobald man mit einem sehr großen Radius puffern möchte, so dass die Maßstabsunterschiede dabei nicht vernachlässigbar sind. Man müsste dann die Pufferung in Schritten durchführen und zwischendrin die Geometrien erneut aufteilen…
gdal_maskbuffer unterstützt im Moment nur isotrope Maßstabsunterschiede – proj4 würde auch Anisotropien unterstützen, die Abstandsfunktion in CImg deckt aber nur den isotropen Fall ab. Auch die periodischen Randbedingungen in x-Richtung sind nicht unterstützt. Beide Werkzeuge habe ich nur mit der Mercator-Projektion (EPSG:3857) getestet, sie sollten jedoch auch für andere Projektionen funktionieren. gdal_valscale lädt im Moment auch noch unnötigerweise das gesamte Bild auf einmal in den Speicher, was bei großen Bildern zu Problemen führen könnte (Wenn ich mich recht erinnere ist GDALDataset::RasterIO() nicht in der Lage, Blöcke größer als das 32bit-Limit zu bearbeiten). Intern verwenden beide Werkzeuge die pj_factors()-Funktion aus proj4, um die benötigten Skalierungsfaktoren zu ermitteln, welche man für einzelne Koordinatenpunkte auch vom proj-Werkzeug ausgeben lassen kann.