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

 

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

Traiter les dalles SRTM avec GDAL

Cette procédure décrit la méthode pour télécharger et prétraiter les dalles SRTM v4.

Téléchargement :

se connecter sur le serveur CGIAR dans la partie SRTM Data Search and Download et choisir une ou plusieurs dalles de 5°x5°

on obtient une archive contenant 4 fichiers (ex. readme.txt ; srtm_41_05.hdr ; srtm_41_05.tif ; srtm_41_05.tfw)

Mosaïque :

Déplacer tous les fichiers dans un seul répertoire.
Puis lancer la commande :
gdal_merge.py -n -32768 -o out_filename input_files

- n -32768 : pour fixer les valeurs nodata
out_filename : nom du fichier en sorite
input_files : liste de fichiers en entrée

Plus d’options GDAL sur http://www.gdal.org/index.html

Découpage d’une zone

gdal_translate -projwin ulx uly llx lly) src_dataset dst_dataset
ulx uly llx lly sont les coordonnées projetées de la zone à découper, ex. en WGS84 sur la Crête 23 36 27 34
src_dataset : fichier en entrée
dst_dataset : fichier en sortie

Plus d’options GDAL sur http://www.gdal.org/index.html

Ombrage (ou estompage)

gdaldem hillshade -s 111120 input_dem output_hillshade

- s 111120 : ratio pour les unités verticales dans le cas où les coordonnées sont en degrés (e.g. Lat/Long WGS84) et les altitudes en mètres.

Plus d’options GDAL sur http://www.gdal.org/index.html

Dallage avec GDAL

Générer des tuiles avec gdal2tiles.py

création d’un tuilage avec cet utilitaire gdal ici

Description d’une procédure complète pour créer un tuilage et le mettre au format KML de Google Earth ici

Re-tuiler des tuiles existantes

utilitaire ici

exemple pour générer les tuiles en JPEG, compressées à un taux de 80, à partir d’une image ECW

gdal_retile.py -of JPEG -ps 2000 2000 -s_srs EPSG:27572 --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL --config JPEG_QUALITY_OVERVIEW 80 -tileIndex TI_ortho_2003.shp -tileIndexField TileName -targetDir /tmp/tileortho/ /mnt/srvstk0_geomor/Data_SIG/OHM_Mine/Fonds_Reference/Ortho_2003/Ortho_2003_L2e.ecw

voir la documentation sur le site GDAL (rubrique Utilities)

Gestion des dalles raster avec GDAL

Ce document décrit différents scripts basés sur GDAL en Shell pour des dalles raster.
L’exemple d’utilisation est la récupération de données auprès de l’IGN, au format TIFF ou ECW

 

Afficher l’information sur une image raster

gdalinfo src_dataset

Affiche l’information sur un jeu de données raster. Permet notamment d’afficher la référence spatiale et l’étendue des données, le nombre de pixel, la taille du pixel. Documentation officielle gdalinfo

Affecter une référence spatiale à un ensemble de dalles

gdal_translate -a_srs srs_def src_dataset dst_dataset

Affecte le système de coordonnées « srs_def » à l’image en sortie. On utilisera la syntaxe EPSG:n (ex. EPSG :2154 pour du Lambert93 RGF93, 27572 pour du NTF Lambert II Etendu). ex. pour du NTF- Lambert II Etendu :

gdal_translate -a_srs EPSG:27572 im_src.tif im_dst.tif

Script pour un lot de dalles :

#!/bin/bash
for FILE in *.tif
do
BASE=`basename $FILE .tif`
NEWFILE=${BASE}_c.tif
gdal_translate -of GTiff -a_srs EPSG:2154 $FILE $NEWFILE
done

options de gdal_translate

  • compresser une image : l’option -co COMPRESS=DEFLATE permet de compresser l’image TIFF sans perte de qualité. L’option -co TILED=YES permet de spécifier que l’on fait un tuilage.
  • codage de l’image : l’option -co NBITS=1 permet de coder l’image en sortie à 1 bit. Documentation officielle gdal_translate

    Création d’une mosaique au format « Virtual Raster » VRT

    gdalbuildvrt mosaic.vrt *.tif

Crée un mosaïque « mosaic.vrt » à partir d’un ensemble d’une collection d’images tif. Ce format vrt est très léger et semblable aux algorithmes d’Er-Mapper, à savoir que les données d’origine sont conservées et inchangées. Le vrt est simplement un catalogue d’images. Toutes les images de la collection doivent avoir le même système de coordonnées. Pour éviter les zones noires là où il n’y a pas de données on rajoutera une option « -addalpha ». De cette manière, avec un éditeur d’image gérant le canal alpha (ex. ArcMap), on aura de la transparence dans les zones sans données source, et de l’opacité dans les zones avec données source. Ex. :

gdalbuildvrt -addalpha mosaic.vrt *.tif

Documentation officielle gdalbuildvrt

Construire des pyramides

Cette fonction permet d’accélérer l’affichage en créant des aperçus selon les niveaux de zoomes choisis. C’est l’équivalent de la création de pyramides dans ArcGIS.

gdaladdo -ro data_src 2 4 8 16 32

L’option « -ro » permet de créer un aperçu en tant qu’image externe et non inclus dans l’image tiff elle-même. Attention, si vous n’utilisez pas l’option -ro, les pyramides seront construites dans l’image d’origine et elle pourra être corrompue.
Pour compresser le plus possible en JPEG, on utilisera ce type d’options :

gdaladdo -ro --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config JPEG_QUALITY_OVERVIEW 80 data_src 2 4 8 16

Documentation gdaladdo

gdal_merge : mosaiquage

Mosaïque automatique d’un lot d’images. Toutes les images doivent être dans le même système de coordonnées et avoir le même nombre de bandes. Les images peuvent se chevaucher, et avoir des résolutions différentes. Dans les zones de chevauchement, les valeurs prendront celles de la dernière image.

gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*
[-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-v] [-pct]
[-ul_lr ulx uly lrx lry] [-n nodata_value] [-init "value [value...]"]
[-ot datatype] [-createonly] input_files

Exemple pour des dalles SRTM en Geotiff :

gdal_merge.py -o srtm_global.tif *.tif

Aide sur gdal_merge.py

gdal_translate

Avec l’utilitaire gdal_translate, les rasters peuvent être convertis dans différents formats, et certaines opérations comme le découpage, le rééchantillonnage, et la mise à l’échelle peuvent être réalisées durant le traitement.

Exemple : découpage d’un MNT SRTM sur une zone géographique, mise à “nodata” des valeurs -32768 et affectation du système de coordonnées WGS84

gdal_translate -a_nodata -32768 -a_srs EPSG:4326 -projwin 23.4 35.8 26.5 34.7 srtm_global.tif srtm_crete.tif

Aide de gdal_translate