Bezugsrahmenwechsel: Transformation von Lidardaten
24 June 2015
Im Frühjahr 2014 liess der Kanton Solothurn eine LiDAR-Befliegung über das gesamte Kantonsgebiet durchführen. Die Daten wurden, wie von uns gewünscht, im Bezugsrahmen LV03 geliefert. Nicht die sinnvollste Entscheidung… Der Bezugsrahmenwechsel steht definitiv vor der Tür und falls wir wollen, dass die Rohdaten in Zukunft noch verwendet werden, müssen wir wohl oder übel auch die 1'066 LAZ-Dateien transformieren.
Für die Transformation stehen - je nach Genauigkeit der Daten - verschiedene Transformationsmethoden zur Verfügung. Eine davon ist die Transformation mittels NTv2-Gitter. Anstelle der Dreiecke wird hier ein regelmässiges Gitter verwendet. Die Differenz zwischen der Transformation mittels Dreiecksvermaschung und NTv2-Gitter ist im Kanton Solothurn sehr klein:
Die Abbildung zeigt die Differenzen eines 100 m x 100 m Rasters. In drei sehr kleinen Gebieten wird die maximale Abweichung von 25 mm erreicht. Der Mittelwert ist 1.3 mm, der Median bloss 0.6 mm. Die Genauigkeit der LiDAR-Daten liegt im einstelligen bis tiefen zweistelligen Zentimeterbereich. Für die Transformation genügt das NTv2-Gitter also vollkommen.
Eine Bedingung an die transformierten Daten ist, dass sie ebenfalls wieder in Quadratkilometer-Kacheln mit «schönen» Boundingbox-Koordinaten vorliegen:
Das Bild zeigt exemplarisch vier Kacheln. Im linken Teil sind diese vier Kacheln im Bezugsrahmen LV03 gezeichnet. Im rechten Teil (Bezugsrahmen LV95) sind die transformierten Kacheln grau gezeichnet. Diese «hässlichen» Koordinaten sind aber für die Kacheleinteilung nicht erwünscht, sondern die «schönen» runden Koordinaten. Das bedeutet, dass die graue Kachel mit dem Gebäude neu in vier Kacheln zu liegen kommt. Anders formuliert: Das Gebäude liegt in einer neuen Kachel, die Daten aus vier alten Kacheln beinhaltet.
Die Software muss also etwas Ähnliches wie GDAL mit seinem vrt-Treiber anbieten. Dann kann man genau gleich vorgehen wie bereits mit den Orthofotos.
PDAL ist ein relativ neues Projekt. Es ist eine Bibliothek für "translating and manipulating point cloud data". Was GDAL für 2D-Pixel-Daten heute ist, will PDAL für multidimensionale Punkte werden.
Zuerst muss ein Tileindex für sämtliche Kacheln erzeugt werden:
find /home/stefan/mr_candie_nas/Geodaten/ch/so/agi/hoehen/2014/lidar/ -maxdepth 1 -iname "*.laz" | pdal tindex tileindex.gpkg -f "GPKG" --a_srs EPSG:21781 --t_srs EPSG:21781 --fast_boundary --stdin --lyr_name tileindexDer Befehl sucht alle Dateien mit der Endung *.laz und schreibt die Boundingbox in eine GeoPackage-Datei. Mit der Option --a_srs EPSG:21781 wird den LiDAR-Daten ein Koordinatensystem zugewiesen. Ob die LiDAR-Daten bereits mit einem Koordinatensystem versehen sind, lässt sich mit dem pdal info-Befehl herausfinden:
pdal info --metadata LAS_592228.lazliefert:
{
"filename": "LAS_592228.laz",
"metadata":
{
"comp_spatialreference": "",
"compressed": true,
"count": 14347366,
"creation_doy": 312,
"creation_year": 2014,
"dataformat_id": 1,
"dataoffset": 329,
"filesource_id": 0,
"global_encoding": "AAA=",
"header_size": 227,
"major_version": 1,
"maxx": 592999.98999999999,
"maxy": 228999.98999999999,
"maxz": 2134.52,
"minor_version": 2,
"minx": 592000,
"miny": 228000,
"minz": 918.08000000000004,
"offset_x": -0,
"offset_y": -0,
"offset_z": -0,
"project_id": "00000000-0000-0000-0000-000000000000",
"scale_x": 0.01,
"scale_y": 0.01,
"scale_z": 0.01,
"software_id": "TerraScan",
"spatialreference": "",
"system_id": ""
},
"pdal_version": "1.0.0.b1 (git-version: 1dce22)"
}Hat die Datei bereits ein Koordinatensystem, sieht der Output so aus:
{
"filename": "LAS_592228.laz",
"metadata":
{
"comp_spatialreference": "PROJCS[\"CH1903 / LV03\",GEOGCS[\"CH1903\",DATUM[\"CH1903\",SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],TOWGS84[674.4,15.1,405.3,0,0,0,0],AUTHORITY[\"EPSG\",\"6149\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4149\"]],PROJECTION[\"Hotine_Oblique_Mercator_Azimuth_Center\"],PARAMETER[\"latitude_of_center\",46.95240555555556],PARAMETER[\"longitude_of_center\",7.439583333333333],PARAMETER[\"azimuth\",90],PARAMETER[\"rectified_grid_angle\",90],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",600000],PARAMETER[\"false_northing\",200000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Y\",EAST],AXIS[\"X\",NORTH],AUTHORITY[\"EPSG\",\"21781\"]]",
"compressed": true,
"count": 14347366,
"creation_doy": 224,
"creation_year": 2015,
"dataformat_id": 3,
"dataoffset": 2181,
"filesource_id": 0,
"global_encoding": "AAA=",
"header_size": 227,
"major_version": 1,
"maxx": 592999.98999999999,
"maxy": 228999.98999999999,
"maxz": 2134.52,
"minor_version": 2,
"minx": 592000,
"miny": 228000,
"minz": 918.08000000000004,
"offset_x": 0,
"offset_y": 0,
"offset_z": 0,
"project_id": "00000000-0000-0000-0000-000000000000",
"scale_x": 0.01,
"scale_y": 0.01,
"scale_z": 0.01,
"software_id": "PDAL 1.0.0.b1 (e412bd)",
"spatialreference": "PROJCS[\"CH1903 / LV03\",GEOGCS[\"CH1903\",DATUM[\"CH1903\",SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],TOWGS84[674.4,15.1,405.3,0,0,0,0],AUTHORITY[\"EPSG\",\"6149\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4149\"]],PROJECTION[\"Hotine_Oblique_Mercator_Azimuth_Center\"],PARAMETER[\"latitude_of_center\",46.95240555555556],PARAMETER[\"longitude_of_center\",7.439583333333333],PARAMETER[\"azimuth\",90],PARAMETER[\"rectified_grid_angle\",90],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",600000],PARAMETER[\"false_northing\",200000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Y\",EAST],AXIS[\"X\",NORTH],AUTHORITY[\"EPSG\",\"21781\"]]",
"system_id": "PDAL"
},
"pdal_version": "1.0.0.b1 (git-version: 1dce22)"
}Die Option --t_srs beschreibt das Koordinatensystem des zu erzeugenden Tileindexes. Mit --fast_boundary wird die Boundingbox nicht aus den Daten selbst eruiert, sondern es werden die Koordinaten aus den Metadaten verwenden. Und ja, es ist massiv schneller.
Anschliessend an die Erstellung des Tileindexes, kann mittels Pythonskript, einer For-Schleife und dreier PDAL-Befehlen das gewünschte Resultat berechnet werden:
#!/usr/bin/python
# -*- coding: utf-8 -*-
import os.path
from osgeo import ogr, osr
import os
import sys
ogr.UseExceptions()
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"
LAYERNAME = "tileindex"
TEMPDIR = "/tmp/"
TILEINDEX = "/home/stefan/tmp/lidar/tileindex.gpkg"
OUTDIR = "/home/stefan/tmp/"
BUFFER = 2
SCALE = 0.01
gpkg = ogr.Open(TILEINDEX)
lyr = gpkg.GetLayerByName(LAYERNAME)
for feat in lyr:
cmd = "rm " + os.path.join(TEMPDIR, "*.las")
os.system(cmd)
cmd = "rm " + os.path.join(TEMPDIR, "*.laz")
os.system(cmd)
filename = feat.GetField("location")
print "*** " + os.path.basename(filename) + " ***"
env = feat.GetGeometryRef().GetEnvelope()
min_x = int(env[0] + 0.001)
min_y = int(env[2] + 0.001)
max_x = int(env[1] + 0.001)
max_y = int(env[3] + 0.001)
# Create new file name.
filename_lv95 = "LAS_" + str(min_x/1000 + 2000) + "_" + str(min_y/1000 + 1000) + ".laz"
filename_lv95 = os.path.join(OUTDIR, filename_lv95)
# Get a slightly larger tile. Two meters buffer is enough. We crop afterwards.
# Reason: 620'000 -> 2'620'000.65
bounds = "(["+str(min_x-BUFFER)+", "+str(max_x+BUFFER)+"], ["+str(min_y-BUFFER)+", "+str(max_y+BUFFER)+"])"
cmd = 'pdal tindex --merge ' + TILEINDEX + ' --lyr_name tileindex --bounds "' + str(bounds) + '" --t_srs EPSG:21781 ' + os.path.join(TEMPDIR, 'tmp.las')
os.system(cmd)
# Now change reference frame with ntv2.
cmd = 'pdal translate --a_srs "'+S_SRS+'" --t_srs "'+T_SRS+'" -i ' + os.path.join(TEMPDIR, 'tmp.las') + ' -o ' + os.path.join(TEMPDIR, 'tmp_lv95.las')
os.system(cmd)
# Set better/nicer EPSG:2056. The ntv2 transformation lacks of datum name.
# Crop to nice bbbox.
bounds = "(["+str(min_x+2000000)+", "+str(max_x+2000000)+"], ["+str(min_y+1000000)+", "+str(max_y+1000000)+"])"
cmd = 'pdal translate --a_srs EPSG:2056 --t_srs EPSG:2056 --writers.las.format="1" --bounds "' + bounds + '" -z -i ' + os.path.join(TEMPDIR, 'tmp_lv95.las') + ' -o ' + os.path.join(OUTDIR, filename_lv95)
os.system(cmd)Das Grundprinzip resp. -vorgehen steht bereits hier. Interessant sind bloss die drei PDAL-Aufrufe:
Zeilen 48 - 50: Es wird aus dem Tileindex eine etwas (BUFFER = 2 Meter) grössere Kachel ausgeschnitten. Es muss nur soviel mehr sein, dass wir anschliessend an die Transformation garantiert «schöne» LV95-Koordinaten ausschneiden können. Mit pdal tindex und der Option --merge können Daten, die in einem Tileindex vorliegen, zusammengefügt werden. Zusätzlich gibt es die Option --bounds mit der nur ein ganz bestimmtes Rechteck gemerged resp. ausgeschnitten werden kann. Leider - so wie ich es verstanden habe - werden jeweils zuerst sämtliche vom Rechteck betroffene Kacheln komplett zusammengefügt und erst anschliessend ausgeschnitten. Das schlägt sich einerseits in der Geschwindigkeit und andererseits im RAM-Verbrauch nieder. In Normalfall sind ja neun Kacheln betroffen, die zuerst zusammengefügt werden müssen. Das braucht circa 9 GB RAM.
Zeilen 53 - 54: Nach dem Zusammenfügen und Ausschneiden können die Daten mit dem NTv2-Grid transformiert werden. Dazu verwendet man den Befehl pdal translate.
Zeilen 58 - 59: Zu guter Letzt müssen die transformierten Daten auf die «schönen» Koordinaten zurechtgeschnitten werden. Zudem kann in diesem Schritt das Koordinatensystem in den Metadaten sauber gesetzt werden (--a_srs EPSG:2056 und --t_srs EPSG:2056). Sonst stehen da die hässlichen PROJ4-Strings der NTv2-Transformation drin. Für das Zuschneiden kann der Befehl pdal translate mit der Option --bounds verwendet werden.
Das Resultat sind jetzt 1'066 LiDAR-Kacheln mit «schönen» Boundingbox-Koordinaten in LV95.
Die ganze Rechnerei dauert aber ein wenig. Liegen die Daten lokal auf der SSD vor, dauert das Prozedere für eine Kachel circa 150 Sekunden. Mein Setup war aufgrund des Platzbedarfs leider ziemlich übel: Die Daten liegen auf einem NAS, das via WLAN an den Computer angebunden ist. Gerechnet wird in einer virtuellen Maschine, in der das NAS mit sshfs gemounted ist. So rechnete das fast 4 Tage vor sich hin. Also fast doppelt so lange.
Hätte man jedoch genügend RAM, könnte man den Prozess leicht parallelisieren. Man muss bloss das Skript x-mal anstossen und die For-Schleife dürfte nicht den gleichen Tileindex als Input haben.
Der Prozess selbst konnte leicht optimert werden, indem beim Zwischenschritt nicht mit der Option -c komprimiert wurde. Ebenfalls gewinnbringend wirkt sich der Einsatz von --bounds anstelle --polygon aus. Will man also bloss ein Rechteck ausscheiden und nicht ein beliebiges Polygon reicht --bounds allemal und pdal muss in diesem Fall keine teuren GEOS-Operationen verwenden.