LiDAR HD brut - Distinguer sol et sursol avec pdal#

Date de publication initiale : 24 juin 2024

Prérequis#

• PDAL
• bash
• GDAL/OGR

Intro#

Le relevé LiDAR (Light Detection and Ranging) est une technique de plus en plus utilisée pour la création de modèles numériques de précision. Cependant, les données LiDAR de part leur volumétrie et leurs spécificités, peuvent être difficiles à manipuler et à interpréter, en particulier lorsqu'il s'agit de distinguer le sol et le sursol.

Heureusement, des outils tels que PDAL (Point Data Abstraction Library) et GDAL (Geospatial Data Abstraction Library) aident à traiter et à analyser les données LiDAR pour créer des MNT ou des MNS.
Dans cet article, je vais vous expliquer comment j'ai utilisé ces deux outils pour distinguer le sol et le sursol à partir des données LiDAR brutes de l'IGN, et fournir des exemples de code pour vous aider à démarrer dans la manipulation de ces données.

Commenter cet article

Processus global#

Voici un schéma récapitulatif de la procédure mise en place.

graph TD
    A[LiDAR HD] -->G(Décompression)
    G --> C(Colorisation)
    B[Ortho] --> C
    C --> R(Filtrer le bâti)
    R --> D(Sol)
    R --> E(Sursol)
    D --> F(Raster)

Un environnement de travail : config.env#

Avant de se lancer, il est bon de vous parler du fichier de configuration que vous devrez adapter à votre organisation et qui sera utilisé par la suite pour télécharger et traiter les images. On y définit le répertoire de travail et différentes variables nécessaires à la bonne exécution des scripts.

Voici le fichier config.env à adapter :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# REPERTOIRE DE TRAVAIL
REPER='/home/....'

# REPERTOIRE DE STOCKAGE DES LOGS
REPER_LOGS='logs'

# REPERTOIRE D ENTREE
REPER_IN='data_in'

# REPERTOIRE TEMPORAIRE
REPER_TMP='data_tmp'

# REPERTOIRE DE SORTIE
REPER_OUT='data_out'

# PARAMETRES OGR
ENCODAGE='UTF-8'

Consulter le fichier de configuration

Un script pour orchestrer les différentes étapes#

Maintenant, entrons dans le vif du sujet avec la présentation des différentes étapes qui me permettent de traiter les données LiDAR brutes et de distinguer le sol et le sursol. Je vais décrire ici chaque étape et fournir des exemples de code pour vous aider à comprendre le processus de traitement global.

Consulter le fichier "maître"

La décompression#

Les fichiers bruts ayant été livrés au format 7-zip, la première étape consiste à extraire les 4 dalles contenues dans chaque fichier 7z.

# Installation de p7zip : sudo apt install p7zip
7z x "$REPER/data_in/$base.7z" -o$REPER'/data_tmp/un7z'

Colorisation des images#

Dans cette étape, nous allons utiliser pdal pour affecter les informations colorimétriques de l'image aérienne sur chaque point LiDAR car cela facilite l'interprétation et améliore le rendu.

Pour arriver au résultat, nous devons créer un pipeline de traitement pdal (enchainement des étapes à réaliser) au format json que l'on pourra ensuite appeler pour lancer une commande pdal pipeline.

Ici on définit :

1. une variable d'entrée générique : input.laz
2. le filtre de colorisation qui va nous permettre de récupérer la valeur RGB des pixels sur chacun des points du nuage : filters.colorization auquel on va ajouter la variable de l'image à utiliser
3. une variable de sortie output.laz

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

{
  "pipeline":[
    "input.laz",
      {
        "type":"filters.colorization",
        "raster":"ortho.jp2"
      },
      "output.laz"
  ]
}

Consulter le fichier

Ensuite dans le script "maître", je vais appeler mon pipeline au format json pour lancer la commande en définissant les variables correspondant aux fichiers à utiliser :

• readers.las.filename : correspond au fichier en entrée qui est associé au input.laz du pipeline
• filters.colorization.raster : correspond à l'image à utiliser pour la colorisation qui est associée au ortho.jp2 du pipeline
• writers.las.filename : correspond au fichier en sortie qui est associé au output.laz du pipeline

1
2
3
4

pdal pipeline 1_colorize.json
  --readers.las.filename=$REPER'/data_tmp/un7z/'${base}'/'$filename_LiDAR \
  --filters.colorization.raster="/home/utilisateur/Documents/traitement_LIDAR/data_ortho/34-2021-0"$ref_ortho_x"-"$ref_ortho_y"-LA93-0M20-E080.jp2" \
  --writers.las.filename=$REPER'/data_tmp/colorisation/'${base}'/'$base_lidar'_color.laz'

Distinguer le sol et le sursol#

Dans l'étape suivante, nous allons procéder à l'identification du sol et du sursol. Nous allons utiliser pdal pour filtrer les points LiDAR en utilisant une donnée vectorielle des bâtiments et différents filtres mathématiques pour isoler le sol et le reste des éléments de sursol.

Identifier les bâtiments#

Pour me "faciliter le travail" et pour vous présenter comment l'utilisation d'une donnée vectorielle peut être mobilisée pour classifier un nuage de points, j'ai décidé d'utiliser une couche vectorielle des bâtiments sur laquelle je vais ajouter une colonne où je vais attribuer un attribut de classification pour les distinguer (ici j'ai utilisé la valeur 7).

Voici comment j'utilise gdal pour créer la colonne que j'ai nommé classif :

1
2
3
4
5
6
7

ogr2ogr \
    -f "ESRI Shapefile" \
    -progress -skipfailures -overwrite \
    -dialect SQLITE \
    -nlt POLYGON \
    -sql "SELECT ST_Union(st_buffer(ST_MakeValid(geometry),0)) as geometry, 7 AS classif FROM batiments_clip WHERE st_IsValid(st_buffer(ST_MakeValid(geometry),0))" \
    batiments_classif.shp batiments_clip.shp

Ensuite le même principe que pour la colorisation, on va créer un pipeline de traitement :

1. une variable d'entrée : input.laz
2. on va assigner la valeur 0 à l'attribut classification de l'ensemble des points du nuage pour s'assurer qu'il n'y ait pas eu d'affectation antérieure
3. on va utiliser un filtre permettant de travailler sur la superposition d'une donnée vectorielle avec des données LiDAR : filters.overlay
4. en définissant notre volonté de classifier "dimension":"Classification"
5. en définissant le fichier à utiliser
6. la colonne stockant la valeur de classification à attribuer préalablement définie
7. une variable de sortie output.laz

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

{
  "pipeline":[
  "input.laz",
    {
      "type":"filters.assign",
      "assignment":"Classification[:]=0"
    },
    {
        "type":"filters.overlay",
        "dimension":"Classification",
        "datasource":"batiments_classif.shp",
        "layer":"batiments_classif",
        "column":"classif"
    },
    "output.laz"
  ]
}

Consulter le fichier

On exécute le pipeline en définissant nos données colorisées en entrée sur laquelle devra se faire la classification.

1
2
3
4

pdal pipeline 2_pipeline.json \
  --verbose 4
  --readers.las.filename=$REPER'/data_tmp/colorisation/'${base}'/'${filename_lidar%%.*}'_color.laz' \
  --writers.las.filename=$REPER'/data_tmp/filter/'${base}'/'$base_lidar'_filter.laz'

Identifier le sol#

On va maintenant utiliser l'algorithme de classification Simple Morphological Filter (SMRF) pour distinguer le sol en excluant les points préalablement classifiés comme "bâtiment".

Le pipeline de traitement est défini par :

1. une variable d'entrée : input.laz
2. un filtre d'exclusion : filters.range pour ignorer les points dont l'attribut classification est différent de 7 afin d'ignorer les bâtiments
3. un filtre de classification filters.smrf pour isoler les éléments du sol
4. un filtre de sélection filters.range pour conserver tous les points ayant été classifiés comme du sol par l'algorithme SMRF avec la valeur 2 (valeur conventionnelle utilisée pour le sol)
5. une variable en sortie output.laz

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

{
  "pipeline":[
    "input.laz",
    {
        "type":"filters.range",
        "limits":"Classification![7:7]"
    },
    {
        "type":"filters.smrf",
        "returns":"only",
        "cell":1,
        "slope":0.15,
        "scalar":1.2,
        "threshold":0.2,
        "window":22.0,
        "ignore":"Classification[7:7]"
    },
    {
        "type":"filters.range",
        "limits":"Classification[2:2]"
    },
    "output.laz"
  ]
}

Consulter le fichier

Toujours sur le même principe, on définit les variables d'entrée et de sortie pour la bonne exécution :

1
2
3
4

pdal pipeline 3_ground.json \
  --verbose 4 \
  --readers.las.filename=$REPER'/data_tmp/filter/'${base}'/'${filename_lidar%%.*}'_filter.laz' \
  --writers.las.filename=$REPER'/data_tmp/ground/'${base}'/'$base_lidar'_ground.laz'

Identifier le sursol#

Ensuite, on relance l'algorithme de classification Simple Morphological Filter (SMRF) en ajustant les paramètres pour isoler le sursol.

Le pipeline de traitement est défini ainsi :

1. une variable d'entrée : input.laz
2. un filtre d'assignation de la valeur 6 pour les bâtiments (valeur conventionnelle utilisée pour les bâtiments)
3. un filtre de classification filters.smrf pour isoler les éléments du sursol en ignorant les points classés 6
4. un filtre d'exclusion : filters.range pour ignorer les points dont l'attribut classification est 2 (points correspondant au sol)
5. un filtre d'assignation de la valeur 1 pour les éléments du sursol (valeur conventionnelle utilisée pour les éléments non classifiés)
6. une variable en sortie output.laz

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

{
  "pipeline":[
    "input.laz",
    {
        "type":"filters.assign",
        "assignment":"Classification[7:7]=6"
    },
    {
        "type":"filters.smrf",
        "cell":1,
        "slope":0.15,
        "scalar":1.2,
        "threshold":0.45,
        "window":22.0,
        "ignore":"Classification[6:6]"
    },
    {
      "type":"filters.range",
      "limits":"Classification![2:2]"
    },
    {
        "type":"filters.assign",
        "assignment":"Classification[0:1]=1"
    },
    "output.laz"
  ]
}

Consulter le fichier

On appelle le pipeline et on définit les variables :

1
2
3

pdal pipeline 4_non_ground.json
  --readers.las.filename=$REPER'/data_tmp/filter/'${base}'/'${filename_lidar%%.*}'_filter.laz' \
  --writers.las.filename=$REPER'/data_tmp/no_ground/'${base}'/'$base_lidar'_no_ground.laz'

Combinaison sol/sursol#

Voici le résultat après avoir chargé mon fichier de sol et mon fichier de sursol en leur attribuant une couleur par classe d'affectation ou en jouant avec les valeurs RGB.

Fusionner les dalles#

Maintenant que toutes les dalles ont été traitées (4 dalles dans un fichier 7zip), je décide de fusionner séparément les dalles de sol et les dalles de sursol afin de couvrir une zone plus grande et d'avoir deux fichiers distincts.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

{
  "pipeline":[
        {
            "tag" : "las1",
            "type" : "readers.las"
        },
        {
            "tag" : "las2",
            "type" : "readers.las"
        }
        ,
        {
            "tag" : "las3",
            "type" : "readers.las"
        }
        ,
        {
            "tag" : "las4",
            "type" : "readers.las"
        },
        {
            "type" : "filters.merge"
        },
        "placeholder.laz"
  ]
}

Consulter le fichier

On paramètre maintenant la fusion des différentes dalles.

1
2
3
4
5
6
7

pdal pipeline 5_merge.json
  --verbose 4 \
  --writers.las.filename=$REPER'/data_tmp/ground/'${base}'_ground.laz' \
  --stage.las1.filename=$las_ground1 \
  --stage.las2.filename=$las_ground2 \
  --stage.las3.filename=$las_ground3 \
  --stage.las4.filename=$las_ground4

Créer un raster du sol (MNT)#

La dernière étape consiste à générer un raster à partir du nuage de points classifié comme du sol en définissant :

1. une variable d'entrée : input.laz
2. les paramètres de création du raster à l'aide de gdal :
   • Format : GTiff
   • La méthode d'interpolation
   • La résolution
   • le fichier en sortie output.tif

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

{
  "pipeline":[
    "input.laz",
      {
        "type": "writers.gdal",
        "gdaldriver":"GTiff",
        "output_type":"mean",
        "resolution":"1",
        "filename":"output.tif"
        }
  ]
}

Consulter le fichier

On appelle le pipeline et on définit les variables :

1
2
3

pdal pipeline 6_ground_raster.json
  --readers.las.filename=$REPER'/data_tmp/ground/'${base}'_ground.laz' \
  --writers.gdal.filename=$REPER'/data_tmp/ground_raster/'${base}'_ground_raster.tif'

Conclusion#

Cet article a vocation à montrer les possibilités et le fonctionnement global de PDAL.

La méthode proposée pour distinguer le sol et le sursol à partir de données LIDAR est à repositionner dans le contexte de la livraison des fichiers bruts. Celle-ci est perfectible si on compare aux fichiers classifiés livrés par l'IGN mais elle a l'avantage d'être facilement adaptable dans un autre contexte et suffisamment robuste pour pouvoir traiter de vastes emprises.

Auteur·ice#

Florian Boret#

Géomaticien/cartographe, je suis arrivé dans le monde de la géomatique en suivant un cursus « professionnalisant » (BTS Géomètre-Topographe, Licence pro GGAT, Master SIGAT). J’ai ensuite travaillé dans un bureau d’études spécialisé dans la production de données d’occupation du sol et puis pour des raisons personnelles je me suis expatrié quelques années au Sénégal où je me suis lancé comme géomaticien indépendant (DATA\WAX).

Depuis mon retour en France, je suis en charge du SIG de la communauté d'agglomération Lunel Agglo.

En dehors de ces expériences, j'ai aussi régulièrement initié de petits projets personnels iGeo-Topo, GIS-Blog.fr, osm2igeo, osm2igeotopo. Aujourd'hui, c’est avec plaisir que j’interviens également comme contributeur de GeoTribu.

Licence Beerware #

Ce contenu est sous licence Beerware (Révision 42).
Les médias d'illustration sont potentiellement soumis à d'autres conditions d'utilisation.

Commentaires