Images satellite et Leaflet

J'étais en train d'imaginer les prochaines randonnées après le déconfinement. La vue de ma terrasse sur le Massif du Mercantour me montrait que la limite d'enneigement, en ce début mai, était vers 2.300 mètres. Mais est-ce que la Baisse de Prals était accessible ? Est-ce que le Col de Salèse était encore sous la neige ?

Les satellites Sentinel-2 de l'ESA passant tous les 5 jours au-dessus de nous, je me suis intéressé à leurs images, et à comment les afficher sur une carte avec l'application Javascript Leaflet. Récit et démonstrations.

Images Sentinel-2

Les images prises par les satellites de l'Agence Spatiale Européenne (ESA) sont disponibles gratuitement, mais encore faut-il les trouver et les récupérer.

Nous avons opté pour l'outil PEPS du CNES. Créer un compte est rapide, tout comme son explorateur qui permet de sélectionner les dates et les types d'images souhaités. En dessinant une zone rectangulaire sur le Mercantour, en spécifiant la collection "SENTINEL-2 tuilés" et un intervalle de dates récentes, apparaissent des polygones roses qui représentent les fichiers disponibles.

En cliquant sur un polygone, une icône d'information permet de vérifier visuellement qu'il n'y avait pas trop de nuages lors de cette prise vue, puis une autre permet de télécharger le fichier. Un moment plus tard (le fichier peut faire 700 Mo), l'ouvrir avec votre outil ZIP préféré, naviguer jusqu'au dossier GRANULE puis IMG_DATA, et extraire le fichier dont le nom se termine par TCI (True Colour Image).

Le fichier au format JPEG 2000 a une taille de 10980 x 10980 pixels couvrant un carré de 100km, donc 1 pixel représente 10 mètres, on parle de résolution de 10 mètres. Il peut être visualisé directement dans un éditeur d'images comme par exemple Affinity Designer.

Dernière chose à récupérer, les coordonnées géographiques de l'image, c’est-à-dire la position de leurs coins supérieurs gauche et inférieurs droit dans l'espace, exprimées en longitude et latitude.

Si l'on ne dispose pas (pas encore, cela va venir plus bas) d'un outil capable d'extraire les métadonnées d'un fichier JPEG 2000, le plus rapide est de les récupérer dans le INPIRE.xml trouvé aussi dans le fichier ZIP :

<gmd:EX_GeographicBoundingBox>
  <gmd:westBoundLongitude>
    <gco:Decimal>6.495928592903788</gco:Decimal>
  </gmd:westBoundLongitude>
  <gmd:eastBoundLongitude>
    <gco:Decimal>7.883505557566801</gco:Decimal>
  </gmd:eastBoundLongitude>
  <gmd:southBoundLatitude>
    <gco:Decimal>43.23826255332707</gco:Decimal>
  </gmd:southBoundLatitude>
  <gmd:northBoundLatitude>
  <gco:Decimal>44.24783068396879</gco:Decimal>
  </gmd:northBoundLatitude>
</gmd:EX_GeographicBoundingBox>

Affichage sous forme d'image dans Leaflet

Il est assez simple de superposer de manière sommaire le fichier extrait à une carte dans Leaflet, comme montré dans l'exemple ci-dessous (cliquer sur le bouton HTML pour voir le code source). Comme les navigateurs ne reconnaissent pas le format JPEG 2000, il faut convertir l'image vers le format JPEG (le format PNG génère des fichiers beaucoup trop gros). Mettre finalement les coordonnées géographiques de l'image dans la variable bounds, et le résultat s'affiche :

Mais il y a un important problème de précision, que l'on perçoit bien en zoomant par exemple sur le bord de mer de Vintimille, et notamment sur son port circulaire : l'image satellite est décalée de quelques centaines de mètres par rapport à la carte.

Le problème vient d'une incohérence entre le système de coordonnées utilisé par l'ESA et celui utilisé par les cartes en ligne. Plus précisément, l'ESA utilise, comme les systèmes GPS, le système WGS84, aussi nommé EPSG:4326. Alors que les cartes comme OpenStreetMap ou Google Maps utilisent le système EPSG:3857. Une bibliothèque comme Leaflet sait utiliser les deux systèmes, mais pas en même temps. Il va donc être nécessaire de reprojeter (ou déformer) l'image satellite vers le système de coordonnées EPSG:3857.

Reprojection de l'image pour Leaflet

L'outil de référence pour le traitement d'images géographiques est GDAL. Si son installation peut être ardue, selon utilisation pour notre reprojection est assez simple. En repartant du fichier JPEG 2000 d'origine, qui contient les informations sur le système et les coordonnées géographiques, voici la ligne de commande pour l'outil GDALWARP pour la reprojeter en EPSG:3857 :

gdalwarp -t_srs EPSG:3857 T32TLP_20200429T102559_TCI.jp2 sentinel-demo2.tiff

Vous noterez que l'on est obligé de créer un fichier TIFF, car GDALWARP sait lire mais pas écrire des fichiers JPEG 2000.

Pour extraire les coordonnées géographiques du fichier résultant, l'on dispose maintenant de l'outil GDALINFO. Il est nécessaire de demander le résultat au format JSON pour avoir des degrés décimaux, au lieu du format degrés, minutes et secondes.

gdalinfo -json sentinel-demo2.tiff

Dans le résultat de cette commande, les coordonnées des coins sont en première et troisième position, donc 44.2478307,6.4959286 et 43.2382709,7.888669 :

"wgs84Extent":{
  "type":"Polygon",
  "coordinates":[
    [
      [
        6.4959286,
        44.2478307
      ],
      [
        6.4959286,
        43.2382709
      ],
      [
        7.888669,
        43.2382709
      ],

Il n'y a plus qu'à convertir le fichier TIFF en JPEG, modifier la variable bounds, et voici le résultat, la carte et l'image sont parfaitement superposées :

Seul inconvénient, le fichier JPEG faisant plus de 20 Mo, les déplacements dans l'image sont saccadés, même avec un ordinateur rapide. La solution, celle retenue dans toutes les applications de cartographie, est le tuilage de notre image satellite.

Tuilage de l'image

Ce processus consiste à découper l'image de départ de 10980 x 10980 pixels en une multitude de petites images de 256 x 256 pixels, et cela pour chaque niveau de zoom auquel nous souhaitons exploiter l'image. Nous avons déterminé de manière empirique que les niveaux de zoom souhaitables allaient de 8 à 16, et avons ainsi paramétré l'outil GDAL2TILES :

gdal2tiles.py --webviewer=none --no-kml -e --zoom=8-16 --processes=4 sentinel-demo2.tiff dossier_tuiles

Le traitement est long, environ une heure sur une machine Linux virtuelle hebergée d'entrée de gamme, et l'arborescence résultante comprend près de 90.000 fichiers pour un volume de 1,5 Gigaoctets. Nous attendons avec impatience la future version 3.1.0 de GDAL qui nous permettra de générer des tuiles de 512 x 512 pixels, et donc de réduire drastiquement le nombre de fichiers.

Voici le résultat, dont les tuiles sont servies par la machine Linux précitée. En bonus, un petit contrôle en bas à gauche qui permet de modifier l'opacité de l'image satellite.