Tuto : créer une image satellite à partir de données Landsat avec Python
Que ce soit pour étudier la végétation, évaluer la température de la surface terrestre, cartographier les effets d’une inondation ou tout simplement pour espionner la Corée du Nord, on a toujours une bonne raison d’utiliser des données satellitaires. Et ça tombe bien : plusieurs programmes d’observation de la Terre mettent leurs données en accès libre.
Dans ce tutoriel, nous allons voir :
- Comment accéder aux données du programme américain Landsat
- Comment utiliser ces données pour produire une image composite “vraies couleurs” de la surface avec Python
A propos : Callendar est une start-up spécialisée dans le développement de solutions innovantes pour l’évaluation des risques climatiques. Conscients du défi que représente l’adaptation au changement climatique, nous nous efforçons de partager notre expertise au travers d’outils gratuits ou de tutoriels comme celui-ci.
Comment télécharger des données Landsat ?
L’accès aux données du programme Landsat peut se faire via le site Earth Explorer. Pour cet exemple, nous allons chercher des enregistrements réalisés au-dessus de la région parisienne pendant l’été 2019.
Le satellite Landsat 8 réalise 700 prises de vue (ou “scènes”) par jour dans plusieurs bandes de fréquences de l’ultraviolet à l’infrarouge en passant évidemment par le visible. Il couvre l’ensemble de la planète avec une résolution de 30 mètres tous les 16 jours… Et toutes ces données sont accessibles gratuitement !
Voici comment procéder :
- Pour commencer, il va vous falloir créer un compte. Vous pouvez explorer les données sans mais vous en aurez besoin pour les télécharger à la fin.
- Le premier onglet, “Search criteria”, vous permet de sélectionner la zone d’intérêt. De nombreuses méthodes s’offrent à vous, la plus simple consiste à cliquer directement sur la carte à l’endroit qui vous intéresse, dans notre cas : Paris.
- Un peu plus bas dans le même onglet, vous pouvez sélectionner la période de temps qui vous intéresse. N’hésitez pas à prendre une période large : le satellite passe deux fois par mois environ mais la couverture nuageuse peut rendre certains passages inexploitables. Dans notre cas, la sélection va du 1er juin au 31 août 2019 :
- Vous pouvez maintenant passer à l’onglet suivant, “Data Sets”, pour choisir les données que vous allez télécharger. Vous allez voir que le choix est vaste… Mais ce qui nous intéresse, c’est seulement les données Landsat 8 OLI/TIRS :
- Vous pouvez maintenant aller à l’onglet suivant, “Additional Criteria”, notamment pour écarter automatiquement les images trop nageuses, par exemple avec une couverture superieure à 20%.
- Il est maintenant temps de passer à l’onglet “Results”, vous y trouverez la liste des scènes qui correspondent à votre recherche. Vous pouvez voir la zone couverte par chaque scène en cliquant sur l’icone en forme de pied :
- Si la scène vous convient, cliquez sur l’icone caddie pour l’ajouter à votre commande. Une fois la sélection terminée, cliquez sur “Item basket” à droite au-dessus de la carte (vous avez besoin de posséder un compte et d’être identifié). Vous accédez au récapitulatif de votre commande :
- Cliquez sur “Proceed to checkout”. Vous recevrez un mail de confirmation et un nouveau mail lorsque les données seront prêtes à être téléchargées. Pas d’impatience : l’attente peut durer plusieurs heures…
- Lorsque vous recevez un mail avec l’objet “USGS ESPA product request order number ************** available for download”, cliquez sur le lien. Sur la page, cliquez sur le lien qui se trouve sous “Order ID”, vous pouvez finalement télécharger vos scènes une par une en cliquant sur “Download” :
Comprendre les données Landsat
Bravo, vous avez téléchargé vos premières données satellitaires. Mais qu’avez-vous récupéré exactement ? Regardons ça de plus près.
Chaque scène se présente sous la forme d’un fichier compressé d’environ 400Mo. Vous pouvez décompresser ce fichier par exemple avec Winrar ou 7zip.
Une fois décompressé, vous trouvez 10 fichiers .tif et un fichier .xml. Le fichier xml vous donne toutes les informations que vous pouvez souhaiter sur la scène : coordonnées, dates, etc.
Les fichiers TIFF sont des images tramées (ou raster). Chacun de ces fichiers contient des données différentes. Inutile de rentrer dans les détails pour l’instant (s’il vous intéresse vous pouvez les trouver ici), il vous suffit de savoir que :
- la bande 2 (le fichier dont le nom se termine par “band2.tif”) enregistre le rayonnement émis par la surface terrestre à une longueur d’onde comprise entre 450 et 510 nanomètres, ce que notre œil verrait comme du bleu,
- la bande 3 enregistre le rayonnement entre 530 et 590 nanomètres, ce qui correspond à la couleur verte,
- la bande 4 se situe entre 640 et 670 nm, c’est-à-dire dans le rouge.
L’addition de ces trois bandes correspond au spectre visible par l’œil humain. En utilisant ces trois fichiers, nous allons donc pouvoir construire une image qui nous semble naturelle, un peu comme ce que nous aurions vu en étant à bord du satellite.
Mais avant ça, essayons de comprendre un peu mieux ce que sont ces fichiers TIFF. En réalité, ce ne sont pas des objets très exotiques, vous pouvez les voir de deux façons :
- Un fichier .tif est une image numérique sans aucune compression. Vous pouvez d’ailleurs ouvrir ces fichiers avec des logiciels de traitement d’image grand public comme Xnview ou paint. Dans ce cas, vous verrez quelque chose qui ressemble à une image en noir et blanc, les zones claires représentent les surfaces qui émettent beaucoup dans la bande de fréquence que vous avez ouvert, les zones sombres celles qui émettent moins.
- Autre façon de le voir : un fichier .tif est une matrice ou un tableau contenant des chiffres. Chaque chiffre correspond au rayonnement émis dans le bande de fréquence du fichier par une surface de terrain de 30 mètres par 30 mètres, plus ce chiffre est élevé plus le rayonnement est important, 0 correspondant à un rayonnement nul.
Créer une image en couleurs naturelles avec Python
Comme nous l’avons vu, les bandes 2, 3 et 4 correspondent à la lumière telle que notre œil la perçoit avec, respectivement, le bleu (longueur d’onde entre 450 et 510 nanomètres), le vert (530–590 nm) et le rouge (640–670 nm). Ces données peuvent être utilisées pour composer une image proche de ce qu’aurait vu un œil humain s’il avait été à la place des capteurs de Landsat 8. C’est ce que nous allons faire maintenant avec Python.
En plus des grands classiques, nous allons utiliser rasterio et surtout earthpy. Earthpy est une librairie développée par l’université du Colorado qui va grandement simplifier la manipulation de nos données.
import os
from glob import glob
import matplotlib.pyplot as plt
import rasterio as rio
import earthpy.spatial as es
import earthpy.plot as epLa première étape consiste à récupérer la liste des images pour les différentes bandes de fréquence et à la classer, de façon à ce que les bandes visibles soient en 2e, 3e et 4e position dans la liste.
#Dossier dans lequel se trouvent les données
path = os.path.join("data")#Création de la liste des chemins vers les fichiers .tif à utiliser
bandeliste_landsat_data_band = glob(path + "/*band*.tif")#Classement
liste_landsat_data_band.sort()
Nous allons maintenant “empiler” les données pour chaque bande de fréquence dans un fichier unique. C’est simple comme bonjour grâce à earthpy :
#Nom et emplacement du fichier à créer
stack_path = os.path.join("data", "stack.tif")#Création et enregistrement du raster
es.stack(liste_landsat_data_band, stack_path)
Remarquez que nous aurions pu ne prendre que les 3 bandes de fréquences qui nous intéressent. Mais nous aurons peut-être besoin des autres bandes plus tard : maintenant notre fichier global est enregistré, nous pourrons l’utiliser sans repasser par cette étape.
On peut même supprimer les données originales, histoire de laisser respirer le disque dur:
#Suppression des données originales
for f in liste_landsat_data_band:
os.remove(f)L’étape suivante consiste à ouvrir le fichier que l’on vient de créer avec rasterio :
#Nom et emplacement du fichier (a priori inchangé)
stack_path = os.path.join("data", "stack.tif")#Ouverture
with rio.open(raster_path) as src:
raster = src.read()
Dernière étape, créer une image à partir des bandes visibles. Là encore, earthpy rend l’opération incroyablement simple :
ep.plot_rgb(raster,
rgb = (3, 2, 1))plt.show()
L’argument rgb correspond aux bandes à utiliser dans l’ordre red, green, blue/rouge, vert, bleu. Il s’agit des 4e, 3e et 2e bandes mais comme l’indexation commence à zéro : rgb = (3, 2, 1).
Quelques arguments complémentaires peuvent être utiles :
- title vous permettra de rajouter un titre
- figsize pour fixer la taille de l’image
- stretch = True pour corriger la luminosité
Par exemple :
ep.plot_rgb(raster,
rgb = (3, 2, 1),
stretch = True,
figsize = (79, 80),
title = "Ile de France\n(Image composite, Landsat)")plt.show()
Et voilà le travail :
Évidement sur une scène de 185 kilomètres par 180, on ne voit pas grand chose au premier coup d’œil. Pour cela, vous pouvez enregistrer cette image et zoomer sur les points qui vous intéressent, par exemple :
Vous voyez certainement déjà des applications que ces images peuvent avoir dans l’aménagement du territoire, le suivi des grandes infrastructures, l’agriculture ou la foresterie et certainement beaucoup d’autres domaines… Et nous n’avons utilisé que 3 des 10 fichiers de données disponibles : beaucoup d’informations sont encore inexploitées !
