Gros rasters et Compression

Lors de l’écriture d’un raster, surtout en THR (très haute résolution), il est important de compresser au maximum les images pour avoir des fichiers qui soient gérables. Il est important d’avoir en tête que certaines compression se font avec une perte de qualité, et d’autres sans perte. Si l’on veut simplement faire de la visualisation, on peut se permettre un peu de perte, par contre pour du traitement ou pour des données type MNT, on ne voudra aucune perte d’information.

Toutes ces histoires de compression sont bien expliquées sur ce blog http://blog.cleverelephant.ca/2015/02/geotiff-compression-for-dummies.html

Tuilage

Possibilité de tuiler (daller) pour permettre de gérer des trop grosses images. C’est parfois indispensable, par exemple si l’on passe au dela de 20 à 30000 pixels de large pour une image.

On pourra consulter pour plus d’infos cette page https://sigeo.cerege.fr/?p=116

Le BigTIFF et les tuiles internes

le GeoTIFF permet la gestion de très grosses images avec un système de tuilage interne. Tout est donc stocké dans 1 fichier TIFF.

Pour les images >= 4Go, il faudra les écrire en BigTiff.

Parfois cela est réalisé automatiquement par l’application qui va écrire l’image. Si l’on utilise GDAL, il est préférable de le signifier explicitement avec l’option “-co BIGTIFF=YES”

pour forcer l’écriture des tuiles en interne, c’est avec l’option -co TILED=YES

exemple de commande avec ces 2 options: gdal_translate im_in.tif im_out.tif -co BIGTIFF=YES -co TILED=YES

Compression sans perte de qualité

par exemple pour un MNT ou de l’imagerie que l’on veut garder intacte, j’utilise la compression avec DEFLATE

ex . gdal_translate dem_in.tif dem_out.tif -co COMPRESS=DEFLATE

Compression avec perte de qualité

on peut utiliser la compression JPEG pour réduire efficacement la taille des images multi-canaux en 8 bits (ne fonctionne pas pour les images 16 ou 32 bits) :

gdal_translate \
-co COMPRESS=JPEG \
-co PHOTOMETRIC=YCBCR \
-co TILED=YES \
input.tif output.tif

Le stockage en JPEG est encore plus performant dans l’espace de couleurs YCBCR.

Pour des images en 16 ou 32 bits, on utilisera LZW ou PACKBITS pour la compression avec perte de qualité.

Compression des pyramides raster

Les pyramides permettent d’accélérer l’affichage en créant des aperçus selon les niveaux de zoomes choisis. Ici la perte de qualité n’est pas un problème car c’est simplement pour l’affichage, et la donnée originale n’est pas altérée.

Les pyramides peuvent être stockées dans un fichier séparé (format OVR), ce qui est en général préférable pour plus de clarté, ou en interne dans un GeoTiff.

exemple de pyramides pour de l’imagerie 8 bits :

gdaladdo \
--config COMPRESS_OVERVIEW JPEG \
--config PHOTOMETRIC_OVERVIEW YCBCR \
--config INTERLEAVE_OVERVIEW PIXEL \
-r average \
-ro \
im_input.tif \
2 4 8 16 32

pour un MNT (16 ou 32 bits en général), on fera plutôt avec du DEFLATE ou du LZW

gdaladdo --config COMPRESS_OVERVIEW LZW dem_input.tif 2 4 8 16 32 64

Références

Compression pour dummies :  http://blog.cleverelephant.ca/2015/02/geotiff-compression-for-dummies.html

Format GeoTifff : https://gdal.gloobe.org/gdal/formats/gtiff.html

Format Geotiff : https://www.gdal.org/frmt_gtiff.html

gdaladdo pyramides : https://www.gdal.org/gdaladdo.html

 

Traitement MicMac d’images satellites

Workflow complet pour traiter des images satellites (multi-) stéréoscopiques

avec MicMac pour le traitement de photogrammétrie, et GDAL et OTB pour les post-traitements

Génération d’images Fusionnées (Pan-Sharpen) ortho-rectifiées avec MicMac

Le process est le suivant, utilisable avec le script https://github.com/JulesFleury/JulesFleury/blob/master/Script_Decomposition_Band.sh

#0- réaliser la fusion du bundle d’origine avec otb_bundletoperfectsensor

#1-
split Pan-sharpen bands
avec otbcli_splitimage.....

#2-
copie et renommage des orientations initiales et ajustées de micmac dans le dossier Ori-XXX-adj (cf. script)

#3-
#relancer MicMac Malt pour générer les ortho sur les différents canaux, plusieurs options sont possibles

  • a- si Malt n’a pas encore tourné

#mm3d Malt Ortho "Pattern de toutes les images" "dossier_orientations" DirMEC="dossier de la MEC des panchromatiques" InMNT="Pattern des panchromatiques" DoMEC=1 DoOrtho=1 ImOrtho="Pattern des différentes bandes multispectrales" DirOF="dossier des orthophotos" EZA=1

  • b- si Malt a déjà tourné sur les panchro

#mm3d Malt Ortho "Pattern de toutes les images" "dossier_orientations" DoMEC=0 DoOrtho=1 ImOrtho="Pattern des différentes bandes multispectrales" DirOF="dossier des orthophotos pour le multispectral"
#exemple : mm3d Malt Ortho "IMG.*.TIF" Ori-RPC-d0-adj DirMEC="MEC-Malt-Zf4" DoMEC=1 ImMNT="IMG.*(00[123]).TIF" ImOrtho="IMG.*_b.*.TIF" DirOF="Ortho_Fusion" DoOrtho=1 ZoomF=4 EZA=1

#4-
#Rassembler les différents canaux des images pan-sharpen
#otbcli_ConcatenateImages ....

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