Projections et conversions rasters

Je suis confronté à des problèmes de projections et de format d'images rasters. Le format "raster" est un des piliers du  domaine de la géomatique. L’image raster est constituée de petits rectangles pixels géo-localisés. A chaque pixel est associé une ou plusieurs valeurs appelées « bande raster » décrivant les caractéristiques de l'espace. Il peut s’agir par exemple dans une image couleur à 3 bandes de l'intensité lumineuse des trois couleurs: rouge, vert, bleu. Mais les bandes d’une image raster peuvent contenir des attributs très variés que l’on souhaite représenter à très grande échelle : altitude, exposition au risque d’inondation, population…

Les fonds de cartes tuilées sont le résultat de l’accolage de multiples images appelées aussi dalles rasters, exactement comme les carreaux de votre salle de bain mais sans les jointures. En version "tuilée", les images rasters géo-référencées constituent le support des fonds de cartes  des nombreux cartes et globes interactifs qui fleurissent sur le web.

Le problème :

L'IGN m'a transmis une série de plusieurs dizaines de milliers de ses images rasters. Le vénérable institut produit et maintient en particulier des photos aériennes issues de la BDOrtho et des plans cadastraux rasters de la BDParcellaire. L'IGN fournit ces rasters en projection Lambert 93 (L93 : EPSG 2154) et en format d'image géo-référencées .jp2 ;  format compressé JPEG2000 adapté aux applications internet. Certains SIG ou logiciel "maison" ne sont cependant pas paramétrés pour lire des rasters en projection L93 et sur le format JPEG2000.  Mon objectif est donc de convertir les images de l'IGN dans l'ancienne projection officielle française Lambert 2 Etendu (L2E : EPSG 27572 ou mieux IGNF:LAMB2C) et au format raster plus classique géoTif.

Ce problème est très spécifique et pas folichon ! Je note ici cependant mes solutions car j'ai beaucoup tâtonné et pour tout dire "galéré" pour trouver une solution acceptable. Il y a de grandes chances que les pistes exposées n'intéressent que moi, mais sait-on jamais !  D'autres travaux de conversions de formats et de re projections de rasters peuvent en bénéficier.... Pour les aspects théoriques sur le géoréféncement, la géodésie et les systèmes de projection, je renvoie vers le très instructif dossier du CEREMA (ex CERTU). Pour les aspects pratiques, c'est plus bas.

Pré-requis techniques :

Sous Windows, tout se passe en ligne de commande DOS. Il faut installer sur votre PC :

• Obligatoirement GDAL (Geospatial Data Abstraction Library) : le "couteau Suisse" du géomaticien. Ce  logiciel est open source et se substitue très bien aux équivalents payants (FME, ArcGISS....) pour ceux qui manipulent le mode commande DOS. L’installation est expliquée sur le site OSGeo4W. On peut aussi se référer à ce mode d'emploi. En respectant scrupuleusement ces consignes d'installation, il faut absolument faire attention à : 
• bien installer une version de GDAL et de python cohérentes ;
• ajouter à la variable d'environnement  PATH (accessible Panneau de configuration+Système+parametres système avancés+Variable d'environnement + Variables système)  les chemins d'accès à votre installation GDAL (sur mon PC: C:\Program Files (x86)\GDAL) ,  python et scripts python ;
• bien configurer les nouvelles variables d'environnement pour GDAL en fonction de la localisation de votre installation : GDAL_DATA (chez moi C:\Program Files (x86)\GDAL\gdal-data),  PROJ_LIB et GDAL_DRIVER_PATH .
• Optionnellement, le logiciel KAKADU La version gratuite de Demo téléchargeable ici suffit pour gérer les conversions dont nous aurons besoin.
• J'ai contrôlé les résultats des conversions sous QGIS et MapInfo.
• Pour la production de masse, j'utilise SAS pour générer des fichiers .bat. Des commandes Python, R ou même Excel font bien l'affaire aussi. Je ne détaille cependant pas vraiment ce point dans ce post.

Pour avancer dans la suite de ce post, il faut donc être à l’aise avec les SIG et aussi l’écriture et l’automatisation de commande DOS.

1/ Conversion dalle à dalle de L93.jp2 en  L2E.tif :

GDAL possède une panoplie de commandes de traitement de données SIG, il gère de nombreux format raster/vecteurs et de nombreux systèmes de projections. La librairie GDAL est parfaitement documentée en anglais et en français.

Tout se passe avec les deux commandes gdalwarp et gdal_translate. Gdalwarp fait l'essentiel du travail en passant d'une projection source (-s_srs) vers la projection cible (-t_srs) en convertissant le format vers le géotiff (-of GTiff ). Il faut absolument mentionner que la transparence de la dalle cible (-dstnodata) sera gérée via un code RGB au triple 255 (valeur du blanc) ou triple 0 (valeur du noir)

c:

cd c:\Program Files\GDAL

gdalwarp -overwrite -s_srs EPSG:2154 -t_srs EPSG:27572 -of GTiff -co TILED=YES -r average -dstnodata "255 255 255" Raster_L93.jp2 RasterV0_L2E.tif

gdal_translate  -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 -of GTiff RasterV0_L2E.tif Raster_L2E.tif

del  RasterV0_L2E.tif

2/ Conversion dalle à dalle L93.jp2 en L2E.jp2

J'ai cherché longtemps dans les dédales de la géomatique sur le web la solution de ce problème apparemment bénin. Il ne l'est pas. La piste est formulée dans cet échange de mail entre géomaticiens anglosaxon.

Pour ce changement de projection, bizarrement il va falloir repasser par un raster intermédiaire .tif en projection cible L2E. C'est nécessaire car le changement de projection sur des tuiles de format JP2000 génère des artefacts en bordure de tuiles qui provoquent une difficulté d’assemblage. Contrairement au format géotiff, le format d’hyper compression JPEG2000 ne gère pas la transparence (nodata) dans les 3 bandes rasters RGB classiques. Une solution est donc l'ajout d'une 4° bande raster qui comprend l’information de transparence nodata aux bordures des tuiles. Les applicatifs qui utilisent ces tuiles doivent donc gérer la transparence via la 4° bande raster.  

2.1 en utilisant GDAL et son driver intégré JP2OpenJPEG

2 étapes :

2.1.1 Re-projection dans un autre format L93.jp2  vers L2E.tif :

gdalwarp -overwrite -s_srs EPSG:2154 -t_srs EPSG:27572 -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 -co TILED=YES -r average -dstnodata "255 255 255" -dstalpha Raster_L93.jp2 Raster_L2E.tif

Cette opération se passe bien ; les Raster_L2E.tif  s’assemblent bien pour tous les niveaux de zoom.

Note : la projection  EPSG:27572 est mal gérée dans QGIS/GDAL avec un décalage de quelques mètres en superposition de données vecteurs de relevé terrain en projections L2E. Il est  donc préférable de  spécifier une projection IGNF qui correspond à un "vrai" L2E
gdalwarp -overwrite -s_srs EPSG:2154 -t_srs IGNF:LAMB2C -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 -co TILED=YES -r average -dstnodata "255 255 255" -dstalpha Raster_L93.jp2 Raster_L2E.tif

Notez que nous ajoutons une 4° bande raster alpha pour la gestion des pixels en transparence (-dstalpha).

2.1.2 Rranscription de format:  L2E.tif à L2E.jp2 :

gdal_translate -a_srs EPSG:27572 -r average  -of JP2OpenJPEG Raster_L2E.tif Raster_L2E.jp2

Cette opération convertit bien le format mais nous obtenons cependant de légers artefacts entre les tuiles pour les petites échelles de zoom ;   un liseré blanc s’affiche entre les tuiles. Pour les grandes échelles (>1/5000°), le liseré disparaît… La compression du format jp2 est de mauvaise qualité pour les petites échelles.  A date, je n'ai pas trouvé de solution "simple" pour résoudre ce problème pour une conversion dalle à dalle. Je renvoie vers le point 3/ pour une solution où l’accolage des dalles re projetée est propre.

Avec QGIS on vérifie le résultat en précisant bien la propriété de transparence dans la 4° bande raster :

Ouverture pas du tout OK des images tuilées avec QGIS sans prendre en compte la transparence dans la 4° bande 

Ouverture presque OK des images tuilées avec QGIS avec  paramétrage de la transparence dans la 4° bande raster (avec léger liseré entre les tuiles pour les petites échelles de zoom  surligné en rouge)

2.2/ Une variante avec le driver Kakadu

Si vous souhaitez pouvoir générer des dalles jp2 avec une très forte compression (pour usage web intensif des dalles), Kakadu peut être utile sinon passez ce paragraphe.

Le driver KAKADU autorise de multiples options de compression, cependant aucune ne résout le problème des liserés entre tuiles. Attention aussi à créer un tif sans compression (ne pas passer la commande GDAL_translate) car l'opération suivante kdu_compress ne reconnait que des GeoTiff non compressés.

GdalWarp va faire le travail via l'option -dstnodata

c:

cd c:\Program Files\GDAL

gdalwarp -overwrite -s_srs EPSG:2154 -t_srs IGNF:LAMB2C -of GTiff -co TILED=YES -r average -dstnodata 0 Raster_L93.jp2 Raster_L2E.tif

Pour la conversion du TIF en JP2 par Kakadu, on peut appliquer les paramètres recommandés par l'IGN en p. 21 de ce document. Cela donne :

c:

cd c:\Program Files (x86)\Kakadu

kdu_compress -i Raster_L2E.tif -o .Raster_L2E.jp2 -rate - Clevels=8 Clayers=12 Creversible=yes Cprecincts={256,256},{256,256},{128,128} Corder=RPCL Sprofile=PROFILE1

ou encore en version avec haute compression (Taille divisée par 10)

kdu_compress -i Raster_L2E.tif -o .Raster_L2E.jp2 -rate 1.2 Clevels=8 Clayers=12 Cycc=yes Creversible=no Cprecincts={256,256},{256,256},{128,128} Corder=RPCL ORGgen_plt=yes ORGtparts=R Sprofile=PROFILE1

3/  Reprojection en bloc à l’aide d’une mosaïque des dalles sources

Cette solution produit les résultats les plus convaincants et est à privilégier. On utilise uniquement des commandes GDAL que l’on peut automatiser via SAS ou Python. Il s'agit d'assembler les dalles sources pour créer une mosaïque qui sera re-projetée d’un bloc et redécoupée selon le modèle des dalles sources. On évite ainsi tous les artefacts de jointure entre les tuiles.  Attention cette solution est fortement consommatrice de ressources machines. Pour la  re-projection des dalles ortho ou parcellaires de l’IGN, on boucle le traitement en quatre étapes suivant pour chaque agglo ou au maximum par département.

3.1/ Transcription (avec compression) de toutes les sources L93.jp2 en L93.tif

Cette transcription passe par un format tif intermédiaire,  par n commandes (une par tuile source) du type :

gdal_translate  -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 -of GTiff  Raster_L93.jp2 Raster_L93.tif

3.2/ Assemblage des dalles L93.tif sous forme d'une mosaïque virtuelle de dalles :

gdalbuildvrt -hidenodata -addalpha - Index_L93TIF.vrt *.tif

Si votre rectangle assemblé de tuiles n'est pas complet,les paramètres -hidenodata et -addalpha permettent de gérer la transparence des dalles manquantes. Aux trois bandes classiques du raster (une pour chaque couleur Red, Green, Blue), nous ajoutons une 4° bande "alpha" dont la valeur est égale à 255 pour les pixels colorés et 0 pour les pixels transparents. Toutes les tuiles créées à partir de cet index possèdent alors cette information de "nodata" sur la 4° bande facilement réinterprétable par les SIG pour la gestion de la transparence.

3.3 / Re-projection de la mosaïque virtuelle en projection L2E pour sortir une grosse image tif.

Nous allons construire une grosse tuile facile à découper. Pour cela, il faut connaitre la taille du rectangle de tuiles : nombre de lignes et de colonnes. Pour l’agglomération de Roanne par exemple, mon échantillon de la base ortho de l'IGN comprend 332 tuiles chacune de 2000 pixels *2000 pixels réparties sur 19 lignes et 28 colonnes, soit une mosaïque de 532 tuiles potentiels. Le rectangle mosaïque comprend donc des parties vides où nous n’avons pas de sources. Les dalles seront vierges et pourront être éliminées a posteriori. Avec la commande gdalwarp, on force la taille de la mosaïque à 19*2000=38000 pixels en ligne  et 28*2000=56000 pixels en colonne. Du fait de la re-projection, la résolution de la dalle de sortie est nécessairement légèrement différente de celle des dalles d’entrée. Il faut prévoir des ressources sur un disque ; pour l’agglomération de Roanne, la mosaïque de sortie pèse 2,6 GO. Les temps de traitement sont OK et permettent de générer plusieurs dizaines de milliers de dalles si nécessaire.

gdalwarp -overwrite -s_srs EPSG:2154 -t_srs IGNF:LAMB2C -of Gtiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9  -co TILED=YES -ts 56000 38000 Index_L93TIF.vrt  Mosaique_L2E.tif

Note : en automatisation, la récupération des informations sur la taille des mosaïques passe en amont par la commande gdalinfo. On récupère la taille de la mosaïque virtuelle pour alimenter la mosaïque physique. On peut aussi vérifier que tout se passe bien pour mosaïque physique en sortie.

Gdalinfo Index_L93TIF.vrt >info_ Index.txt

Gdalinfo Mosaique_L2E.tif >info_Mosaique.txt

3.4/ Conversion de format .tif vers .jp2 et  redécoupe simultanée de l’image mosaïque par dalle

En fait, nous devons générer n commandes translate en bouclant sur la largeur et la longueur de la tuile par pas de 2000 pixels.

Avec une succession de commandes du type :

gdal_translate  -of JP2OpenJPEG -srcwin 0 0 2000 2000 Mosaique_L2E.tif Dalle_1_1.jp2

gdal_translate  -of JP2OpenJPEG -srcwin 2000 0 2000 2000 Mosaique_L2E.tif Dalle_2_1.jp2

gdal_translate  -of JP2OpenJPEG -srcwin 0 2000 2000 2000 Mosaique_L2E.tif Dalle_1_2.jp2

gdal_translate  -of JP2OpenJPEG -srcwin 2000 2000 2000 2000 Mosaique_L2E.tif Dalle_2_2.jp2

….

Pour cette partie (/4), j’indique un exemple de code SAS qui interprète l'indexation ligne / colonne sur mes sources de l'IGN et qui boucle pour régénérer et redécouper toutes ces tuiles selon ce même système d'indexation et de nommage. Il faut bien mettre au point ce paramétrage de bouclage pour que l’on ait une bonne correspondance entre les tuiles sources et le résultat.

/*** génération de la découpe de l'image mosaïque  ***/

/*** liste des noms des dalles sources pour récupération du mode d'indexatio, je lis en entrée des dalles nommée selon le schéma suivant

42-2013-0765-6550-LA93-0M50-E100.jp2

**/

 data WORK.liste    ;

       infile 'C:\Mosaic\JP2In_Roanne\liste.txt'  ;

         format l_A c_A $4. nom $100. l c 4.;

       input

                 enrg $72. ;

       ;

         nom=substr(enrg,37,36);

         l_A=substr(enrg,45,4);

         c_A=substr(enrg,50,4);

         l=input(l_A,4.);

         c=input(c_A,4.);

     drop enrg;

     run;

proc means data=liste min max;

var l c;

output out=minimax;

run;

Data minimax; set minimax;

if _Stat_="MIN"  then do; CALL SYMPUT('Minl', l);CALL SYMPUT('Minc', c);end;

else if _Stat_="MAX" then do; CALL SYMPUT('Maxl', l);CALL SYMPUT('Maxc', c);end;

run;

/*** boucle pour générer les commandes de redécoupe de la mosaïque physique ***/

%macro nomf;

CGDAL="C:"; output;

CGDAL="C:\Mosaic\JP2In_Roanne"; output;

%do i=&Minl %to &Maxl ;

%do j=&Minc %to &Maxc ;

%let ib=%eval(&i-&Minl);%let ib2=%eval((&i-&Minl)*2000);

%let jb=%eval(&j-&Minc);%let jb2=%eval((&j-&Minc)*2000);

/***42-2013-0765-6550-LA93-0M50-E100.jp2**/

CGDAL="gdal_translate -of JP2OpenJPEG -srcwin "!!tranwrd(right(put(&ib2,$4.)),' ','0')!!" "!!tranwrd(right(put(&jb2,$4.)),' ','0')!!" 2000 2000 Mosaique_L2E2.tif L2Ejp2\42-2013-"!!"0"!!trim(left(&i))!!"-"!!trim(left(&j))!!"-LA2E-0M50-E100.jp2";

output;

%end;

%end;

/*** on ote les dalles vierges (DOS) **/

CGDAL="pushd C:\Mosaic\JP2In_Roanne\L2Ejp2";output;

CGDAL="for %%F in (*.jp2) do (if %%~zF lss 5000 del %%F)";output;

CGDAL="gdalbuildvrt  mosaique_L2Ejp2.vrt *.jp2";output;

CGDAL="popd";output;

CGDAL="pause";output;

run;

%mend;

data bat  ;

length CGDAL $800.;

%nomf;

run;

run;

/*** écriture du fichier de commande .bat à exécuter sous dos ***/

 data _null_;

     file 'C:\Mosaic\JP2In_Roanne\CGDAL_DecoupTranslateMosaique.bat' lrecl=800;

     set  WORK.Bat   end=EFIEOD;

        format CGDal $800. ;

        put CGDal $ ;

      run;

Bravo à ceux, certainement nombreux, qui sont arrivés jusque-là. Les autres survivront cependant très bien sans avoir lu ce post !