Capturer des ressources WMS dans un fichier raster avec GDAL

Introduction

Le protocole WMS (Web Mapping Service) est un protocole officiel de l'OGC qui permet de télécharger de l'information géographique directement sous forme d'images. Ça a des avantages (les ressources sont en ligne et on a rien à stocker sur sa machine) et également des inconvénients dont l'impossibilité d'utiliser le fichier bitmap pour faire des tirages papiers de grande taille (A0) ou encore, des délais de téléchargement qui peuvent être parfois longs et peu réactifs suivant l'étendue sur laquelle on travaille.

Pour pallier à ce problème, il est parfois bon de vouloir récupérer la dalle raster (bitmap) sur une machine et de travailler directement avec. Toutefois, si vous n'avez pas accès directement aux fichiers, vous pouvez demander au service WMS de faire ce travail pour vous. Après tout, il est fait pour produire des images !

Cet article présente une méthode pour récupérer un bon gros raster à partir d'un service WMS. J'ai utilisé les ressources de mon travail pour tester cette méthode. Elle n'a pas la prétention d'être la solution de référence mais elle propose tout ce que j'ai pu trouver sur le sujet qui donne un résultat final digne de ce nom.

Bien entendu, nous allons utiliser des logiciels libres ! Il y a sans doute des logiciels privateurs qui font la même chose en mieux mais comme je n'y ai pas accès pour des raisons de coût et de temps, je me base exclusivement sur des logiciels disponibles dans la distribution Debian.

Dans le monde du SIG Libre, la référence pour gérer des fichiers image, c'est GDAL. C'est la bibliothèque (et les outils qui vont avec) que nous allons utiliser dans le cours de cet article.

Ainsi, avant de commencer, je vous invite à lire la documentation de GDAL sur les services WMS.

Nous étudierons comment configurer un accès vers un service WMS. Puis nous verrons quel est l'élément de résolution important à gérer. Enfin, nous effectuerons un exemple de récupération de dalle à l'aide d'un outil fourni par GDAL.

Configuration de l'accès au service WMS

Avant de commencer, il faut indiquer aux outils GDAL comment se connecter au service WMS. Pour cela, il faut disposer des bonnes URL, savoir quelles sont les couches qu'on désire récupérer, retrouver leur nom au sein du service WMS, gérer leurs projections, leurs styles, etc...

Ces paramètres sont trop nombreux pour la ligne de commande et c'est pour cela que GDAL utilise un fichier de définition au format XML, comme indiqué dans la documentation. Voici un exemple de fichier de configuration XML (sans les commentaires qui viendront après).

<GDAL_WMS>
  <Service name="WMS">
    <Version>1.3.0</Version>
    <ServerUrl>http://georef.application.i2/cartes/mapserv?</ServerUrl>
    <CRS>EPSG:2154</CRS>
    <ImageFormat>image/jpeg</ImageFormat>
    <Layers>scan25</Layers>
  </Service>
  <DataWindow>
    <UpperLeftX>1181556</UpperLeftX>
    <UpperLeftY>6174866</UpperLeftY>
    <LowerRightX>1184158</LowerRightX>
    <LowerRightY>6173064</LowerRightY>
    <SizeX>2400</SizeX>
    <SizeY>2400</SizeY>
  </DataWindow>
  <OverviewCount>12</OverviewCount>
  <BlockSizeX>512</BlockSizeX>
  <BlockSizeY>512</BlockSizeY>
  <Projection>EPSG:2154</Projection>
  <BandsCount>3</BandsCount>
  <ClampRequests>false</ClampRequests>
</GDAL_WMS>

Globalement, on voit que ce fichier est séparé en trois parties:

  • La première (Service) permet de définir l'accès au serveur ainsi que la (ou les) couche(s) qui est (sont) demandée(s).
  • La seconde (DataWindow) précise la taille du fichier final (SizeX et SizeY) ainsi que l'étendue qu'on souhaite récupérer, exprimée dans les coordonnées de la couche.
  • La troisième partie (sans nom) contient des éléments de paramétrage général de la requête.

Ce que j'ai indiqué dans le fichier sont les balises minimales à insérer avant de pouvoir lancer une requête. Bien entendu, il y a d'autres options. Vous pouvez sans problème lire la documentation spécifique pour en connaître l'exhaustivité.

Vous devez néanmoins faire attention à un détail particulier, celui qui concerne les paramètres dans la partie DataWindow, comme expliqué ci-dessous...

Quelques problèmes à régler

Un des principaux problèmes à gérer réside dans "le seuil de zoom". En effet, le serveur WMS ne retourne une image que si elle respecte une certaine résolution. Par exemple, sur un scan25, il est inutile de demander une échelle inférieure à 10000: c'est le seuil minimal, au dessous, il faut fabriquer des pixels artificiels. Ça n'a aucun intérêt pour le serveur et pour le client.

C'est pourquoi, la majorité des services WMS ne servent des images qu'à condition qu'elles respectent une certaine étendue et une certaine échelle. Tout dépend de la configuration du serveur et ce paramétrage est lié à chaque couche servie. Il faut donc le tester à l'avance.

En effet, demander une image à un service WMS avec une résolution trop faible donnera un raster flou (vu d'en haut) alors que le demander avec une résolution trop importante retournera une image vide ! Il faut donc trouver la bonne valeur: celle qui permet le maximum de précision et qui se trouve donc juste au dessus du seuil qui retourne des images vides.

Dans GDAL, l'information la plus simple à gérer à ce niveau est le "Pixel Size". C'est ce que vous allez devoir découvrir pour télécharger votre dalle. Cette taille de pixels est dépendante des éléments de configuration de votre fichier XML, notamment, dans la partie . En fonction des paramètre d'étendue (les UpperLeft et LowerRight surtout) et de taille de dalle (SizeX et SizeY), la taille de pixel va varier.

Il faut déjà voir comment accéder à cette information qui doit être calculée puisqu'elle n'est pas un paramètre du fichier XML. Nous allons utiliser l'outil gdalinfo qui est écrit pour donner le plus d'informations sur un fichier raster.

$ gdalinfo wms.xml
Driver: WMS/OGC Web Map Service
Files: georef.xml
Size is 21300, 26000
Coordinate System is:
PROJCS["RGF93 / Lambert-93",
    GEOGCS["RGF93",
    DATUM["Reseau_Geodesique_Francais_1993",
            SPHEROID["GRS 1980",6378137,298.257222101,
         AUTHORITY["EPSG","7019"]],
    TOWGS84[0,0,0,0,0,0,0],
    AUTHORITY["EPSG","6171"]],
    PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4171"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",49],
PARAMETER["standard_parallel_2",44],
PARAMETER["latitude_of_origin",46.5],
PARAMETER["central_meridian",3],
PARAMETER["false_easting",700000],
PARAMETER["false_northing",6600000],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","2154"]]
Origin = (1171625.395620000082999,6182216.030849999748170)
Pixel Size = (2.649430311267601,-2.315194487692287)
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 1171625.396, 6182216.031) (  8d44'33.86"E, 42d35' 6.91"N)
Lower Left  ( 1171625.396, 6122020.974) (  8d41'24.05"E, 42d 2'44.39"N)
Upper Right ( 1228058.261, 6182216.031) (  9d25'37.28"E, 42d32'46.30"N)
Lower Right ( 1228058.261, 6122020.974) (  9d22' 5.04"E, 42d 0'25.16"N)
Center      ( 1199841.828, 6152118.503) (  9d 3'25.00"E, 42d17'47.53"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 10650x13000, 5325x6500, 2663x3250
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 10650x13000, 5325x6500, 2663x3250
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 10650x13000, 5325x6500, 2663x3250

On voit que notre élément "Pixel Size" apparaît: Pixel Size = (2.649430311267601,-2.315194487692287)

Dans notre cas, il vaut 2,64 et -2.31.

Vous allez donc devoir procéder par tatonnements pour trouver la bonne valeur de Pixel Size, c'est à dire, la valeur à partir de laquelle le serveur WMS retourne une image non vide.

Pour y parvenir, sachant que l'étendue est fixe (on veut bien récupérer la couche entre tels et tels points), c'est sur les paramètres SizeX et SizeY qu'il faut jouer. Vous devrez donc modifier votre fichier XML pour diminuer ou augmenter ces valeurs, lancer gdalinfo pour otbenir la taille de la valeur Pixel Size puis lancer une commande de récupération de dalle pour voir si, au final, vous obtenez une dalle vide ou non. Sachant que le temps de récupération est parfois long, je vous conseille de travailler sur une petite étendue pour trouver les bonnes valeurs.

A noter que la valeur de Pixel Size pour les données Y sera souvent négative (c'est le cas avec le référentiel Lambert93). Seule la valeur absolue compte.

Dans la pratique, j'ai trouvé que seule la première valeur de Pixel Size compte pour le serveur. A force de tatonner, j'ai trouvé que la valeur limite pour la couche scan25 du service WMS que j'interroge est de 2,64.

Récupérer une dalle

Maintenant que nos réglages sont complets, nous pouvons utiliser les outils GDAL. Dans la gamme, gdal_translate est un programme qui est spécialisé dans le passage d'un format à un autre. C'est donc lui notre coeur de cible.

Voici quelques exemples de commandes de récupération à partir du même fichier de configuration WMS:

# Récupérer au format GeoTIFF en compression DEFLATE de niveau 9.
# Si la taille de la dalle l'exige, le format sera en BigTiff (>4Go)
gdal_translate -of GTiff -co COMPRESS=DEFLATE -co BIGTIFF=IF_NEEDED  -co ZLEVEL=9 georef.xml georef.tiff

# Récupérer au format PNG en compression maximale (niveau 9). On
# génère un fichier world pour le calage.
gdal_translate -of PNG -co WORLDFILE=YES -co ZLEVEL=9 georef.xml georef.png

Si tout se passe bien, vous devez récupérer vos fichiers de dalles correctement projetés.

Pour finir: assembler des dalles

En théorie, vous pouvez à l'aide d'une simple requête WMS, même sur une grande étendue, récupérer une dalle dans le format voulu. Toutefois, en pratique, vous n'aurez peut-être pas le temps d'attendre que toute la dalle soit constituée pour commencer à travailler. Car dans certains cas, le processus peut prendre 4 ou 5 heures !

Le plus simple est donc de préparer plusieurs fichiers XML sur des emprises différentes et d'assembler les dalles obtenues une fois les téléchargements terminés.

GDAL dispose également de bons outils pour faire ce travail. Une méthode simple consiste à réaliser un assemblage virtuel (production d'un fichier VRT qui indique comment les dalles sont assemblées) puis de transformer (via gdal_translate), ce fichier VRT vers une seule dalle. Cette méthode demande deux lignes de script et fonctionne assez bien sauf lors de problèmes de palettes de couleur.

    # Assemblage virtuel de 3 dalles dans le fichier assemblage.vrt
$ gdalbuildvrt assemblage.vrt dalle1.tiff dalle2.tiff dalle3.tif

# Transformation du fichier VRT en fichier GeoTIFF:
gdal_translate -of GTiff -co TILED=YES -co COMPRESS=DEFLATE -co BIGTIFF=IF_NEEDED -co ZLEVEL=9 assemblage.vrt georef.tiff

Au final, on récupére un seul fichier georef.tiff qui contient tout ce dont nous avons besoin. Il est certes lourd mais il doit s'ouvrir sans difficultés dans QGis. Rien ne vous empêche d'ajouter d'autres options à gdal_translate pour obtenir un fichier moins gros (options de compression) ou avec d'autres métadonnées.

Petit bonus: faire des pyramides

Avant de passer à la conclusion et maintenant que nous disposons d'une vraie grande dalle, nous voulons logiquement l'utiliser. Si vous le faîtes directement, vous allez vous rendre compte que c'est très long. Particulièrement si vous utilisez une petite échelle (1/500000 par exemple). C'est normal: pour afficher l'image à cette hauteur, il faut "scanner" une grande partie du fichier et pendant ce temps, la machine calcule !

Pour éviter ce désagrément, on peut utiliser une astuce prévue dans GDAL: les "pyramides". Ce sont des miniatures de l'image de la dalle qui peuvent être injectées directement dans la dalle et qui permettent d'afficher rapidement l'image de la dalle a des échelles réduites.

L'outil dédié à la création de ces miniatures est gdaladdo. Sa commande est assez simple: il suffit généralement de préciser combien on veut de niveau de pyramide. En règle générale, j'en place 4 (2 4 8 16) ce qui permet une excellente réactivité générale de la dalle.

Un petit exemple:

# ajout de 4 pyramides dans le fichier avec la méthode average: gdaladdo -r average dalle.tiff 2 4 8 16

Conclusion

Nous avons maintenant suffisamment de matière pour réaliser ce qu'on peut qualifier de capture WMS. La méthode que je vous ai présentée est très empirique. Dans l'idéal, il faudrait un moyen de récupérer facilement la résolution maximale de la couche demandée (taille des pixels avec la meilleure qualité possible). De plus, un moyen de calculer, avec précision, la taille en pixels du fichier raster en fonction de l'étendue demandée et de cette résolution maximal serait une bonne avancée. Je pense qu'on peut scripter ces besoins, notamment en utilisant la bibliothèque GDAL.

Toutefois, la récupération de données raster via WMS pour constituer des grandes dalles n'est pas forcément un travail récurrent: WMS est fait pour gérer de l'affichage de grandes dalles à la volée. Il vaut donc mieux investir sur de bonnes performances sur ce service que de se tourner vers du téléchargement massif.

En attendant, nous avons tous un moyen de faire avec...

Posted ven. 03 ao�t 2012 20:30:00 Tags:

Ce sera Toulouse...

Je change de job. Je travaillerai toujours pour l'Etat et toujours pour le ministère en charge de l'agriculture (ces dernières années, le nom du ministère a bien dû changer en moyenne tous les deux ans donc voilà pourquoi j'utilise cette dénomination). Néanmoins, je rejoins l'administration centrale qui a la particularité de disposer "d'antennes" situées en province, notamment, son centre informatique. Il est situé à Toulouse (Auzeville-Tolosane plus précisément) et héberge également d'autres services. Je vais donc y travailler pendant au moins trois ans avant, éventuellement, de voir ailleurs.

Je passe donc d'un boulot multicarte et plus axé sur le métier des agents à un emploi de spécialiste technique dans les systèmes d'information. Sans être donc révolutionnaire, le changement, est bien pour maintenant !

J'ai passé une semaine début juillet pour découvrir la ville, et surtout, pour trouver un logement. Voilà mes premières impressions...

D'abord, ce que j'ai remarqué en amont de mon arrivée, c'est que Toulouse est une ville extrèmement bien décrite dans OpenStreetMap. A vrai dire, lors de mes recherches et explorations, je n'ai pas du tout utilisé GoogleMaps. De plus, grâce à OSRM, j'avais plutôt une bonne estimation de mes distances de trajet et les noms des rues sont vraiment bien saisis. Chapeau aux OSMeurs locaux qui ont donc bien travaillé.

Ensuite, les transports en commun me semblent plutôt bien organisés et bien desservis, y compris dans les communes de l'agglomération toulousaine. Je compte bien me passer de ma voiture pendant mon séjour sur place. A vrai dire, le lieu que je vais habiter est proche de toute commodités à 400m à pied de mon lieu de travail, même s'il est plus éloigné du centre-ville toulousain. En dehors de ça, le métro avec ses deux lignes est vraiment un concept intéressant. Ça reste presque à taille humaine et il est assez peu bondé. Un truc qui m'a frappé c'est la politesse des gens avec les chauffeurs de bus: tout le monde dit bonjour lors de la montée mais tout le monde salue également le chauffeur lorsqu'il descend. A Angers, ce dernier point ne se fait pas très souvent. A Toulouse, c'est systématique !

Les quelques jours passés dans la ville me donne l'image d'une cité comme étant un Paris à taille réduite. Certes, il y a du monde, de grands transports en commun, de grands ensembles mais je n'y sens pas la lassitude et la monotonie du peuple parisien.

Après quelques recherches, je constate qu'il y a des tas de choses à faire dans cette ville. Je compte bien consacrer une part non négligeable de mon temps personnel à découvrir les richesses de cette ville. D'abord, il y a un cinéma qui sort de l'ordinaire: Utopia. Ils proposent un service inédit de récupération sur clef USB de films sans DRM pour la modique somme de 5€ ! Qui dit grande ville dit également Hackerspace. C'est le cas à Toulouse avec le TetaLab. J'ai hâte de voir ce qui s'y passe, même si ça semble loin de là où j'habiterai.

D'ailleurs, ça me fait penser que ça fait des années que je cherche à prendre un abonnement ADSL auprès d'un fournisseur associatif (FDN par exemple). Je pense que ça va être possible à Toulouse grâce à l'association Tetaneutral, un FAI associatif de la fédération FDN. A voir si c'est techniquement possible !

A l'origine, j'avais programmé cette semaine de congés pour aller aux RMLL à Genève mais, malheureusement, j'ai dû la réaffecter à la dernière minute pour cette quête logement... Ce sera donc pour la prochaine fois. A noter que j'ai également complètement loupé la DebConf 11, y compris les retransmissions vidéos. Je me console toutefois grâce aux efforts de l'équipe Debian en charge de la "capture" vidéo qui a bien travaillé.

Posted jeu. 30 ao�t 2012 19:32:30 Tags: