PDAL Workshop

Beispiel: Schutzwald Katzberg

Lage

Karte

Information

Bitte starte die Conda-Konsole und aktiviere die PDAL Umgebung bevor du mit dem Beispiel beginnst.

conda activate pdal

Wechsle anschließend mit 'cd' in das WS-Verzeichnis.

1. Merge

Vereinigen aller Dateien mit der Endung 2014.laz im WS-Verzeichnis

Pipeline e01_merge.json
[
    "las_a_2014.laz",
    "las_b_2014.laz",
    "las_c_2014.laz",
    "las_d_2014.laz",    
    {
        "type": "filters.merge"
    },
    "merge2014.laz"
]
Aufruf
pdal pipeline e01_merge.json

2. Crop

Beschneiden der Datei merge2014.laz auf ein Rechteck: PLL(640007;5583845), PUR(640087;5583934)

Pipeline e02_crop.json
[
    "merge2014.laz",
    {
      "type": "filters.crop",
      "polygon": "Polygon ((640007 5583845,640087 5583845,640087 5583934,640007 5583934, 640007 5583845))"
    },
    "crop2014.laz"
  ]
Aufruf
pdal pipeline e02_crop.json

3. 3D-Map View

Betrachte die Datei crop2014.laz in QGIS

  1. QGIS starten

  2. Projekt als 'pdal-ws' im WS-Ordner speichern

  3. Datei crop2014.laz zur Karte hinzufügen

  4. Layer Properties aufrufen

  5. Tab Symbology aufrufen

  6. Dialog zeigt die Klassifikation nach Classification Attribut

    Details

  7. Tab: Stats aufrufen

  8. Dialog zeigt die Attribute und Classification Statistik

    Details

  9. Tab: Elevation öffnen

  10. Damit in der 3D Ansicht auf die Scene gezoomt werden kann, ändern wir den Elevation Offset von 0 auf -450m

    Details

  11. Neue 3D Kartenansicht anlegen mit Menu 'View/3D Map Views/New 3D Map View'

  12. 3D Map View zeigt Layer

  13. Configuration öffnen

  14. Shadows aktivieren

    Details

  15. Navigation: Zoom in/out über Maus-Wheel-Rad

  16. 3D Map view zeigt LAS mit Schatten

    Details

Information

Um das Handling von großen Point Cloud Daten zu verbessern, erzeugt QGIS beim erstmaligen Lesen einer LAS/LAZ Datei eine Cloud Optimized Point Cloud Dateiopen in new window mit der Endung .copc.laz im Verzeichnis der Originaldatei. Falls sich die Originaldatei ändert, muss der Cache manuell gelöscht werden.

Information

Weitere Möglichkeiten zur Darstellung von 3D-Daten werden im QGIS-Handbuch erläutert:

4. Punktdichtekontrolle

Erzeugen eines Hex-Bin Vektor-Layers zur Dichtekontrolle mit Hilfe der Density Anwendung

a. Rufe die Hilfe zur pdal density Anwendung auf

Aufruf
pdal density --help

b. Erzeuge einen GeoJSON Hex-Bin Density Layer mit 5m Kantenlänge

Aufruf
pdal density --lyr_name density --ogrdriver GeoJSON --edge_length 5 crop2014.laz crop2014hex5.geojson

c. Öffne QGIS und erzeuge eine 'Graduated Symbology Legend' auf dem Feld 'crop2014hex5.count' mit der Klassifizierungsmethode 'Natural Breaks (Jenks)' und sieben Werteklassen

QGIS Karte

5. Eliminieren von Ausreißern

Löschen von Punkten mit einer Höhe < 450m

Pipeline e05_outlier.json
[
    "crop2014.laz",
    {
      "type": "filters.expression",
      "expression": "(Z>=450)"
    },
    "crop2014g.laz"
]
Aufruf
pdal pipeline e05_outlier.json

Notiz

Eine weitere Möglichkeit zur Detektion und Klassifizierung von Punkten mit niedrigen Höhen ist die Anwendung des ELM Filters nach Chen.
filters.elmopen in new window

6. Categorized Symbology

Lade die Datei classification.geojson in QGIS und erzeuge eine 'Categorized Symbology' auf dem Feld 'name'

QGIS Karte

7. Vector Overlay mit Filter

Pipeline
a. Overlay von crop2014g.laz mit dem Polygon aus der Datei classification.geojson
b. Extraktion aller Punkte innerhalb der Schutzwälder: (Classification==33)

Pipeline e07_overlay_expression.json
[
  "crop2014g.laz",
  {
    "type": "filters.overlay",
    "dimension": "UserData",
    "datasource": "classification.geojson",
    "layer": "classification",
    "column": "name"
  },
  {
    "type": "filters.expression",
    "expression": "UserData==33"
  },
  "overlay2014.laz"
]
Aufruf
pdal pipeline e07_overlay_expression.json
QGIS Karte

8. Ground Klassifikation

Pipeline

Pipeline e08_ground.json
[
    "crop2014g.laz",
    {
      "type": "filters.assign",
      "value": "Classification = 0"
    },
    {
       "type":"filters.smrf",     
       "cell": 0.5
    },
    {
       "type":"filters.expression",
       "expression":"Classification==2"
    },
    {
    "type":"writers.las",
    "filename":"ground2014g.laz"
    }    
 ]
Aufruf
   pdal pipeline e08_ground.json
QGIS Karte

9. DGM

9.1 Direkter Export

a. GDAL-Export

Aufruf
pdal translate ground2014g.laz ground2014g.tif --writers.gdal.resolution=0.5

b. Anzeige in QGIS

Karte mit NODATA Lücken

c. Füllen der NODATA Werte über Interpolation in QGIS

Dialog Raster/Analysis/Fill nodata

d. Kontrolle in QGIS

Details

Karte ohne NODATA Lücken

3D Map View

9.2 Triangulation

Pipeline aus fiters.delaunayopen in new window, filters.facerasteropen in new window und writers.rasteropen in new window

Pipeline e09_tin.json
[
    "ground2014g.laz",
    {
        "type": "filters.delaunay"
    },
    {
        "type": "filters.faceraster",
        "resolution": 0.5,
        "width": 160,
        "height": 160,
        "origin_x": 640007,
        "origin_y": 5583845
    },    
    {
        "type": "writers.raster",
        "filename":"ground2014t.tif"
    }
]
Aufruf
pdal pipeline e09_tin.json
QGIS Karte

10. Baum Klassifikation und Analyse

Klassifikation der Bäume in einer Pipeline aus filters.hag_delaunayopen in new window, filters.sortopen in new window, filters.litreeopen in new window,filters.expressionopen in new window und writers.lasopen in new window

Pipeline e10_tree.json
[
    "crop2014g.laz",
    {
        "type": "filters.hag_delaunay"
    },
    {
        "type": "filters.sort",
        "dimension": "HeightAboveGround",
        "order": "DESC"
    },
    {
        "type": "filters.litree",
        "min_points": 50,
        "min_height": 10.0,
        "radius": 20.0
    },
    {
        "type": "filters.expression",
        "expression": "(Classification!=2) && (Classification!=30) && (ClusterID > 0)"
    },
    {
        "type": "writers.las",
        "filename": "tree2014.laz",
        "extra_dims": "ClusterID=uint8,HeightAboveGround=float32",
        "minor_version": 4
    }
]
Aufruf
# Achtung: Die Berechnung dauert ca. 2 min. bis ca. 150 Bäume erkannt werden ....
pdal -v 8 pipeline e10_tree.json
QGIS

2D Karte
2D Karte
3D View
3D View

Pipeline zur Berechnung der maximalen Höhe über Grund für jeden Baum. Die Pipeline verwendet filters.pythonopen in new window zur Ermittlung der Höhe und schreibt das Ergebnis über writers.ogropen in new window in eine GeoJSON-Datei.

Pipeline e10_tree_stats.json
[
    "tree2014.laz",
    {
        "type": "filters.python",
        "script": "lib.py",
        "function": "max_z_cluster",
        "module": "lib"
    },
    {
        "type": "writers.ogr",
        "filename": "tree2014.geojson",        
        "attr_dims": "ClusterID, HeightAboveGround"
    }
]
Python Datei `lib.py`
import numpy as np

def max_z_cluster(ins, outs):           
    max = {}
    # Index des Punktes mit max z je Cluster
    it = np.nditer(ins["Z"], flags=['f_index'])
    for z in it:    
        c = ins['ClusterID'][it.index]
        if c not in max:
            max[c] = it.index
        else:            
            if ins['Z'][max[c]] < z:
                max[c] = it.index

    # Ausgabe über Maske
    outs['Mask'] =  np.zeros(len(ins["Z"]), dtype = bool)        
    for index in max.values():
        outs["Mask"][index] = True    
                
    return True
Aufruf
pdal pipeline e10_tree_stats.json
QGIS

2D Karte mit Stats zur Höhe

11. DOP als RGB Farbwert

Zuweisen von Farbwerten aus einer Raster-Datei mit filters.colorizationopen in new window

Pipeline e11_dop.json
[
    "merge2014.laz",
    
    {
        "type":"filters.colorization",
        "raster":"dop2022.tif"
    },    
    {
        "type":"filters.expression",
        "expression":"(RED>0) && (Classification!=30)"
    },
    {
        "type":"writers.copc",
        "filename":"dop2014.copc.laz"
    }
]
Aufruf
pdal pipeline e11_dop.json
QGIS 3D View

3D View
3D View

12. Schutzwalddarstellung

Skalieren von Farbwerten mit Hilfe von filters.assignopen in new window

Pipeline e12_dop_sample.json
[
    "dop2014.copc.laz",
    {
        "type":"filters.overlay",
        "dimension":"Classification",
        "datasource":"classification.geojson",
        "layer":"classification",
        "column":"name"
    },    
    {
        "type":"filters.assign",
        "value": [
            "Green = Green * 1.5 WHERE Classification == 33",            
            "Blue = Blue * 1.5 WHERE Classification == 33"            
        ]
    },
    {
        "type":"filters.assign",
        "value": [
            "Green = 255 WHERE Green > 255",            
            "Blue = 255 WHERE Blue > 255"           
        ]
    },
    "dop2014s.copc.laz"
    
]
Aufruf
pdal pipeline e12_dop_sample.json
3D View

3D View

Last Updated: