Bezugsrahmenwechsel: Transformation von Rasterdaten #1
15 November 2014
Mit dem NTv2-CHENyx06-Datensatz lassen sich beliebige Datenformate (Vektor und Raster) transformieren. So unterstützt z.B. QGIS seit der Version 2.2 diese Möglichkeit der Datumstransformation. Die Genauigkeit dieser Methode ist im Kanton Solothurn nur marginal schlechter als die strenge Transformation mit dem FINELTRA-Algorithmus. Für die allermeisten Darstellungen am Bildschirm oder für die Transformation von Rasterdaten reicht die Genauigkeit völlig aus.
Ein NTv2-Datensatz kann in verschiedenen WMS-Servern auch dazu verwendet werden Transformationen on-the-fly durchzuführen (z.B. QGIS Server, GeoServer). Will man aber trotzdem die Rasterdaten nicht nur on-the-fly transformieren, sondern sie auch physisch vorliegen haben, kann man das mit einem kleinen Skript und GDAL machen. Oftmals reicht für nicht-hochauflösende Rasterdaten z.B. eine Translation um +2 Mio / +1 Mio. Eine andere Variante ist die strenge Transformation der Metadaten («Worldfile»). Für hochauflösende Rasterdaten funktioniert das nur noch bedingt.
Viele der Rasterdaten liegen in einem nahtlosen Mosaik vor (z.B. Orthofotos). Nach der Transformation sollen weiterhin keine Lücken und Überlappungen zwischen den einzelnen Kacheln sichtbar sein. Um das zu erreichen wird zuerst eine Shapedatei mit den einzelnen Kacheln als Polygon erstellt. Im Kanton Solothurn sind das 1km2-Kacheln:
Dieser sogenannte Tileindex wird mit gdaltindex erstellt:
gdaltindex -write_absolute_path ortho2014.shp *.tifDabei werden sämtliche Tiff-Dateien in dem Verzeichnis berücksichtigt. Mit der Option -write_absolute_path wird im Attribut location in der Shapedatei ortho2014.shp der absolute Pfad der Tiff-Dateien gespeichert.
Jetzt ist es Zeit für das kleine Pythonskript, dass schlussendlich nichts Anderes macht als für jede dieser Kacheln ein gdal_warp-Befehl mit passenden Parameter zusammenzustellen und auszuführen:
#!/usr/bin/python
# -*- coding: utf-8 -*-
from osgeo import ogr, osr
import os
S_SRS = "+proj=somerc +lat_0=46.952405555555555N +lon_0=7.439583333333333E +ellps=bessel +x_0=600000 +y_0=200000 +towgs84=674.374,15.056,405.346 +units=m +units=m +k_0=1 +nadgrids=./chenyx06/chenyx06a.gsb"
T_SRS = "+proj=somerc +lat_0=46.952405555555555N +lon_0=7.439583333333333E +ellps=bessel +x_0=2600000 +y_0=1200000 +towgs84=674.374,15.056,405.346 +units=m +k_0=1 +nadgrids=@null"
ogr.UseExceptions()
shp = ogr.Open("mosaic/ortho2014.shp")
layer = shp.GetLayer(0)
for feature in layer:
infileName = feature.GetField('location')
geom = feature.GetGeometryRef()
env = geom.GetEnvelope()
minX = int(env[0] + 0.001 + 2000000)
minY = int(env[2] + 0.001 + 1000000)
maxX = int(env[1] + 0.001 + 2000000)
maxY = int(env[3] + 0.001 + 1000000)
outfileName = str(minX)[0:4] + str(minY)[0:4] + "_12_5cm.tif"
outfileName = os.path.join(os.path.dirname(infileName), "lv95", outfileName)
cmd = "/usr/local/gdal/gdal-dev/bin/gdalwarp -s_srs \"" + S_SRS + "\" -t_srs \"" + T_SRS + "\" -te " + str(minX) + " " + str(minY) + " " + str(maxX) + " " + str(maxY)
cmd += " -tr 0.125 0.125 -wo NUM_THREADS=ALL_CPUS -co 'PHOTOMETRIC=RGB' -co 'TILED=YES' -co 'PROFILE=GeoTIFF'"
cmd += " -co 'INTERLEAVE=PIXEL' -co 'COMPRESS=DEFLATE' -co 'BLOCKXSIZE=512' -co 'BLOCKYSIZE=512'"
vrt = "/opt/Geodaten/ch/so/kva/orthofoto/2014/rgb/12_5cm/ortho2014rgb.vrt"
cmd += " -r bilinear " + vrt + " " + outfileName
os.system(cmd)
cmd = "/usr/local/gdal/gdal-dev/bin/gdal_edit.py -a_srs EPSG:2056 " + outfileName
os.system(cmd)
cmd = "/usr/local/gdal/gdal-dev/bin/gdaladdo -r nearest --config COMPRESS_OVERVIEW DEFLATE --config GDAL_TIFF_OVR_BLOCKSIZE 512 " + outfileName + " 2 4 8 16 32 64 128"
os.system(cmd)Zeilen 4 - 5: Benötigte Module werden geladen.
Zeilen 7 - 8: Es ist möglich neben EPSG-Codes für Transformation mit GDAL auch eigene proj4-Definitionen als Parameter zur übergeben. Bei der Verwendung von NTv2-Gittern ist das ein bisschen hakelig. Es gibt viele Varianten. Viele davon führen zu falschen Resultaten, zwei davon führen zu richtigen Resultaten. Abhängig davon welcher NTv2-Datensatz (chenyx06a.gsb oder chenyx06etrs.gsb) verwendet wird. Langer Rede kurzer Sinn: Das hier ist eine Kombination, die korrekt transformiert.
Zeile 12 - 13: Die Shapedatei mit der Kacheleinteilung wird geöffnet und ein Layer wird angelegt.
Zeilen 15 -16: Die For-Schleife bearbeitet jede Kachel einzeln. Zuerst wird das Attribut location ausgelesen. Es beinhaltet den kompletten absoluten Pfad der einzelnen Tiff-Datei. Davon verwenden wir aber später nur das Verzeichnis, da wir in ein Unterverzeichnis die resultierenden Tiff-Dateien speichern wollen.
Zeilen 17 - 18: Die Geometrie der Kachel wird ausgelesen und die Envelope (Tupel mit den min/max Koordinatenwerten) berechnet.
Zeilen 20 - 23: Aus der Envelope lesen wir die minimalen X- resp. Y-Werte. Die resultierenden Kacheln sollen wie die Ausgangskacheln «schöne» Kilometerkacheln sein. Dazu wird einfach 2 Mio. resp. 1 Mio dazugerechnet. Weil der int-Befehl abrundet (?), muss den Koordinatenwerten ein kleiner Wert dazu addiert werden, um keine falschen Ganzzahleswerte zu bekommen.
Zeilen 25 - 26: Aus den Koordinatenwerten wird der Dateinnamen der neuen Kacheln bestimmt.
Zeilen 28 - 33: Hier geschieht die Zauberei. Als Ausgangsdatei wird nicht die einzelne Kacheln verwendet sondern immer die vrt-Datei (Diese muss u.U. vorgängig noch erstellt werden). GDAL sieht jetzt nur eine einzelne, flächendeckende Rasterdatei. Somit muss nur noch der Ausschnitt gewählt werden, der aus dieser einzelnen Rasterdatei ausgeschnitten werden soll. Dieser Ausschnitt wurde ja bereits in den Zeilen 20 - 23 ermittelt. Das Resultat wird wiederum in einer Tiff-Datei gespeichert.
Zeilen 35 - 39: Anschliessend werden in der neuen GeoTiff-Datei noch korrekte Metadaten zum Koordinatensystem gespeichert und es werden interne Overview gerechnet.
Das Resultat sind einzelne Orthofotokacheln, die sich weder überlappen noch sind Lücken zwischen den Kacheln vorhanden. Die Übergänge sind wie im Ausgangsmaterial nicht sichtbar. Einzig am Perimeterrand sind kleine schwarze Ränder sichtbar. Diese entstehen weil im Ausgangsmaterial in diesem Bereich keine Daten vorhanden sind. Liegt das Ausgangsmaterial im Bezugsrahmen LV95 vor und wird in den Bezugsrahmen LV03 transformiert, dürfte dieses Phänomen nicht auftreten.
Der Transformationsprozess dauerte für 384 Kacheln circa 1h 45m. Ein Qualitätsverlust aufgrund des Resamplings ist nicht feststellbar.