De nombreux fonds cartographiques existent et sont largement utilisés par la communauté géomaticienne. OpenStreetMap, Google Maps, Bing Maps et bien d’autres ornent la plupart des cartographies et sont souvent présentes sur les interfaces web. Basées spatialement sur la mappemonde, ces données géographiques sont par nature des mosaïques d’images satellites, de photo-aériennes, ou encore de données vectorielles qui peuvent représentées une grande variété d’objets géographiques comme des routes, des fleuves, du bâti, des cultures, etc. Ces représentations cartographiques se rapprochent dans la majorité des cas de la réalité sur le terrain mais on observe certains fonds stylisés comme CartoDB voire insolite pour la gamme des fonds Stamen. Vous êtes plutôt adepte d’une cartographie réelle ou imaginaire?  « Le monde de la réalité a ses limites ; le monde de l’imagination est sans frontières », dixit Jean-Jacques Rousseau. Et l’une des vertus de l’informatique est de pouvoir créer son univers en toute liberté comme le fonds cartographique -non détaillé mais en toute modestie : du grand art!- produit dans cet article. A partir d’une image quelconque, l’objectif est de générer une mappemonde corrigée géométriquement dans le logiciel SIG QGIS (v 2.18.9). La donnée spatiale géoéréférencée en sortie est ensuite stockée dans GeoServer et mise en évidence sur une cartographie dynamique grâce à l’API d’OpenLayers (v4.1.1).

1. Géoréférencement de l’image dans QGIS

Une image est par définition un tableau ou une matrice de pixels caractérisés par une valeur numérique qui représente l’information d’une grandeur physique. En géomatique, les principaux pré-traitements d’une image avant leur intégration dans un SIG, le calcul d’indices, etc, sont les corrections atmosphériques, radiométriques et géométriques. Trois types de corrections géométriques existent : l’ortho-rectification, le redressement et le géoréférencement d’une image.

Le redressement a pour objectif d’éliminer toutes les distorsions d’une image (géométrie non régulière) provoquées par la prise de vue (capteur, angle de vue, etc) et l’environnement (courbure de la terre par exemple). Le géoréférencement consiste à associer chaque pixel à une coordonnée géographique afin de spatialiser l’image -à priori redressée- dans une projection cartographique donnée. Ces deux opérations suivent les mêmes étapes de correction décrits ci-après mais les processus pris en compte dans la transformation géométrique diffèrent (cf. 1.1.2). Quant à l’ortho-rectification, c’est un processus plus lourd utilisé en télédétection afin que l’image satellite corrigée apparaisse de manière identique à une prise de vue verticale ou « au Nadir ».

1.1. Principes du géoréféncement

1.1.1. Les points d’intérêt

Pour le géoréférencement, il est nécessaire de s’appuyer sur une référence spatiale qui présente des points communs avec l’image à caler. Bien évidemment, on ne considère pas l’ensemble des pixels de l’image mais on choisit quelques points stratégiques, aussi appelés points d’intérêt, points de contrôle, points d’amer ou encore points remarquables. Pour des images satellites, on peut cibler des croisements de routes, de manière plus générale des ruptures dans les changements d’occupation des sols. Les coordonnées géographiques des points d’intérêt peuvent être issues de relevés terrain (relevés GPS) ou d’une image déjà corrigée et géoréférencée. Il faut noter que plus le nombre de points d’intérêt est élevé et plus le calage de l’image sera précis dans l’espace.

1.1.2. Les algorithmes de transformation géométrique

Les autres pixels de l’image à corriger sont également associés à des coordonnées de manière automatisée par des algorithmes de transformation qui provoquent une distorsion géométrique plus ou moins importante. Cela va de la simple translation ou rotation jusqu’à la courbure du pixel et le choix de l’algorithme dépend donc de l’état et de la nature de l’image d’origine. Les algorithmes les plus utilisés sont les transformations polynomiales d’ordre variable qui génèrent des degrés différents de distorsion. Ainsi, le polynôme du 1er ordre est une transformation linéaire (ou affine) qui concerne les translations, rotations et mise à l’échelle de l’image; ce qui peut suffire pour un géoréférencement. Le polynôme du 2e ordre, plus utilisé pour le redressement d’images, modifie la géométrie avec la prise en compte des courbes. Enfin, les modifications géométriques apportées par le polynôme d’ordre 3 sont plus importantes et concernent des données nécessitant une rectification importante. On peut par exemple évoquer les images avec de forts dénivelés pour lesquelles les transformations géométriques pourront s’appuyer sur des Modèles Numériques de Terrain (MNT) et un nombre élevé de points d’intérêt. D’autres algorithmes existent comme la transformation linéaire, la transformation d’Helmert ou encore celle de Thin Plate Spline (TPS).

1.1.3. Erreur liée à la correction géométrique

La correction géométrique de l’image est caractérisée par des écarts moyens entre les coordonnées géographiques réelles du point d’intérêt saisies et celles calculées qui définissent le Root Mean Square (RMS) ou Erreur Quadratique Moyenne (EQM). Pour simplifier, le RMS est la racine carrée de la moyenne des carrés des écarts entre les points calculés et les points de calage. La distance entre chaque couple de points (points d’intérêt souhaités et points calculés par l’algorithme) définit une erreur résiduelle qui s’exprime en unité de référence (par exemple en mètre, en pixel, etc). Ainsi, une valeur du RMS faible signifie que la variance des écarts résiduels est faible et donc que la correction géométrique est correcte. Pour diminuer ce coefficient, chaque erreur résiduelle de la série des points doit être analysée et les points de calage avec une distance résiduelle élevées sont modifiées.

1.1.4. Le ré-échantillonage des pixels corrigés

Lors de cette transformation géométrique, l’image produite ne possède pas la même matrice que l’image d’origine et il est donc nécessaire de déterminer les valeurs numériques des pixels de l’image corrigée : on parle de ré-échantillonage. Différentes méthodes de ré-échantillonage existent :

• La méthode du plus proche voisin : Les pixels de l’image corrigée prennent la valeur du pixel le plus proche dans l’image initiale. Les valeurs numériques de la nouvelle matrice sont donc peu altérées par rapport à celle d’origine, ce qui est préférable pour les valeurs de luminance par exemple.
• La méthode par interpolation bilinéaire : on affecte au nouveau pixel la moyenne pondérée par la distance des quatre pixels environnants (2×2) de l’image d’origine.
• La méthode par convolution cubique ou interpolation bicubique : elle calcule la moyenne pondérée des 16 pixels (4×4) de l’image originale les plus proches du nouveau pixel. Cette méthode entraîne un lissage des valeurs de la nouvelle image.

Pour réaliser le géoréférencement des données images (ou raster), QGIS met à disposition un module assez simple d’utilisation qui est documenté sur ce lien : https://docs.qgis.org/2.8/fr/docs/user_manual/plugins/plugins_georeferencer.html.

1.2. Module de géoréférencement dans QGIS

Ici, l’image à caler spatialement est de format JPEG et représente une mappemonde -je vous avais prévenu : du grand art!- (figure 1).

Figure 1 : Image de format JPEG à géoéréférencer

N’ayant pas de points d’amer de référence sur toute la planète, on s’appuie sur une donnée spatiale déjà calée. Dans QGIS, on installe le plugin OpenLayers qui permet d’afficher Bing Maps, Google Maps, etc. Ici, on utilise Bing Maps Aerial affiché dans la projection cartographique WGS 84 /Pseudo Mercator (EPSG 3857).

Dans l’onglet Raster > Géoréférencer > Géoréférencer, le module de géoréférencement s’ouvre et on importe l’image à caler. L’image importée n’a pas de projection cartographique associée, il est important de choisir la même projection que la donnée d’appui. On ajoute alors des points de contrôle sélectionnés sur l’image d’origine puis les coordonnées géographiques des points sont fournies depuis le canvas, c’est-à-dire depuis les sélections sur l’image de référence Bing Maps. La figure 2 montre la table des points de contrôles que l’on peut enregistrer, et les distances résiduelles en pixel entre les points sources de l’image à caler et les points de destinations sélectionnés sur Bing Maps. Dans l’exemple, les points de contrôle ont été pris rapidement, ce qui entraîne des résidus importants. Le point 7 à l’est de la Russie est éliminé car il a un résidu de plus de 1900 pixels.

Figure 2 : Table des points de contrôle avec les distances résiduelles en pixels.

On configure la transformation géométrique (figure 3) : ici, on choisit l’algorithme polynomiale d’ordre 2 avec la méthode de ré-échantillonage du plus proche voisin. L’image est transformée dans la projection cartographique WGS 84 / Pseudo Mercator et la donnée raster en sortie est enregistrée en format GeoTiff.

Figure 3 : Configuration de la transformation géométrique de l’image à caler.

L’image géoréférencée en sortie peut être superposée à l’image d’appui afin d’observer les erreurs de calage. Sur la figure 4, on note que les erreurs importantes se trouvent dans les zones géographiques où aucun point de contrôle n’a été saisi (par exemple en Eurasie). L’idée est donc d’ajouter d’autres points d’appui de manière homogène sur toute la cartographie afin de diminuer les erreurs résiduelles.

Figure 4 : Superposition de l’image géoréférencée à l’image d’appui.

QGIS nous a permis de géoréférencer une image représentant la mappemonde, il est désormais temps de la publier sur internet.

2. Affichage du raster tuilé en webmapping

Dans l’état, le raster géoréférencé ne peut pas être publié directement sur une interface web via l’API d’OpenLayers car les formats images GeoTiff et Tif ne sont pas pris en charge directement. De plus, les données raster sont en général des fichiers lourds, il n’est donc pas judicieux de les publier telles quelles. Il faut alors décomposer l’image en dalles (ou tuiles) pré-calculées selon la projection et l’étendue de la carte : on parle de tuilage des données spatiales. Pour réaliser ce procédé, l’utilisation d’un serveur cartographique est indispensable.

2.1. Intégration de la donnée dans GeoServer

GeoServer prend en charge plusieurs formats de données raster dont le GeoTiff. On créé un espace de travail dédié à la donnée puis un entrepôt pour lequel on choisit le répertoire de la donnée (figure 5).

Figure 5 : Création de l’entrepôt Mappemonde

Ensuite, dans l’onglet couches, on ajoute la donnée en l’éditant (nom, titre, description) et en calculant l’emprise géographique selon le système de projection géographique choisi (figure 6).

Figure 6 : Edition de la couche SIG avec le calcul de l’emprise géographique.

Basiquement, le raster géoréférencé est prêt à être prévisualisé par exemple dans OpenLayers. Dans l’onglet « Prévisualisation des couches », on affiche la donnée appelée « georef » de l’espace de travail « test » par le protocole WMS via une URL de type :

http://localhost/geoserver/test/wms?service=WMS&version=1.1.0&request=GetMap&layers=test:georef&styles=&bbox=-1.9587787078002207E7,-8934317.725420557,2.1594936455664076E7,1.940261280468178E7&width=768&height=528&srs=EPSG:3857&format=application/openlayers.

Toutefois, comme évoqué précédemment, on veut afficher le fonds cartographique en tuiles. GeoServer implémente le service Web Map Tile Service  (WMTS v 1.0.0) grâce au logiciel GeoWebCache qui permet de cacher (mettre en mémoire) et de tuiler les couches SIG. Dans l’onglet « Cache des tuiles », cinq grilles de tuilage sont mises à disposition à l’installation du serveur cartographique dans les projections 4326 et 900913 (figure 7) et elles ne peuvent pas être modifiées. L’unité du pixel pour l’EPSG 4326 est le degré et pour l’EPSG 900913 le mètre.

Figure 7 : Grille de tuilage EPSG:900913 interne à GeoServer.

Chaque couche SIG ajoutée dans le serveur cartographique est par défaut cachée et tuilée selon les grilles de tuilage disponibles. En résumé, la couche est découpée en dalles d’une certaine dimension et résolution spatiale calculée selon l’étendue de la carte, l’échelle et la projection choisie.

2.2. Affichage du fonds cartographique en WMTS avec OpenLayers

Dans ce paragraphe, on utilise OpenLayers pour afficher l’image raster via le protocole WMTS. Deux classes sont importantes à appréhender : ol.source.WMTS qui permet de charger la source de la donnée tuilée depuis le serveur cartographique et ol.tilegrid.WMTS qui définit la grille de tuilage pour cette même source.

Dans la déclaration de la source, les paramètres requis sont donc la grille de tuilage tileGrid, la couche layer, le style style (non nécessaire avec GeoServer) et le nom du jeu de la matrice du tuiles matrixSet. Sur la figure 7 ci-dessus, on note par exemple que ce dernier est nommé ESPG:900913.

Dans tileGrid, doivent être paramétrés au minima l’étendue extent ou l’origine origin de la grille de tuilage, un tableau des tailles du pixel resolutions en adéquation avec le tableau des IDs (ou niveau) du jeu de la matrice matrixIds  (on peut se référer à la figure 7 pour observer ces éléments). Il est à noter que si l’origine ou les origines de la grille de tuilage ne sont pas paramétrées, l’origine doit être définie à partir du coin en haut à gauche de l’étendue de la carte.

Pour le script, on s’appuie sur l’exemple intitulé WMTS de l’API d’OpenLayers mis à disposition sur le site et le jeu de grille de tuilage interne de GeoServer EPSG:900913 qui est défini pour le système de projection EPSG:900913 avec une largeur et une hauteur de tuile de 256 pixels et le mètre comme unité (1 mètre correspondant à une unité).

Après les chargements des fichiers CSS et JavaScript d’OpenLayers et la déclaration de la carte dans le body (div map), on écrit le code JavaScript et on commence par définir les options origin, resolutions et matrixIds de la grille de tuilage ol.tilegrid.WMTS.

On définit l’origine par le coin en haut à gauche de l’étendue de la carte pour la projection ESPG:900913.

// Déclaration de la projection
var projection = ol.proj.get('EPSG:900913');
// Déclaration de l'étendue de la carte
var projectionExtent = projection.getExtent();
// Déclaration de l'origine en haut à gauche de l'étendue de la carte
var originTopLeft = ol.extent.getTopLeft(projectionExtent);

La taille du pixel (en mètre) est calculée en fonction du nombre de tuiles de 256 pixels lié à la largeur de l’étendue de la carte pour chaque niveau. Par exemple, au niveau de zoom 0, on a une largeur de la mappemonde égale à 40 075 016 m, soit pour 1 dalle de 256 pixels, une taille de pixel de 156 543 m environ.

// Déclaration du tableau de la taille des pixels
var resolutions = new Array(31);
// Déclaration du tableau des ids (Noms) du jeu de la matrice
var matrixIds = new Array(31);
// Unité de la tuile = Largeur de l'étendue de la carte / Taille de la tuile
var size = ol.extent.getWidth(projectionExtent) / 256;
// Pour chaque niveau de zoom 
for (var z = 0; z < 31; ++z) {
	// Intégration de chaque nom du jeu de matrice
	matrixIds[z] = "EPSG:900913:" + z;
	//Taille du pixel =  Unité de la tuile / Nombre de tuiles sur une largeur pour chaque zoom
	resolutions[z] = size / Math.pow(2, z);
	//On retrouve les données du jeu de la matrice pour chaque niveau de zoom
	console.log('Niveau : '+z+ '| Taille du pixel : '+resolutions[z]+' | Nom : '+matrixIds[z]+' | Tuiles : '+Math.pow(2, z)+'x'+Math.pow(2, z) );
}

On a donc défini les paramètres de la grille de tuilage (figure 8), ici elle correspond à la grille interne EPSG:900913 de GeoServer (figure 7).

Figure 8 : Déclaration des paramètres de la grille de tuilage dans OpenLayers.

On peut donc déclarer la grille de tuilage  WMTS :

// Déclaration de la grille de tuilage
var tileGrid = new ol.tilegrid.WMTS({
	origin: originTopLeft,
	resolutions: resolutions,
	matrixIds: matrixIds
});

La source et la couche WMTS par la classe ol.layer.Tile sont ensuite initialisées :

// Source WMTS
var sourceWMTS = new ol.source.WMTS({
	attributions: '© <a href="http://www.geomatick.com">Geomatick</a>',
	// Service wmts par GeoWebCache dans GeoServer
	url: 'http://localhost/geoserver/gwc/service/wmts',
	// Notre magnifique mappemonde géoréféencée
	layer: 'test:georef',
	// Le nom du jeu de matrices de tuiles
	matrixSet: 'EPSG:900913',
	format: 'image/png',
	projection: projection,
	// grille de tuilage
	tileGrid: tileGrid,
	// GeoServer gère le style configuré par défaut dans le serveur
	//style: '...',
	// Couverture du monde à l'horizontal
	wrapX: true
})
// Couche WMTS
var coucheWMTS = new ol.layer.Tile({
	opacity: 1.0,
	source: sourceWMTS
});

Enfin, la map est initialisée :

var map = new ol.Map({
		controls: ol.control.defaults({
		attribution : true,
	}),
	target: 'map',
	layers: [coucheWMTS],
	view: new ol.View({
		center: [-11158582, 4813697], //900913
		zoom: 3
	})
});

On peut désormais ajouter d’autres couches SIG à notre magnifique fond cartographique (figure 9)!

Figure 9 : Représentation du raster géoréférencé en webmapping.

En conclusion, cet article a permis d’expliquer comment construire un fond cartographique tuilée -ici minimaliste- en webmapping à partir d’une seule image quelconque. Le géoréférencement du raster a été réalisé grâce au logiciel QGIS et son plugin géoréferencer. L’utilisation de Bings Maps a permis de sélectionner des points d’appui pour générer la transformation géométrique de l’image. Ensuite, cette dernière a été stockée dans GeoServer afin de générer le tuilage de la donnée SIG. Le service WMTS est fourni grâce à l’implémentation de GeoWebCache dans le serveur cartographique. Plusieurs grilles de tuilage internes sont déjà paramétrées, on a utilisé la grille EPSG:900913, compatible avec Google Maps et Bings pour générer notre fonds cartographique sur OpenLayers.

Concernant Leaflet, à ce jour, il est possible de charger des données tuilées et un plugin WMTS pour afficher des fonds IGN a été développé mais il ne semble pas qu’une classe WMTS soit fournie par l’API. En perspective, on consacrera un article spécifique pour le tuilage et la mise en mémoire des données spatiales afin notamment d’afficher des fichiers lourds.