Capturer des ressources WMS dans un fichier raster avec GDAL🔗

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.

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…

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ètres 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 tâtonnements 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 obtenir 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.

À 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 tâtonner, j'ai trouvé que la valeur limite pour la couche scan25 du service WMS que j'interroge est de 2,64.

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 cœur 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.

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.

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

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…