Lidar : Filtrage du sol et Création de MNT

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 :

  1. réassigner tous les points à la classe 0 car on suppose que la classification d’origine n’est pas bonne
  2. utilisation de l’algo ELM pour assigner les points considérés comme du bruit à la classe 7
  3. faire un filtrage statistique pour les points abbérants
  4. utiliser le filtre SMRF sur les points restants pour classer les points
  5. 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:

AWK traitement de fichier texte

La commande unix awk permet des manipulations avancées sur les chaînes de caractères.

awk [options] [programme] [fichier]

L’option la plus courante est -F : séparateur de champ

La structure du programme est la suivante :
'motif1 { action1 }
motif2 { action2 } …'

motif (en anglais pattern) est la condition pour réaliser l’action. On utilise le plus souvent une expression régulière.

Exemple 1 : condition

avec un fichier ASCII “input.xyz” de la forme suivante (données Lidar) :

13.83535621 42.12214528 401.940 32 2
13.83534713 42.12213870 401.870 28 2
13.83533818 42.12213222 401.670 33 2
13.83532879 42.12211645 401.950 30 1
13.83533763 42.12212258 401.630 28 2
13.83534682 42.12212933 401.830 29 2
13.83535625 42.12213609 402.250 25 1

on veut ne conserver que les lignes où la 5ème colonne est différente de 1

awk '{ if ( $5!=1 ) print $0 }' input.xyz > output.xyz

et voilà! c’est tout

Exemple 2 : utilisation dans un script

on a un fichier texte avec différents espacements entre colonne, parfois un espace, parfois plusieurs. C’est le cas des fichiers en sortie du logiciel Circé IGN qui sont formattés vraiment bizarrement. Il contiennent parfois un et parfois plein d’espaces entre les champs, et des lignes d’entête qui commencent par *! , voici un extrait d’un fichier en sortie de Circé.

*!COORDONNÉES: projetées
*!ELLIPSOÏDE: GRS 1980
*!PROJECTION: LAMBERT-93
*!REPERE VERTICAL: IGN69
*!
*! id E N H geod.prec. vert.prec.
point1 830488.865    6252042.364 0.558 -1.167384     542.6 no info 5 / 10 cm
point2 831546.415    6251863.871 6.043 -1.176810     544.3 no info 5 / 10 cm

Voici la manip pour reformatter cela simplement avec awk au sein d’un script.

#! /bin/bash

awk '$1!~/*!/ {$1=$1;print $0}' OFS=, > out.csv

pour utiliser le script, on tape simplement: cat input.txt | ./script.sh

encore une fois, c’est tout!

Quelques explications :

  • la variable $1 correspond au premier champ de l’enregistrement,
  • la variable $0 correspond à tout l’enregistrement,
  • la variable définie par OFS correspond au séparateur de champ voulu en sortie, ici on définie une virgule
  • L’action {$1=$1;print $0}effectue d’abord un calcul de champ avec $1=$1, ici on affecte au premier champ la valeur du premier champ. Donc on ne fait rien dans ce calcul, mais cela a pour effet de ré-écrire tout l’enregistrement avec le séparateur de champ voulu (par défaut un espace), et ainsi cela supprime tous les espaces superflus entre les champs et en début et fin d’enregistrement. Et ensuite avec print $0  on écrie tout l’enregistrement dans la sortie standard.
  • Le motif $1!~/*!/ sélectionne tous les enregistrements dont le premier champ ne contient pas à la chaîne de caractères *!

Exemple 3 : supprimer tous les espaces superflus

awk '{$1=$1}1' input.csv > output.csv

  • L’action {$1=$1} voir Exemple 2 ci-dessus
  • Le 1 situé après les accolades {} est un motif correspondant à Vrai, et lorsqu’il n’y a pas d’action située après le motif c’est l’action par défaut print qui est executée, donc chaque ligne est écrite dans la sortie standard.

Exemple 4 : chaîner les commandes en bash

exemple tiré de https://www.the-art-of-web.com/system/imagemagick-watermark/

ls -1 photos/process_*.jpg | awk -F\/ '{print "composite -gravity SouthEast -dissolve 50% -resize 400 watermark.png photos/"$(NF)" watermarked/"$(NF)}' | sh

Cette commande permet d’ajouter une trame en surimpression (watermark) sur toute une série d’images de façon automatique avec l’outil “composite” (un des outils d’Imagemagick). awk sert ici à écrire la commande de “composite” avec tous les fichiers renvoyés par la commande “ls …”.

Cet exemple est intéressant car on voit bien l’intérêt de awk pour scripter des tâches. On chaîne ici trois commandes avec des pipes. La première liste les photos correspondant à un certain pattern, la seconde écrit une commande “composite” et la troisième exécute cette commande “composite”.

Quelques exemples tirés de la page Wikipedia

  • awk '{print $0}' fichier : affiche toutes les lignes de fichier (idem que cat fichier).
  • awk '/2/ {print $0}' fichier : affiche toutes les lignes où le caractère 2 est présent (idem que grep '2' ref.txt).
  • awk '$1~/2/ {print $0}' fichier : affiche toutes les lignes où le caractère 2 est présent dans le premier champ.
  • awk '{print NR ":", $0}' fichier : affiche le contenu de fichier, mais chaque ligne est précédée de son numéro.
  • awk -F: '{print $1}' /etc/passwd : renvoie la liste des utilisateurs (idem cut -d : -f 1 /etc/passwd).
  • awk 'BEGIN {FS = ":"}{print $1}' /etc/passwd : idem que la précédente commande
  • awk '{s=s+$1} END {print s, s/NR}' fichier : écrit la somme et la moyenne de tous les nombres de la première colonne de fichier.
  • awk '/Motif1/ , /Motif2/' fichier : écrit toutes les lignes contenues dans le fichier entre le Motif1 et le Motif2.

References:

Emprise d’un raster

Conversion en vecteur de l’emprise valide (sans nodata) d’un raster

remarque générale: les fichiers rasters intermédiaires produits peuvent être écrits en format VRT (raster virtuel). Cela permet d’avoir de tous petits fichiers.

Etape 1 : Conversion en 1 bit

exemple avec un fichier in.tif en entrée, et écriture d’un fichier out.vrt en sortie

gdal_translate in.tif out.vrt -ot byte -of VRT -co NBITS=1

l’option -co NBITS=1 permet de passer l’image en binaire

pour écrire un fichier en sortie au format VRT, utiliser l’option -of VRT

Etape 2 : Création d’un canal alpha

gdalwarp -dstalpha out.vrt out1.vrt -co NBITS=1

Le canal alpha est créé avec toutes les données en nodata grâce à l’option -dstalpha

on peut aussi prendre en compte une valeur représentant le nodata dans le fichier en entrée grâce à l’option -srcnodata value

Bien vérifier le raster out1.vrt en sortie pour voir dans quel canal se trouve la bande alpha à utiliser par la suite.

Etape 3: Vectorisation du raster

gdal_polygonize.py out1.vrt -b 1 -f "ESRI Shapefile" out2.shp

l’option -b 1 permet de prendre ce qui est dans la bande 1

Methode Alternative : avec gdal_calc

on calcule un raster avec une expression logique là où il y a des valeurs

gdal_calc.py -A input --outfile=tmp.tif --type=Byte --calc="isnan(A)" 

on convertit en polygone le raster:

gdal_polygonize.py tmp.tif -f "ESRI Shapefile" output

Conversion ASCII XYZ en MNT

On suppose dans cette procédure que nos données constituent un maillage régulier et qu’il n’est pas nécessaire de faire d’interpolation. Le fichier ASCII XYZ est composé de trois colonnes X, Y et Z, ex.

880000 6340000 828
880025 6340000 831
880050 6340000 833
880075 6340000 835

Procédure avec ArcGIS

1-Changer l’extension du fichier en .txt ou .csv

2- Modifier le fichier schema.ini

Plus d’infos sur ce fichier schema.ini à http://www.tectonics.caltech.edu/gi…

Ex. de fichier schema.ini correspondant au fichier input.csv ci-dessus : [input.csv]

ColNameHeader=False
Format=Delimited( )
Col1=X Double
Col2=Y Double
Col3=Z Double
  • si le fichier contient une ligne d’entête, on aura ColNameHeader=True
  • si le séparateur est la virgule on aura Format=CSVDelimited
  • si le séparateur est la tabulation on aura Format=TabDelimited
  • si le séparateur est l’espace on aura Format=Delimited( )
  • si le séparateur est | on aura Format=Delimited(|)
  • … etc avec d’autres séparateurs

Si le fichier n’est pas bien reconnu par ArcCatalog, il est possible que ce soit un problème de codage des caractères. On pourra le modifier en ajoutant la ligne suivante dans schema.ini :

CharacterSet=UNICODE

Changer UNICODE par le codage correct (ANSI, autre …)

3- Vérifier la reconnaissance du fichier par ArcCatalog

  • Vérifier dans ArcCatalog que votre fichier est bien reconnu. Pour cela votre fichier doit s’afficher dans l’onglet Aperçu.
  • Vérifier aussi que le type des champs est correct. Pour cela faire un clic droit sur le fichier et vérifier l’onglet Champs

4- Créer le shapefile de points

  • Pour créer les points, utiliser la fonction d’Arctoolbox « Outils 3D Analyst Depuis un fichier – ASCII 3D vers classe d’entités »

5-Convertir le shapefile de points en raster

  • Utiliser la fonction d’ArcToolbox « Outils de conversion – Vers raster – Entité vers raster »
  • Prendre le champ Shape.Z ou le champ Z comme valeur de raster
  • Définir la taille de cellule selon l’espacement des points en entrée.

Procédure avec Linux et GDAL

ex. conversion d’un fichier XYZ ASCII en TIF et assignation de la référence spatiale RGF93-Lambert93 (Code EPSG 2154)

gdal_translate -a_srs EPSG:2154 input.XYZ output.tif

en batch on pourra faire un script du type :

#!/bin/bash
#Ce script assigne le code EPSG saisi en paramètre aux images en sortie
# Jules Fleury, SIGeo/CEREGE, mars 2011

# saisir le code epsg en parametre

# Assign_SRS_to_TIF.sh

echo “Input spatial reference EPSG (exemple: for WGS84 enter 4326) : ”
read c_EPSG


for FILE in *.XYZ # modifier ici en cas d'autre format
do
BASE=`basename $FILE .XYZ` # modifier ici en cas d'autre format
NEWFILE=${BASE}_c.tif
gdal_translate -a_srs EPSG:$c_EPSG $FILE $NEWFILE
done

Avec Linux et GDAL : cas ou les données en entrée ne sont pas bien triées en X ou Y

solution prise ici : http://osgeo-org.1560.x6.nabble.com/gdal-dev-Opening-gridded-xyz-data-that-is-out-of-order-td5334341.html

if we obtain this kink of error:

> ERROR 1: Ungridded dataset: At line 4125001, change of Y direction

> gdalinfo failed – unable to open ‘DTM_swissALTI3D_XYZ.txt’.

You must sort the data beforehand

– the y position must be sorted properly (ASC or DESC)
Standard linux sort program:
sort -k2 -n -k1 input.xyz -o output_sort.xyz
(2nd column (y) as numeric, then first column (x) -o = output file
puis on relance gdal_translate
gdal_translate input.xyz output.tif -a_nodata -9999
Remplacer par la bonne valeur pour nodata
Si cette méthode ne fonctionne pas, on fera une interpolation comme ci-dessous

Avec Linux et GDAL : cas ou les données ne couvrent pas une dalle entière

On utilise la fonction d’interpolation de GDAL gdal_grid pour ne prendre pour chaque pixel qu’un seul point en entrée. La méthode consiste à premièrement créer un CSV, puis un VRT et utiliser ce fichier comme source de gdal_grid.

ex. gdal_grid -a nearest:radius1=1:radius2=1:nodata=-9999.0 -ot Float32 -a_srs EPSG:2154 -txe 1074000 1079000 -tye 6325000 6330000 -outsize 1000 1000 -l grid grid.vrt grid_out.tif

traitement en batch

1ère étape : conversion des XYZ en CSV avec création de VRT associé

ex. de script

#!/bin/bash

#Lecture de la référence spatiale
echo “Input spatial reference EPSG (exemple: for WGS84 enter 4326, for Lambert93 2154, for L2E 27572) : ”
read c_EPSG

#traitement des fichier en batch
for IFile in *.xyz # modifier ici en cas d’autre format
do
BASE=`basename $IFile .xyz`
OFile=${BASE}.csv

#ajout d’une ligne d’entête “X Y Z”
sed ‘1i\X Y Z’ $IFile > $OFile

#remplacement du caractère espace par une tabulation
TFile=”Temp”
tr ‘ ‘ ‘\t’ < $OFile > $TFile
mv $TFile $OFile


#création du fichier VRT
OFileVRT=${BASE}.vrt
echo "<OGRVRTDataSource>" > $OFileVRT
echo "    <OGRVRTLayer name=\"$BASE\">" >> $OFileVRT
echo "        <SrcDataSource>$OFile</SrcDataSource>" >> $OFileVRT
echo "        <GeometryType>wkbPoint</GeometryType>" >> $OFileVRT
echo "        <LayerSRS>EPSG:$c_EPSG</LayerSRS>" >> $OFileVRT
echo "        <GeometryField encoding=\"PointFromColumns\" x=\"X\" y=\"Y\" z=\"Z\"/>" >> $OFileVRT
echo "    </OGRVRTLayer>" >> $OFileVRT
echo "</OGRVRTDataSource>" >> $OFileVRT
done

2ème étape : conversion des VRT (avec leur CSV) en TIF

on utilise gdal_grid,ex. de script

#!/bin/bash

#Programme pour interpoler des fichiers de points XYZ stockés en CSV/VRT
#selon différents algorithmes et paramètres.
#le fichier grid en sortie est stocké en TIF
#puis rééchantilloné selon la résolution voulue.

#Jules Fleury – SIGéo/CEREGE
#avril 2013

#Lecture de la référence spatiale
echo “Input spatial reference EPSG (exemple: for WGS84 enter 4326, for Lambert93 2154, for L2E 27572) :”
read c_EPSG

#Traitement des fichiers en batch
for IFile in *.vrt # modifier ici en cas d’autre format
do
echo “Processing $IFile”
BASE=`basename $IFile .vrt`
OFile=${BASE}.tif
val_nodata=-9999 #valeur pour nodata

#Traitement du fichier en utilisant l’algorithme Nearest
gdal_grid -a nearest:radius1=0.5:radius2=0.5:nodata=$val_nodata -ot Float32 -a_srs EPSG:$c_EPSG  -l $BASE $IFile $OFile


#ré-échantillonage du fichier en sortie
ressize=10 #valeur de ré-échantillonage
OFileRes=${BASE}_$ressize.tif #fichier ré-échantilloné en sortie
gdalwarp -tr $ressize $ressize -r cubic $OFile $OFileRes
echo "Done"
done

Conversion fichier ASCII XYZ en Shapefile

On fera attention au séparateur de champs du fichier (virgule par défaut, tabulation, espace) et à la présence ou non d’une entête.

2 possibilités, en utilisant ArcGIS ou en utilisant OGR

Avec ArcGIS

1- Renommer votre fichier avec l’extension .csv ou .txt

2- Modifier le fichier schema.ini

Plus d’infos sur ce fichier schema.ini à http://www.tectonics.caltech.edu/gi…

Ex. de fichier schema.ini correspondant au fichier input.csv ci-dessus :

[input.csv]
ColNameHeader=False
Format=Delimited( )
Col1=X Double
Col2=Y Double
Col3=Z Double
  • si le fichier contient une ligne d’entête, on aura ColNameHeader=True
  • si le séparateur est la virgule on aura Format=CSVDelimited
  • si le séparateur est la tabulation on aura Format=TabDelimited
  • si le séparateur est l’espace on aura Format=Delimited( )
  • si le séparateur est | on aura Format=Delimited(|)
  • … etc avec d’autres séparateurs

Si le fichier n’est pas bien reconnu par ArcCatalog, on peut définir le codage des caractères en ajoutant la ligne suivante dans schema.ini : CharacterSet=UNICODE Changer UNICODE par le codage correct (ANSI, autre …)

3- Vérifier la reconnaissance du fichier par ArcCatalog

Vérifier dans ArcCatalog que votre fichier est bien reconnu. Pour cela votre fichier doit s’afficher dans l’onglet Aperçu. Vérifier aussi que le type des champs est correct. Pour cela faire un clic droit sur le fichier et vérifier l’onglet Champs

4- Créer le shapefile de points

Clic droit sur le fichier – Créer une classe d’entités – A partir d’une table XY – Renseigner les champs contenant les coordonnées X, Y et éventuellement Z ; renseigner aussi le système de coordonnées ; spécifier le fichier en sortie – Cliquer sur OK Vérifier votre shapefile.

Avec linux et OGR

Plus d’infos sur http://www.gdal.org/ogr/drv_csv.html et http://gdal.gloobe.org/ogr/formats/…

1- Renommer votre fichier avec l’extension .csv

2- Ajouter la ligne d’entête La ligne d’entête doit contenir les noms des colonnes séparés par le délimiteur. Ex.

X Y Z
869324.598 1783031.225 125.000
869781.445 1783057.845 126.000
869340.038 1782760.764 127.000

Pour cela, on peut soit éditer le fichier avec un éditeur de texte soit utiliser une commande sed dans le shell : Ex. (ajout d’une ligne ’’X Y Z’’ en début de fichier)

sed '1i\X Y Z' input.csv > output.csv

3- Remplacement du caractère Espace par le caractère Tabulation

OGR ne reconnaissant par le caractère Espace comme séparateur, il faut le remplacer par une tabulation ou une virgule. Pour cela on utiliser une commande shell :

tr ' ' '\t' <input.csv >>output.csv

on peut aussi utiliser sed, exemple :

sed 's/[ ]/\t/g' test3.csv > test3_b.csv

4- Lecture d’un CSV contenant des données spatiales

Il est possible d’extraire de l’information spatiale (points) d’un CSV ayant des colonnes X et Y en utilisant le driver VRT (http://www.gdal.org/ogr/drv_vrt.html ) Il faut écrire le fichier VRT associé (input.vrt) au fichier CSV (input.csv)

<OGRVRTDataSource>
<OGRVRTLayer name="input">
<SrcDataSource>input.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>EPSG:27572</LayerSRS>
<GeometryField encoding="PointFromColumns" x="X" y="Y"/>
</OGRVRTLayer>
</OGRVRTDataSource>

et ogrinfo input.vrt renverra :

INFO: Open of `input.vrt'
using driver `VRT' successful.
1: input (Point)

5- Conversion du CSV en Shapefile

ogr2ogr output.shp input.vrt

On peut ajouter un filtre spatial sur les coordonnées en utilisant “-spat”

ogr2ogr -spat 822500 6245000 836500 6264000 -a_srs epsg:2154 Veran_T_total.shp T_total.vrt

Correction différentielle des données GPS ArcPad

Introduction

Cette procédure traite de la correction différentielle des associés au fichier SSF, acquis avec des appareils Trimble Juno SB, équipés des logiciels ArcPad8 et GPSCorrect.

Etapes

0.Le fichier SSF ne contient que des positions et non de l’information attributaire. Les informations attributaires sont contenues dans les fichiers SHP ou les AXF.

Pour réaliser la correction différentielle, vous devrez corriger le fichier SSF et ensuite réaliser la synchronisation entre le fichier corrigé COR et le fichier SHP ou AXF.

1. Copier les données de votre PAD (ex. Trimble Juno SB) sur le poste local. Mettre le fichier .ssf et les shapefiles dans un même répertoire.

2. Pour effectuer la correction différentielle, la méthode est d’utiliser l’assistant de correction de GPS Pathfinder Office disponible par le menu “Outils/Correction différentielle”. Vous obtenez alors un fichier GPSCorrect.cor.

A noter dans cet assistant, l’importance d’inclure dans le fichier “Cor” à la fois les positions corrigées et non corrigées (paramètres disponibles dans l’assistant, à la fenêtre “Paramètres de correction”, en cliquant sur “Changer”, dans l’onglet “Sortir”). Ce paramètre permet d’obtenir une correspondance exact entre le fichier SSF et le fichier COR.

3. Avant de procéder à la synchonisation, assurez-vous que les fichiers SHP (DBF, SHX….ou AXF seul) se trouvent dans un répertoire spécifique, dans lequel se trouvent également les fichiers SSF et COR….

4. Aller ensuite dans le dossier “Outils/Autre/ShapeCorrect” afin de lancer l’assistant de synchronisation.

6. Dans “ShapeCorrect”, cliquer sur “Parcourir” afin de sélectionner les fichiers SHP (ou le fichier AXF). Il n’est pas nécessaire de sélectionner les fichiers SSF et/ou COR, mais ils se doivent d’être présents dans le meme répertoire.

7. C’est à vous de décider si vous souhaitez dans les fichiers résultants les positions corrigées et/ou non corrigées.

8. Cliquer sur “OK” afin de lancer la synchro. Vérifiez le résultat de la synchronisation.