Workflow pour filtrer les points du sol dans des données Lidar et créer un MNT
Ce workflow est un exemple de cas d’étude que nous avons eu à traiter, il est à adapter en fonction des données utilisées
Nous allons utiliser LibLAS, PDAL et Points2Grid
Données en entrée
Nous avons des nuages de points Lidar au format LAS, avec des coordonnées en WGS84 pour le X et Y, et en mètres pour le Z
Affichage des métadonnées (avec libLAS):
lasinfo input.las
Affichage des métadonnées (avecPDAL):
$ pdal info input.las
il faut vérifier si le système de coordonnées est définit, s’il y a un facteur d’échelle selon X, Y ou Z, et toutes autres informations
Supposons que le système de coordonnées n’est pas définit ou mal définit dans les métadonnées du LAS, cela va causer des problèmes dans la suite du traitement
Définition du système de coordonnées
Créer un fichier JSON pour le traitement PDAL, nommé par exemple “assign_scr.json” :
{ "pipeline": [ { "type" : "readers.las", "filename" : "input.las" }, { "type" : "writers.las", "a_srs": "EPSG:4326", "filename" : "output.las" } ] }
Puis lancer la commande:
pdal pipeline assign_scr.json
Maintenant la définition est correcte, mais il est préférable pour les interpolations et les traitements sur les MNT d’être en projection avec les mêmes unités en horizontal qu’en vertical, nous allons donc reprojeter en UTM
Reprojection
créer un pipeline JSON comme suit, nommé par exemple “reproject.json” :
{ "pipeline":[ { "type":"writers.las", "filename":"input.las" } { "type":"filters.reprojection", "out_srs":"EPSG:32633" }, { "type":"writers.las", "filename":"output.las" } ] }
Puis lancer la commande:
pdal pipeline reproject.json
Il faut ensuite faire une classification sur les points
Classification du sol
créer un pipeline JSON comme suit, nommé par exemple “ground_classif.json” :
{
"pipeline":[
{
"type":"filters.assign",
"assignment":"Classification[:]=0"
},
{
"type":"filters.elm"
},
{
"type":"filters.outlier",
"mean_k":8,
"multiplier":3
},
{
"type":"filters.smrf",
"ignore":"Classification[7:7]",
"slope":0.2,
"window":16,
"threshold":0.45,
"scalar":1.2
},
{
"type":"filters.range",
"limits":"Classification[2:2]"
}
]
}
Ce pipeline est bien décrit dans les aides de PDAL, les étapes en sont les suivantes :
- réassigner tous les points à la classe 0 car on suppose que la classification d’origine n’est pas bonne
- utilisation de l’algo ELM pour assigner les points considérés comme du bruit à la classe 7
- faire un filtrage statistique pour les points abbérants
- utiliser le filtre SMRF sur les points restants pour classer les points
- ne récupérer que les points classés comme ground
Puis lancer la commande:
pdal pipeline ground_classif.json
on récupère donc un fichier las ne contenant que les points classés comme sol (classe 2)
il est maintenant possible d’interpoler les points pour créer un MNT
Génération du MNT
l’interpolation peut se faire avec PDAL ou avec Points2Grid (outil utilisé par OpenTopography pour générer des MNT). L’implémentation de PDAL est basée sur celle de Points2Grid, donc au final c’est équivalent je pense.
La version de points2grid que nous utilisons est celle de développement ici
points2grid -i ground.las -o dtm.tif -r 1 --output_format grid --input_format las --interpolation_mode incore --idw --resolution 1 --fill --fill_window_size 5
pour des explications sur les paramètres taper : points2grid -h
Références:
- https://pdal.io/tutorial/ground-filters.html : filtrage sol sur PDAL, avec d’autres exemples encore sur ce site
- https://opentopography.org/otsoftware/points2grid : présentation générale de points2grid sur le site d’OpenTopography
- https://github.com/CRREL/points2grid/ : version de développement de points2grid
Degrés Minutes Secondes ;