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
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
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 !
