Dissémination de données géoréférencées – qgis-server et OpenLayers

Tout possesseur de téléphone portable devient en mesure aujourd’hui de produire de la donnée géoréférencée, qu’il s’agisse de photographies numériques ou de traces de parcours pour les données brutes, ou des données issues d’informations géoréférencées acquises. Nous avons par exemple présenté au FOSS4G-fr [1] et dans ces pages [2] la génération de modèles de bâtiments par les photographies géoréférencées prises par téléphone, et leur intégration dans QGis sur fond de carte OpenStreetMaps (OSM).

Obtenir de telles données est bien, mais les partager avec une communauté d’utilisateurs, c’est mieux. Le mode de transmission le plus pratique pour les volumes de données dont il est question est évidemment internet, ne serait-ce que pour centraliser les données partagées en vue de garantir leur cohérence si plusieurs utilisateurs les manipulent.

Trois protocoles de dissémination ont été mis en place à cet effet [3][4] : WM(T)S, WFS et WCS. Le premier dissémine des images représentant les données traitées, avec une résolution adaptée à l’échelle de la carte de l’utilisateur. Ne fournissant qu’une image des informations transmises, il n’est plus possible pour l’utilisateur de s’approprier ou retraiter lui-même des informations. Les deux autres protocoles disséminent les données brutes, vectorielles ou matricielles. Nous nous proposons d’aborder séquentiellement trois problèmes : un serveur de données pour disséminer les informations qui auront été traitées dans QGis (qgis-server) ; la récupération des informations sur un logiciel dédié tel que Qgis pour le traitement des informations ; et finalement la récupération des données dans un client Web pour affichage comme couche OpenLayers. Afin de ne pas faire durer le suspens plus longtemps, le lecteur est encouragé à consulter qgis.sequanux.org/jmfwms.html et jmfriedt.sequanux.org/reproj.html pour se faire une idée du résultat recherché.

Le contexte de notre étude porte sur le retrait de glaciers à front marin en milieu arctique. Les images des satellites Landsat sont acquises depuis les années 1970 et mises à disposition par l’USGS, par exemple sur http://landsatlook.usgs.gov/viewer.html. Elles sont plus ou moins bien géoréférencées (une petite correction est parfois nécessaire, surtout pour les plus anciennes), mais surtout souffrent d’une résolution médiocre (de l’ordre de 30 m/pixel) compte tenu des standards actuels. Plus de 40 ans d’histoire d’images satellites en font néanmoins une source irremplaçable pour appréhender l’évolution des régions qui vont nous intéresser. Notre travail initial sous QGis a été d’importer toutes ces images, les reprojeter dans un référentiel localement plan (WGS84/UTM33N, le référentiel projeté approprié pour le nord de la Norvège tel que nous en informe http://spatialreference.org/ref/epsg/wgs-84-utm-zone-33n/), s’assurer par quelques points de référence que les images sont convenablement positionnées, et dans le cas contraire les repositionner au moyen de 4 ou 5 points de référence sur le pourtour de la région considérée (presqu’île de Brøgger, Svalbard). Ces informations matricielles sont alors utilisées pour tracer manuellement (création d’un shapefile comprenant des lignes) les fronts marins des glaciers, créant ainsi un jeu de données vectorielles. Par ailleurs, un trajet en avion a été enregistré au moyen du récepteur GPS d’un téléphone portable (logiciel OSMTracker sous Android) et nous désirons insérer la séquence de points au format GPX dans nos cartes. Sous QGis, l’opération est triviale (Vector > GPS Tools > Load GPX File) mais pourrons-nous exporter cette information par le Web et l’inclure dans les pages affichées par un navigateur ? Nous verrons que la réponse n’est pas aussi triviale qu’il pourrait paraître.

1. Le serveur : qgis-server

QGis a tout prévu pour nous : un projet QGis s’exporte dans un format accessible par le Web au travers de qgis-server. Le principal exercice va donc porter sur l’installation de cet outil, puisqu’exporter le projet se résume à Project > Projet Properties et activer WMS, WFS et WCS en cochant toutes les couches (Select All) dans les deux derniers cas.

L’installation de qgis-server ne présente aucune difficulté et est décrite parfaitement sur http://docs.qgis.org/testing/en/docs/user_manual/working_with_ogc/ogc_server_support.html. On notera cependant que bien que le service WM(T)S fonctionne convenablement avec des versions plus anciennes de qgis-server, l’export des couches vectorielles (formats WFS et WCS) ne semble pas opérationnel pour les versions précédent 2.12. Nous expérimenterons donc avec une version au moins supérieure à 2.14 de QGis et de son serveur qgis-server. Sous Debian GNU/Linux, cela signifie éviter Wheezy mais travailler dans la version stable à la date de cette rédaction.

À l’issue de l’installation du serveur comme service d’Apache2 (nécessitant le support FCGI fourni par libapache2-mod-fcgid), le bon fonctionnement du service est validé (dans le cas d’une installation locale) par :

La réponse doit être du type :

Le service étant fonctionnel, divers projets QGis (fichier .qgs) sont placés, avec toutes les données associées aux diverses couches du projet, dans des sous-répertoires de /usr/lib/cgi-bin, en complétant avec une copie du serveur Fast-CGI (qgis_mapserv.fcgi) .

Si le projet QGis autorise l’export des couches au format WM(T)S ou WFS, nous validons la capacité à rapatrier ces couches (dans notre cas le projet et ses fichiers se trouvent dans /usr/lib/cgi-bin/project/landsat de qgis.sequanux.org) par :

La réponse est la liste des couches accessibles sous la forme :

Nous voyons donc que la projection dans laquelle les coordonnées des points formant la couche (ici vectorielle) est renseignée (EPSG:4326 signifie WGS84 en coordonnées sphériques tel que nous l’indique le CRS sous QGis), et l’extension de la couche (autour de 12°E et 79°N).

Ces points sont récupérés – au format GML – par une requête WFS de la forme :

La réponse est :

Nous avons ici les diverses coordonnées qui forment la ligne d’une des couches vectorielles de notre projet QGis, accessible au travers d’une requête HTML. Ces données sont donc accessibles par la communauté, il reste à en faire un bon usage avec les divers outils à notre disposition. Nous présentons quelques exemples ci-dessous.

2. Le client : Qgis

Le premier environnement le plus simple pour accéder à ces données est QGis lui même. Dans la barre des icônes de gauche de QGis, les trois icônes en forme de globe proposent l’accès respectivement aux données au format WM(T)S, WCS et WFS. Dans tous les cas, il nous faut renseigner l’url du serveur, par exemple http://qgis.sequanux.org/cgi-bin/project/landsat/qgis_mapserv.fcgi. Lors de la connexion, la liste des couches accessibles est proposée, et cliquer sur une couche suivi de Add donne accès au jeu de données. La seule subtilité dans cette application tient à ne pas utiliser un serveur trop ancien (avant 2.12) sous peine de voir la liste des couches s’afficher, mais une liste vide de points être retournée par les requêtes WFS. WM(T)S semble avoir toujours bien fonctionné (voir figure 2).

3. Intégration dans une page web : OpenLayers 2

La portée de notre serveur se retreint encore à un public relativement averti puisque utilisateur de QGis. Nombreuses sont les cartes sur le Web qui ne sont pas destinées à un auditoire de géographes mais simplement pour illustrer un concept à des utilisateurs probablement spécialistes dans d’autres domaines. Afin d’insérer des cartes et les couches générées sous QGis, l’environnement de développement en JavaScript OpenLayers semble idéal puisque fournissant la capacité à charger des couches dans tous les formats qui nous intéressent – WM(T)S, WFS ou GPX pour les traces GPS. Nous avons initialement été sensibilisés à OpenLayers comme environnement de développement pour accéder aux couches fournies par l’IGN au travers de Géoportail. Nous allons dans un premier temps nous intéresser à la version stable d’OpenLayers (version 2) qui semble mature, mais verrons que la nouvelle mouture (version 3) offre des fonctionnalités précieuses qui justifient la migration vers l’environnement en cours de développement.

Une carte OpenLayers est formée d’un appel aux bibliothèques JavaScript (<script>...</script>), la création d’une carte remplie d’un fond par défaut – nous choisissons le fond vectoriel OpenStreeMaps (OSM) – et recouvert de diverses couches additionnelles que nous fournirons depuis qgis-server.

Dans son expression la plus simple, une page OpenLayers ressemble à :

Et tant que nous n’insérons que des couches WMS – donc des images (puisque les données, même vectorielles, sont transmises sous forme d’images avec une résolution adaptée au facteur de grossissement), tout se passe bien. Nous pouvons ainsi ajouter une image Landsat accessible sous forme WMS, voir même un des fronts de glacier :

Le lecteur est encouragé à compléter cet exemple en ajoutant un front additionnel, par exemple la couche nommée fronts5, et observer la conséquence d’échanger l’ordre d’affichage des couches. Ces exemples sont fortement inspirés des excellents tutoriaux fournis par http://demo.3liz.com/wfst/wfs.html.

Pour WFS, le problème se corse : le principe de Same Origin impose que le serveur de données soit sur le même site que le serveur Web. Ceci est valable pour toute donnée vectorielle. Afin d’illustrer ce problème, nous proposons trois sources de données : un même fichier GPX stocké sur jmfriedt.sequanux.org, jmfriedt.free.fr, puis qgis.sequanux.org. Le serveur qgis-serveur fournissant les couches WFS est exécuté comme service Apache2 sur qgis.sequanux.org, tandis que la page Web qui fait appel à ce service peut être placée sur chacun de ces trois serveurs. La figure 3 illustre les interdictions d’accès rencontrées, tel qu’indiqué par Firebug : le fichier GPX stocké sur jmfriedt.sequanux.org peut être lu par la page http://jmfriedt.sequanux.org/jmfwfs.html mais est rejeté par la même page stockée sur http://jmfriedt.free.fr/jmfwfs.html. Cette dernière n’a accès qu’au même fichier GPX placé dans le répertoire du serveur jmfriedt.free.fr.

Les exemples de la figure 3 sont obtenus, après la définition de la carte, par l’appel aux couches vectorielles aux formats GPX ou WFS, de la forme :

Ce code fait appel au même fichier GPX stocké sur deux serveurs (afin de valider que seule la couche située sur le même serveur que la page Web appelante est acceptable), et la couche nommée fronts5 fournie au format WFS par qgis.sequanux.org/cgi-bin/project/landsat/qgis_mapserv.fcgi.

Ayant dépassé tous les écueils de configuration, nous avons finalement la satisfaction de voir toutes les couches s’afficher sur la même carte, en ajoutant l’onglet de commandes qui permet d’activer ou désactiver chaque couche afin d’aisément comparer la position des fronts des glaciers et leur évolution dans le temps :

Il n’aura cependant pas échappé au lecteur que nos images satellites semblent fortement déformées (comparer les figures 1 et 4). Il s’agit ici de la conséquence de la projection de Mercator, qui tend à considérablement grossir l’axe des abscisses lorsque la zone considérée s’éloigne de l’équateur. Quoiqu’acceptable jusqu’à des latitudes de ±60°, l’effet devient vraiment significatif au-delà et carrément désagréable par 79°N. C’est pourquoi les projections locales, tentant de trouver un plan localement tangent à la sphère représentant le globe terrestre, ne prétendent pas à s’appliquer à l’ensemble de la Terre mais uniquement à un petit segment de longitude. Malheureusement, OpenLayers 2 ne supporte pas la reprojection à la volée des données bitmap, nous sommes coincés avec la projection imposée depuis Google Maps qui continue à déformer les pôles. Nous allons remédier à ce problème en passant sous OpenLayers 3 et retrouver ainsi les belles cartes de Qgis.

4. À la pointe de la technologie : OpenLayers 3

Dans la course aux fonctionnalités, nous ne cessons de casser ce qui marche pour mettre à jour les nouvelles fonctionnalités. Pourquoi donc casser ce bel exemple OpenLayers2 pour le remplacer par OpenLayers3 ? La déformation excessive des cartes aux latitudes élevées est un défaut bien connu de la projection de Mercator utilisée par défaut par tous les environnements cartographiques de la toile (projection de code 900913 – l’écriture leet de Google – puisque la norme choisie n’a pas été validée par l’instance qu’est l’EPSG). Étant conscients des déficiences de la projection de Mercator, nous voudrions pouvoir sélectionner un mode de projection localement tangent à la région considérée – à savoir Universal Transverse Mercator (UTM). Et nous en arrivons à l’excuse pour passer de OpenLayers 2 à 3 : la version actuelle d’OpenLayers ne permet pas la projection au vol des images, et impose donc une sortie déformée, particulièrement gênante lorsque nous nous éloignons de l’équateur.

Nous reprenons donc toutes nos investigations, mais cette fois dans le contexte d’OpenLayers 3, avec la nécessité pour chaque nouvelle couche de préciser le mode de projection d’entrée (toutes nos données sont stockées en coordonnées sphériques, donc WGS84 de code EPSG:4326. OpenLayers ne connaît pas tous les modes de projection de sortie : nous devons les définir selon nos besoins (voir figure 5).

Heureusement, proj4 se charge de cette opération pour nous, et http://spatialreference.org/ref/epsg/wgs-84-utm-zone-33n/ nous informe de la façon de définir WGS84/UTM33N (la bande qui couvre la Norvège – nous aurions choisi UTM31N pour la France). Ainsi, le mode de projection qui va nous intéresser se définit par :

Et nous pouvons faire appel aux outils de conversion de la forme :

Ici la projection de la carte respecte la norme EPSG (WGS84/UTM33N s’appelle EPSG:32633) que nous avons défini par proj4 auparavant, tandis que les données en entrée (ici issues de la fonction WM(T)S du serveur) sont au format WGS84 en coordonnées sphériques (aussi nommée EPSG:4326).

5. Ajout automatique de toutes les couches WMS d’un projet qgis-server dans OpenLayers2

Une dernière question que nous pouvons nous poser, dans le contexte d’un travail collaboratif où chaque utilisateur est susceptible d’ajouter ses propres couches, tient en la capacité à automatiser la génération d’une page Web contenant toutes les couches fournies par un serveur QGis. Cela signifie donc récupérer la description XML des capacités du serveur, découper cette description pour extraire le nom des couches disponibles, et finalement générer automatique le script JavaScript affichant toutes ces couches. De telles fonctionnalités sont fournies (paquet php-xml dans la distribution Debian) par simplexml. Ainsi, un premier script qui se contente de lister toutes les couches disponibles en recherchant d’abord les fonctionnalités WMS du serveur, puis en étudiant [5]la liste des propriétés (FeatureTypeList) et finalement les propriétés de chacun de ces objets (FeatureType) pour en extraire le nom (FeatureTypeList->FeatureType->Name) ressemble à :

Ayant compris ce mécanisme, il devient simple de générer les champs décrivant chaque couche WMS affichée, en lui attribuant le même nom de couche que celui utilisé dans QGis. Le résultat est :

Ce qui donne la figure 6. Il serait d’une part très fastidieux d’ajouter à la main toutes les couches de ce projet, mais surtout dans cet exemple l’affichage s’adapte automatiquement aux nouvelles couches ajoutées par les divers contributeurs au projet.

6. Ajout d’images orthorectifiées produites par microdrone

Nous avions présenté le flux de traitement permettant de générer une série d’images orthorectifiées et de modèles numériques d’élévation associés avec le souci de comparer ces jeux de données entre eux, mais sans prétention de les positionner dans un référentiel absolu partagé par convention entre plusieurs utilisateurs. Il s’avère [6]que la solution de simplement translater un jeu de données par rapport à l’autre ne suffit plus lorsque la zone considérée s’étend sur plusieurs hectares. Nous avons observé que dans ce cas, le positionnement par GPS (mono-fréquence, telle que fournie sur microdrone DJI Phantom3 Professional) n’est exact qu’à une dizaine de mètres, alors que nous visons un positionnement au décimètre près pour des comparaisons de vols successifs. Il s’avère que vouloir corriger les défauts d’échelle et de position du modèle numérique d’élévation en lui appliquant les corrections permettant de superposer les images orthorectifiées est une approche peu judicieuse, et qu’il vaut mieux insérer dans le flot de traitement l’étape d’exploitation des points de contrôle au sol pour corriger les défauts du modèle de caméra et de leur position (http://forum-micmac.forumprod.com/campari-residuals-differ-from-the-point-cloud-t881.html). Pour notre part, les points de contrôle au sol (GCP) sont acquis a posteriori dans Google Earth – les coordonnées sphériques fournies avec 6 décimales sont précises à 11 cm à l’équateur : 7 points aisément reconnaissables dans les jeux de données sont d’une part identifiée dans Google Earth (bandes de parking, barrière, coin à la base de bâtiments) avec des coordonnées converties de coordonnées sphériques (WGS84) vers un référentiel projeté localement tangent à la Terre (WGS84/UTM31N), et d’autre part leur position (en pixel) identifiée sur les images acquises par drone. Ce couple de fichiers (coordonnées dans l’espace v.s coordonnées sur les photographies) alimente l’outil Campari de Micmac. Après traitement dans ces conditions, nous observons une exactitude de position par rapport aux couches vectorielles de Openstreetmaps mais surtout par rapport à une image aérienne de l’IGN (BD ORTHO, convertie de Lambert93 à WGS84/UTM31N puisque nous avons vu que OpenLayers 2 ne sait pas le faire à la volée) meilleure qu’une vingtaine de centimètres – erreur attribuable à notre manque d’assiduité à pointer les points de contrôle. Ce jeu de données est disponible sur http://qgis.sequanux.org/cgi-bin/project/femto/qgis_mapserv.fcgi.

Conclusion

Nous avons proposé un environnement de travail permettant de disséminer des informations géoréférencées soit vers des utilisateurs de logiciels dédiés de gestion de systèmes d’informations géographiques (SIG) ou vers les API Web associées (OpenLayers). Est-ce que ce mode de dissémination est réellement utilisé en pratique ? L’institut polaire norvégien (NPI) propose ses informations dans ce format tel que décrit sur http://geodata.npolar.no/#basemap-services. Ainsi, sous QGis, configurer l’onglet WM(T)S vers http://geodata.npolar.no/arcgis/rest/services/Basisdata/NP_Ortofoto_Svalbard_WMTS_25833/MapServer/WMTS? donne accès aux informations bitmap que sont les photographies aériennes de très haute résolution (jusquà 16 cm/pixel !) tandis que l’url http://geodata.npolar.no/arcgis/services/CryoClim/glaciers/MapServer/WmsServer? dans l’onglet WFS donne accès aux informations vectorielles que sont les limites des glaciers extraites des images satellitaires SPOT.

Le lecteur est encouragé à reproduire ces expériences sur ses propres zones d’intérêt. Compte tenu de la résolution médiocre de Landsat, le résultat est peu convaincant pour la vallée de Chamonix, mais les glaciers groenlandais (Illulisat près de la baie de Disco) ou la banquise antarctique sont parfaitement appropriés à cette démonstration, voir la mer d’Aral sur http://www.nasa.gov/mission_pages/landsat/news/40th-top10-aralsea.html (un peu pénible à assembler car de superficie supérieure à ce que le site de l’USGS permet d’analyser) ou la mer morte sur http://earthobservatory.nasa.gov/IOTD/view.php?id=77592 et pour laquelle l’extension des bassins d’exploitation du sel est particulièrement visible. Le lecteur saura-t-il calculer la surface ainsi affectée ? Indice : tracer un polygone avec l’outil New Shapefile puis dans la table des attributs, créer un nouveau champ pour chaque polygone avec la variable $area.

Une autre voie d’étude qui semble intéressante dans le cadre de la diffusion sur le Web de données géoréférencées est qgis2web, greffon de QGis actuellement inexploitable sous Debian testing/sid compte tenu d’une dépendance cassée. Finalement, étant donné que toutes les informations pertinentes à une utilisation en surface sont disponibles sur OpenStreetMaps, ne serait-il pas temps de commencer à considérer la cartographie des services sous terrains (eau, gaz, électricité), qui même s’ils sont propriété de leurs exploitants respectifs, sont devenus des services nécessaires au quotidien d’un habitant d’Europe occidentale actuel : chaque kilomètre de route contient plusieurs dizaines de kilomètres de services sous terrains [7]qui ne demandent qu’à être cartographiés à l’occasion des ouvertures de routes.

Références et Notes

[1] PIERROT-DESEILLIGNY M. etFRIEDT J.-M, « La photogrammétrie pour tous : MicMac pour la reconstruction 3D de scènes géoréférencées à partir de photographies numériques », tutorial à FOSS4G-fr 2016, décrit sur jmfriedt.free.fr/foss4g_2016

[2] FRIEDT J.-M. Friedt, TOLLE F., etBERNARD É., « Utilisation de Micmac pour la génération de modèle numérique d’élévation par traitement d’images acquises par microdrone », GNU/Linux Magazine France 191, pp.48–57 (Mars 2016)

[3] ERLE S., GIBSON R., etWALSCH J., « Building the geospatial web », dans « Mapping Hacks – Tips & Tools for Electronic Cartography », O’Reilly (2005)

[4] GIBSON R. et ERLE S., « Google Maps Hacks – Tips & Tools for Geographic Searching and Remixing », O’Reilly (2006)

[5]Les divers champs se déduisent de l’analyse de la sortie de wget -O - "http://127.0.0.1/cgi-bin/mon_projet/qgis_mapserv.fcgi?SERVICE=WFS&VERSION=1.0.0&REQUEST=GetCapabilities"

[6] LISEIN J., PINEUX N., PIERROT-DESEILLIGNY, DEGRÉ A., et LEJEUNE P., « Détection de l’érosion dans un bassin versant agricole par comparaison d’images multidates acquises par drone », Colloque Drones et moyens légers aéroportés d’observation, 26/06/2014 (Montpellier), disponible sur http://orbi.ulg.ac.be/handle/2268/171616

[7] Chaque kilomètre de route à Hong-Kong recouvre 47 km de services enterrés, http://www.uti.hk/media/attachments/2nd_ICUMAS_full_paper.pdf ou http://www.scmp.com/news/hong-kong/article/1647088/small-army-utility-specialists-keeps-hongkongers-safe-pipes-cables, ou en d’autres termes il existe 47 câbles, tuyaux, fibres optiques et autres modes de transmission de fluides et d’informations sous chaque route.

Remerciements

J.-P Culas (CM-Drones, Besançon) a effectué les vols au-dessus du parking de FEMTO-ST. Les membres du forum Micmac (http://forum-micmac.forumprod.com/) ont une fois de plus partagé leurs connaissances pour corriger les défaillances dans les traitements que nous proposions initialement.