Comment accéder en masse aux informations vectorielles d'Open Street Map ?

Open Street Map (OSM) est une carte ouverte à couverture mondiale. Chacun peut puiser dans ses rendus cartographiques, ses outils pour développeurs (API) et/ou la base de données sous-jacente. C’est une source de données localisées de mieux en mieux alimentée, précise et « gratuite ». En France en particulier, soutenu par l’association OSM France et un secteur public très actif sur le développement de cette carte, le contenu est très régulièrement enrichi notamment par les données cadastrales de la DG FIP.

La structure de la base de données sous-jacente est cependant un peu complexe et nécessite un investissement chronophage pour réussir ses extractions. Il faut quelques connaissances en géomatique pour plonger dans le bain de données OSM. Le format du fichier .osm, xml tout en un, est déroutant pour les utilisateurs de SIG et les statisticiens. Le fichier .osm est un format texte balisé (xml) facile à manipuler par les développeurs. Mais il l’est beaucoup moins par les analystes. Les informations ne sont pas tabulées en colonne et les objets géographiques sont, le plus souvent, mal lus par les SIG classiques. QGIS réussit à lire nativement les fichers .osm. Mais il est impossible d’espérer lire et requêter sur un fichier France.osm directement : les structure et volumétrie de données (80 GO pour la France) ne passent pas…   

Ce post détaille une méthode d’extraction des données vectorielles d’OSM et de conversion sous formes de tables couches vectorielles « classiques » facilement interprétables avec un système d’information géographique. On se concentre ici sur les extractions à forte volumétrie.

Il existe de nombreux moyens d’accès aux informations d’OSM.  Les plus faciles et populaires sont :

• L’API Overpass ;  avec une interface web très facile d’emploi, elle permet de faire des extractions ciblées sans connaissance de la structure de la base de données OSM.
• Les shapes produits par Géofabrik sur tous les pays du monde et accessibles en suivant ce lien ou encore sur les grandes aires métropolitaines par Mapzen.

Mais après quelques manipulations de ces outils, nous touchons leurs limites.

Les capacités de l’Api overpass sont vite dépassées lorsque l’on souhaite extraire sur un large territoire. Les transactions qui passent par le web limitent la volumétrie d’extraction. Il est par ailleurs très fastidieux de monter des extractions multi-thématiques et/ou sur plusieurs tag osm avec cette API.

Les shapes de Géofabrik ou Mapzen sont très utiles mais trop généralistes, parfois incomplets (manque d’objets ou d’attributs) et on ne contrôle pas leur contenu. Par ailleurs, les limites de volumétrie des shapes entraînent un découpage régional pour la France et il faut construire des requêtes complémentaires pour obtenir une couche propre pour tout le pays.  

Bref, nous souhaitons passer au stade supérieur ; être autonome dans l’accès à la source, pouvoir gérer de grosses volumétries et mieux comprendre ce que l’on extrait.

1/ Pré-requis :

Je mentionnais plus haut que l’accès à OSM est entre guillemet « gratuit ». La carte est libre mais l’accès à la base de données n’est pas vraiment « grand public ». Il y a quelques pré-requis techniques, le plus souvent documentés en anglais. Il est souhaitable d’avoir une expérience de la géomatique et des SIG pour s’attaquer à ces extractions. On travaille sous windows sans interface logiciel, le langage de commandes MSDOS[1] ne doit pas vous rebuter et la liste des pré-requis est donc un peu longue :

1.1 Il faut d’abord lire et comprendre l’organisation des objets cartographiques de la base de données OSM : le wiki de référence en anglais sur les OSM Map features détaille la nomenclature par catégories (Amenity, Shop, Building, Boundary…)  et sous-catégories des objets de la carte ainsi que les propriétés complémentaires qui peuvent être associées à ces objets (Name, addresses, annotation…)

1.2 Installer le logiciel OSMosis en suivant la procédure décrite dans ce wiki : http://wiki.openstreetmap.org/wiki/Osmosis/Quick_Install_(Windows).

Osmosis utilise java, il faut avoir installé ce logiciel et vérifié que son chemin d’accès est bien reconnu par Osmosis.  

Pour cela, avec votre explorateur Windows, chercher « java.exe » et modifier le fichier de lancement

C:\Program Files (x86)\osmosis\bin\osmosis.bat en conséquence

Sur mon PC, la définition de variable d’accès à java s’écrit :

set JAVACMD=C:\Program Files (x86)\Java\jre1.8.0_111\bin\java

Attention il y a une petite erreur dans le wiki car il ne faut pas de guillemets “ “  pour la définition de la variable javacmd.

N’oubliez pas aussi d’ajouter dans votre variable d’environnement[2] PATH, le chemin d’accès à Osmosis, sur mon PC : "C:\Program Files (x86)\osmosis\bin"

Pour vérifier le bon fonctionnement d’Osmosis , la commande : "C:\Program Files (x86)\osmosis\bin\osmosis.bat" doit renvoyer :

1.3 Installer GDAL pour l’accès à ogr2ogr :

GDAL (Geospatial Data Abstraction Library) est 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 l’invite de commande DOS. L’installation est expliquée sur le site OSGeo4W. On peut aussi se référer à ce mode d'emploi ; attention à bien installer des versions de GDAL et de python cohérentes.

1.4 Télécharger les données brutes d’OSM en format hyper compressé .pbf sur le site de téléchargement de géofabrik. Il faut une bonne fibre Internet et un peu d’espace sur vos disques pour télécharger ces sources (Europe =18.2GB, FR 3.3GB, Germany 2.9 GB, Italie 1.2GB, GB 0.9…). Ces fichiers organisés par continent et pays comprennent l’exhaustivité de la base OSM. Elles sont actualisées tous les jours par Géofabrik.

1.5 Lire aussi l’excellent article en français sur l’extraction de données .osm à l’aide du seul logiciel org2ogr : vraiment utile et très clair. Je suis parti de cette méthode que j’ai adaptée pour l’extraction de gros volumes en ajoutant une étape osmosis en amont de l’étape ogr2ogr.

1.6 Lors de la mise au point et de la validation de mes extractions, je vais toujours sur la carte en ligne OSM dans des zones que je connais bien. Avec l’outil ? j’interroge la carte à un endroit donné et j’ai ainsi accès  à tous les objets avoisinants et tous leurs attributs. En cliquant sur un objet, on voit tous les tags associés. C’est un très bon moyen pour comprendre la structure de cette base.

La description d'un objet polygone OSM

ZAC de Kerfontaine (Pluneret 56) : tags base OSM

2/ Conversion des données osm en couches vectorielles :

Le format texte en clair (xml) .osm est très consommateur d’espace disque. 

La méthode que je préconise part des fichiers ultra compressés .pbf France.pbf =3.3 GB versus 80 GB pour la version décompressée .osm). Osmosis lit directement le fichier compressé et permet de faire des filtres pour sortir un fichier .osm manipulable qui pourra facilement être lu par les SIG ou mieux converti en format géo-tabulaire pour SIG (Shp ou MapInfo files par exemple).

Les objets géographiques OSM sont de types point, ligne polygone ou relation. On oublie les relations, car nous n’en avons a priori pas d’usage avec un SIG.

2.1/ Extraire des points : création d’un fichier France de tous les points d’intérêt de catégorie Shops (magasins) et Amenities (infrastructures)

Avec la commande CD, placez-vous d’abord dans le répertoire  de vos sources .pbf.  L’extraction est réalisée en deux étapes.

2.1.1 Création d’un fichier .osm

osmosis --read-pbf France.pbf --tf accept-nodes shop=* amenity=* --tf reject-relations --used-way --write-xml  FR-ShopsAmenities.osm

Temps de traitement sur mon PC portable 30 mn

Explication des options

--read-pbf option osmosis de lecture d’un fichier osm en format hyper compressé pbf

--tf accept-nodes shop=* amenity=* : tf est l’abréviation de tag Filter, filter tous les nœuds (ie points) des catégorie shop et amenity. Le jocker * signifie que l’on prend toutes les sous-catégories des deux catégories.  

--tf reject-relations: on exclut les relations, à préciser dans toutes nos extractions osmosis sinon on se retrouve avec tous les objets spécifiés par ces relations.

--used-way : Ne sélectionner que les points de polygones identifiés idem à préciser dans toutes les extractions points sinon on extrait des points en pagaille qui ne correspondent pas à notre filtre.
--write-xml  FR-ShopsAmenities.osm : on écrit la sortie dans le répertoire courant un fichier xml .osm

La liste des données (variables tags) extraite par osmosis est définie dans un fichier osmconf.ini. Faites une recherche de ce fichier sur votre PC (chez moi sur C:\Program Files (x86)\GDAL\gdal-data)  éditez-le par exemple avec Notepad+[3]. Pour chaque type d’objet, le fichier liste les attributs extraits par défaut par osmosis. Ajoutez et supprimez des attributs à votre guise en fonction de vos souhaits d’extractions.

Voici ma configuration d’extraction d’attributs pour des objets du type point :

2.1.2 Transformation en fichier SIG

ogr2ogr -f "ESRI Shapefile" FR_shopAmenities_Points.shp -skipfailures -overwrite -gt 65536 -dialect SQLITE  -sql "SELECT  GEOMETRY,osm_id,name, coalesce(shop, amenity) as type, address,brand, operator,other_tags from points" FR-ShopsAmenities.osm -lco ENCODING=UTF-8 FR-ShopsAmenities.osm

-f "ESRI Shapefile" FR_shopAmenities_Points.shp : création d’un fichier shape. Je privilégie le format Shp car ogr2ogr ne gère bien l’encodage des caractères UTF8 via l’option -lco ENCODING=UTF-8 FR-ShopsAmenities.osm

-overwrite : écrasement du précédent fichier de même nom s’il existe

-gt 65536 : option à utiliser pour meilleure gestion de la puissance machine, ne semble pas indispensable cependant

-dialect SQLITE : précise que l’on va utiliser du sql avancé dans la requête qui suit

-sql "SELECT  GEOMETRY,osm_id,name, coalesce(shop, amenity) as type, address,brand, operator,other_tags from points" FR-ShopsAmenities.osm : on passe une requête sql sur la base .osm préalablement constituée via osmosis. La requête figure entre guillemets et précise  la liste des champs à extraire. Ces champs doivent avoir été préextraits en amont via osmosis et doivent donc être listés dans le fichier de configuration osmconf.ini. Notez que tous les tags non spécifiés dans le fichier de configuration sont accessibles en désordre dans un champs other_tags. coalesce(shop, amenity) as type utilise une fonction sql avancée coalesce et indique que l’on inscrit dans le champ type la valeur de la sous-catégorie shop si celle-ci est renseignée et la valeur de la sous-catégorie amenity sinon.

On peut alors charger et manipuler le résultat avec son SIG préféré.

  

2.2 Extraction d’un fichier des contours des zones industrielles et commerciales :
La démarche est très similaire à celle des points sauf que l’on filtre ici des polygones

2.2.1 Création des polygons .osm
osmosis --read-pbf france.pbf --tf accept-ways landuse=retail,commercial,industrial,port shop=department_store,mall,supermarket building=retail,industrial,commercial,warehouse --tf reject-relations --used-node  --write-xml  FR-LanduseWays3.osm

On filtre ici les polygones des sous-categories retail,commercial,industrial,port du tag landuse + les trois sous-catégories department_store,mall,supermarket de la categorie shop + quelques tags de la catégorie building

--used-node  bizarrement il faut préciser que l’on ne sélectionne que les nœuds des polygones de la sélection.  

2.2.2 Définition d’un shape des polygons

ogr2ogr -skipfailures -overwrite -gt 65536 FR_Retail_Poly.shp -sql "SELECT osm_way_id, name, landuse,building,shop,brand,operator from multipolygons" FR-LanduseWays3.osm -lco ENCODING=UTF-8 FR-LanduseWays3.osm

A titre indicatif, voici la trace msdos de ce traitement avec le temps d’extraction :

Et les fichiers résultats :

2.3 Extraire tout le bâti français :

Attention c’est « du lourd », pour la France métropolitaine, le tag « building » (bâti) couvre près des 2/3 de la base soit 55GO sur les 80GO du fichier .osm de référence. Avec cette extraction je souhaite donc « stress tester » la méthode.

osmosis --read-pbf France.pbf --tf accept-ways building=* --tf reject-relations --used-node  --write-xml FR-building.osm  

Log MSDOS (un peu plus de 3h de traitement)

Taille du fichier de sortie :

Puis :

ogr2ogr -f "SQlite" FR_Building_Poly.SQLite -skipfailures -overwrite -gt 65536 -dialect SQLITE  -sql "SELECT GEOMETRY,osm_way_id,name,building from multipolygons " FR-building.osm

Note : cette requête crée un shape énorme (8GO) difficilement exploitable. Il est préférable de se concentrer sur la création de couches de bâti thématiques, ici le bâti qualifié :

ogr2ogr -f "ESRI Shapefile" FR_BuildingQualif_Poly.shp -skipfailures -overwrite -gt 65536 -dialect SQLITE  -sql "SELECT GEOMETRY,osm_way_id,name,buiding,address from multipolygons  where building is not null" FR-building.osm -lco ENCODING=UTF-8 FR-building.osm

2.4 Tous les Parkings (points et polygones)

On extrait séparément la couche de points et la couche de polygones

osmosis --read-pbf France.pbf --tf accept-nodes amenity=parking,motorcycle_parking,parking_entrance,bicycle_parking,bicycle_rental,car_sharing,garages,parking_space tourism=caravan_site  site=parking building=garages,garages landuse=garages --tf reject-relations  --used-way --write-xml  FR_Parking_nodes.osm

osmosis --read-pbf France.pbf --tf accept-ways amenity=parking,motorcycle_parking,parking_entrance,bicycle_parking,bicycle_rental,car_sharing,garages,parking_space tourism=caravan_site  site=parking building=garages,garages landuse=garages --tf reject-relations --used-node  --write-xml  FR_Parking_ways.osm

ogr2ogr -f "ESRI Shapefile" FR_Parking_nodes.shp -skipfailures -overwrite -gt 65536 -dialect SQLITE  -sql "SELECT  GEOMETRY,osm_id,name, coalesce(amenity, tourism, site, building) as type, address,brand, operator,other_tags from points" FR_Parking_nodes.osm -lco ENCODING=UTF-8 FR- FR_Parking_nodes.osm

ogr2ogr -f "ESRI Shapefile" FR_Parking_ways.shp -skipfailures -overwrite -gt 65536 -dialect SQLITE  -sql "SELECT  GEOMETRY,osm_way_id,name, coalesce(amenity, tourism, site, building) as type, address,brand, operator,other_tags from multipolygons" FR_Parking_ways.osm -lco ENCODING=UTF-8 FR- FR_Parking_ways.osm

2.5 Toutes les piscines et terrains de sport :

Attention, les commandes OSMosis et Ogr2Ogr sont parfois un peu capricieuses. Je note une bizarrerie sur certaines requêtes ogr2ogr avec requête option -sql utilisant des fonctions avancées accessibles via l’option -dialect SQLITE. Ceci ne fonctionne pas pour toutes les extractions osmosis (l’ordre SQL  ne reconnait plus la table multipolygons). Dans ce cas il faut simplifier la requête sql en supprimant l’option -dialect SQLIT. Supprimez aussi la mention de l’extraction de l’objet qui est implicite en dehors de SQLite et les appels de la fonction coalesce.

osmosis --read-pbf France.pbf --tf accept-nodes leisure=swimming_area,swimming_pool,water_park,pitch,beach_resort natural=beach  sport=swimming,tennis --tf reject-relations  --used-way --write-xml  FR_Leisure_nodes.osm

osmosis --read-pbf France.pbf --tf accept-ways leisure=swimming_area,swimming_pool,water_park,pitch,beach_resort natural=beach  sport=swimming,tennis --tf reject-relations --write-xml  FR_Leisure_ways.osm

ogr2ogr -f "ESRI Shapefile" FR_leisure_nodes.shp -skipfailures -overwrite -gt 65536 -dialect SQLITE  -sql "SELECT GEOMETRY,osm_id,name,coalesce(leisure,natural,sport) as type,address,access,brand,operator,other_tags from points" FR_leisure_nodes.osm -lco ENCODING=UTF-8 FR_leisure_nodes.osm

ogr2ogr -f "ESRI Shapefile" FR_leisure_ways.shp -overwrite -sql "SELECT osm_way_id,name,leisure,natural,sport,address,access,brand,operator,other_tags FROM multipolygons" FR_Leisure_ways.osm -lco ENCODING=UTF-8 FR_Leisure_ways.osm

En filtrant les sorties sur access= « private », on obtient un fichier de localisation des piscines et cours de tennis privés. La source de ce type de bâti est la DGI n’est cependant pas totalement exhaustive sur ces signes extérieurs de richesse : cf le premier post historique de ce blog piscines à la carte.

Les extractions sur le linéaire routier fonctionnent sur le même principe. L’usage des linaires routiers dans les graphes de calculs de distanciers/itinéraires méritent néanmoins un développement spécifique que j’espère pouvoir poster dans un prochain épisode.  

Si vous êtes arrivé jusqu’ici, je vous tire mon chapeau en guise de conclusion. Vous êtes désormais mûrs pour puiser ce que vous cherchez dans OSM et nous faire vos retours de trucs et astuces d’usage.

---

[1] Accessible par l’invite de commande de windows

[2] Accessible sur Panneau de configuration\Tous les Panneaux de configuration\Système puis propriétés système+variables d’environnement + variable système + modifier la variable PATH et y ajouter le chemin d’accès à votre Osmosis.bat
[3] Si vous avez plusieurs installations Gdal (intégrées avec Qgis par exemple), il y a plusieurs fichiers osmconf.ini il n’est pas évident de savoir quelle version vous utilisez (voir votre variable d’environnement PATH). Dans ce cas, vous pouvez modifier l’un de vos fichiers osmconf.ini et recopier le en annule et remplace de tous les autres. 

Des cartes pour Revivre dans le Monde

Chez Pitney Bowes, nous sommes intervenus en mécénat de compétence pour la mise en place d'un outil de cartographie pour l'association Revivre dans le Monde.

Ce réseau d'associations fondé par des anciens de Danone vient en aide aux exclus, par la collecte et la distribution de produits de première nécessité. Il propose aussi via cette activité de recyclage de ces produits, de venir en aide aux exclus du travail, par la réinsertion.

Revivre dans le monde organise une chaîne de collecte et de recyclage de produits alimentaires d'hygiène et d'entretien auprès de partenaires industriels producteurs et d'agriculteurs volontaires. Les produits donnés par ces partenaires sont stockés et redistribués à des associations d'épicerie sociale et solidaire, des associations humanitaires et caritatives, des centre communaux d'action sociale, des associations d'urgence ainsi que des centres d'accueil et d'hébergement.  Il s'agit aussi de lutter contre le gaspillage agro-alimentaire par la collecte et distribution d'invendus. Des emplois d'approvisionnement, de manutention/stockage ou de logistique sont créés par le réseau et proposés à des chômeurs de longue durée.

En région Île de France et en Auvergne Rhône Alpes, Revivre dans le Monde a monté un projet "tournée village". Un camion distribue des produits de première nécessité à des populations  rurales ou péri-urbaines en grande précarité et éloignées des épiceries sociales. Ces tournées villages sont organisées en lien avec les communes. Le département de l'Essonne et la région Auvergne Rhône Alpes sont pilotes pour ce projet.

L’outil développé par Pitney Bowes "Précari-map" pour Revivre dans le Monde permet de cibler précisément les lieux des actions en cours de ces Tournées Villages dans les deux régions. Il apporte  aussi des réponses étayées à des appels à projets visant ces deux régions.

Voici un lien vers cette cartographie en île de France :

Cliquez sur l'image pour accéder à la carte interactive

Voila une action sociale concrète où un support simple de cartographie interactive est très utile.
Bravo à Blandine Jezequel et aux responsables de Revivre dans le Monde pour ce projet "pro bono"  rondement mené.

Chers Lectrices, Lecteurs, si vous avez des contacts associatifs ou avec des ONG qui ont des projets en lien avec l'information géographique, n'hésitez pas à nous en parler.

Open Street Map : toujours plus précis

Encore un post pour tresser des lauriers au projet Open Street Map (OSM). Depuis 2004, OSM est une collecte et diffusion communautaire et libre de cartographie vectorielle du monde.

Avec l'aide de l'excellent Fréderik Ramm de géofabrik, nous suivons l'évolution de l'alimentation d'OSM.

En l'espace de 2 ans (2014 vs 2016), le réseau routier mondial est passé d'une couverture de  38.9 Millions de kilomètres à 47.7M soit +23%. La croissance est encore plus importante pour les autres linéaires (Rivières, trains,...), pour les différentes zonages de couverture terrestre (Forêt, prairies, lac....), ainsi que pour les points d'intérêt.


OSM : Evolution de la couverture kilométrique du routier et autres linéaires et du nombre d'objets de la cartographie générale

Les utilisateurs d'OSM bénéficient partout dans le monde de ces améliorations. En France par exemple, nous avons noté une très nette amélioration de  la localisation du commerce en centre ville presque comparable désormais à celle de google.

Rue commerçante du centre de Paris
OSM

GoogleMap

Cartographier la planète : la route est longue, mais nous roulons toujours plus nombreux pour OSM.

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 !