GeoExMachina

Blog géomatique

Pytroll au coin du feu (2)

Faire feu de tout bois avec python

D'accord, Pytroll nous a bien aidés à traiter la télémétrie, mais maintenant à part en faire un poster A0 pour un mur de l'appartement, on n'a pas vraiment exploité le potentiel de ces images... Et pour la plupart des traitements, vous n'aurez besoin que de python. Oui, pas forcément besoin de logiciel très élaboré, mieux vaut comprendre ce qu'on fait que de passer par une interface graphique. Votre cerveau reste votre meilleur outil !


Penser python... ?

Et d'ailleurs au risque de vous étonner encore, ça n'est sans doute pas sur votre éditeur de texte préféré que vous passerez le plus de temps... mais dans la biblio ! Oui, restons modestes, à moins de s'attaquer à quelque chose de tout nouveau, vous aurez toujours (heureusement) une montagne d'articles scientifiques capable de vous aider ! Un NDVI pour la végétation, un NDWI pour le stress hydrique... ou un seuillage sur les canaux thermiques (températures en kelvin des canaux 3.9µm et 10.8µm, notés T3.9 et T10.8) pour la détection d'incendies (biblio).

Dans le cas de la télédétection d'incendie avec Meteosat, deux étapes successives sont nécessaires. Il faut d'abord déterminer quels sont les feux potentiels grâce à un premier indice ci-dessous :
  le jour et
  la nuit,
avec  

puis dans un second temps, confirmer ces pixels de feu potentiel par un second indice :


Ces obscures formules vont nous permettre d'examiner deux types de filtre : le filtre absolu et le filtre contextuel !
Vous vous souvenez que l'on avait pu charger les canaux grâce à Pytroll ?
#Data load
try:
	if MSG_FILE_TYPE == 'L' :
		IRchannelList = ['IR_039', 'IR_108']
		global_data.load(IRchannelList, area_extent=aire.area_extent, calibrate=1)
		print "\nDATA LOAD : OK"
	else :
		IRchannelList = ['IR_039','IR_108','IR_120']
		VISchannelList = ['VIS006','VIS008']
		global_data.load(IRchannelList, area_extent=aire.area_extent, calibrate=1)
		global_data.load(VISchannelList, area_extent=aire.area_extent, calibrate=2)
		print "\nDATA LOAD : OK"
except:
	print"\nDATA LOAD FAILED, CHECK DATA AND TIME SLOT."

La variable global_data contient désormais toutes les données des canaux, que l'on peut isoler individuellement sous forme d'array :
array039 = global_data[3.9].data

On obtient alors un tableau en deux dimensions (x et y) pour chaque bande (z). Ce qui peut être schématisé de la manière suivante :



Comme tout tableau on peut obtenir ses dimensions et son type de données :
[H,W]=array039.shape
dtype=array039.dtype

Ainsi on peut créer un nouveau tableau pour accueillir les résultats d'un filtre !
#Initialize potential fire array
potfire=sp.empty((H,W),dtype='int32')
Ici on appelle le tableau résultat potfire pour "potential fire".

Mais alors qu'est-ce qu'un filtre absolu ? Celui-ci est à utiliser lorsqu'on a un seuil fixe comme référence, et qu'on l'applique à chaque valeur de notre tableau, individuellement les uns des autres.

Le premier indice pour déterminer les feux potentiels en est un par exemple :
  le jour et
  la nuit,
avec  

Ici, on va vérifier pour chaque couple de coordonnées x, y (ou i, j dans le code) si :
  • la température dans le canal 3.9µm est supérieure à 300 kelvin le jour ou 290 kelvin la nuit
  • la température dans le canal 10.8µm est supérieure à 290 kelvin
  • la différence de température entre les deux canaux est supérieure à 15 kelvin le jour et 5 kelvin la nuit

  • On va donc utiliser deux boucles for :
    array039=global_data[3.9].data
    array108=global_data[10.8].data
    
    #Initialize potential fire array
    [H,W]=array039.shape
    dtype=array039.dtype
    potfire=sp.empty((H,W),dtype='int32')
    
    #Loop to detect potential fire pixel
    for i in range(0,H):
    	for j in range(0,W):
    		img039=array039[i,j]
    		img108=array108[i,j]
    		delta= img039-img108
    	
    		#Threshold to detect potentiel fire pixel, formula depends on time
    		#Day time 
    		if hh >= day_start and hh< day_end:
    			if (img039>300 and delta>15 and img108>290):
    				potfire[i,j]=2
    			else:
    				potfire[i,j]=1
    		#Night time
    		else :
    			if img039>300 and delta>5 :
    				potfire[i,j]=2
    			else:
    				potfire[i,j]=1
    

    Les feux potentiels sont ainsi enregistrés en valeur 2 dans le tableau de résultat, tandis que les valeurs ne répondant pas aux conditions obtiennent la valeur 1.

    Pourquoi parle-t-on de feux potentiels ?... Parce que la température de surface du globe ne suffit pas à déterminer s'il y a un incendie ou pas. En effet, si vous jetez un oeil du côté du Sahara, le soleil réchauffe tant les sols que les températures de surface dépassent facilement les 300 degrés kelvin dans la bande 3.9µm. Un incendie reste cependant un phénomène localisé (et ici on parle quand même de pixel de 3x3km... c'est dire). Si toute l'étendue du Sahara obtient la valeur 2, pas de panique ! Il faut seulement adapter la détection au contexte !
    Pour ça quoi de plus approprié qu'un filtre contextuel ? On va utiliser une fenêtre mobile pour comparer la valeur d'un pixel en particulier avec les valeurs avoisinantes.



    La fenêtre mobile, c'est un peu comme portique qu'on déplace au-dessus d'une valeur et qui recouvre les valeurs voisines. On va simplement déplacer le portique !

    #True fire array initialize
    fire=sp.empty((H,W),dtype='int32')
    
    #Window width // on détermine la largueur de la fenêtre mobile
    p=window_width
    q=(p-1)/2
    
    #Border processing // on initialise à zéro toutes les valeurs sur les bords, qu'on ne pourra pas traiter
    fire[:q,:]=0
    fire[:,:q]=0
    fire[:-q,:]=0
    fire[:,:-q]=0
    
    #True fire counter initialize
    countf=0
    
    #Loop to detect true fire
    for k in range (q,H-q) :
    	for l in range (q,W-q) :
    		#Si le pixel a été identifié comme un feu potentiel
    		if potfire[k,l]==2 :
    				
    			#Pixel value in [k,l]
    			#on recupère les valeurs de pixel sur les bandes 3.9 et 10.8
    			potf039=array039[k,l]
    			potf108=array108[k,l]
    				
    			#Pixel values in the windows for 3.9um
    			#on récupère les valeurs des pixels dans la fenêtre mobile
    			temp039=array039[k-q:(k+q+1),l-q:(l+q+1)]
    			[h,w]=temp039.shape
    			n=h*w
    			#Mean and mean absolute deviation pixel of the window for 3.9µm channel
    			#on en calcule la température moyenne et la déviation absolue moyenne dans la fenêtre mobile
    			meanpotf039=temp039.mean()
    			devpotf039=sum(sum([abs(x-meanpotf039) for x in temp039]))/n
    			
    			#Delta in [k,l]
    			deltapotf=potf039-potf108
    			
    			#Delta pixel values in the windows
    			#on calcule le delta entre les températures des deux bandes, dans l'ensemble de la fenêtre mobile
    			tempdeltapotf=array039[k-q:(k+q+1),l-q:(l+q+1)]-array108[k-q:(k+q+1),l-q:(l+q+1)]
    			#Delta mean and mean absolute deviation pixel of the window
    			meandeltapotf=tempdeltapotf.mean()
    			devdeltapotf=sum(sum([abs(x-meandeltapotf) for x in tempdeltapotf]))/n
    			
    			#si le pixel concerné répond à toutes les exigences 
    			#(température de surface suffisamment elevée par rapport au contexte)
    			if deltapotf > (meandeltapotf+level_requirement*devdeltapotf) and \
    				potf039 > (meanpotf039+level_requirement*devpotf039):
    				#source Manyangadze 
    				#"Forest fire detection for near real-time monitoring using geostationary satellites"
    				fire[k,l]=array039[k,l]
    				countf+=1
    			else:
    				fire[k,l]=0
    				
    		else:
    			fire[k,l]=0

    Les feux vérifiés (true fire) sont ainsi enregistrés avec la valeur 1 dans le tableau de résultat, tandis que les valeurs ne répondant pas aux conditions obtiennent la valeur 0. Utiliser ce filtre nous aura permis d'éviter les fausses alertes dues à un milieu globalement très chaud ! Sur l'image ci-dessous on observe les feux vérifiés en rouge et les feux potentiels en orange. Les grandes zones oranges remplissaient les conditions du filtre absolu mais, vu la superficie, il est difficile de croire que ce sont des incendies. Ou alors ils sont peut-être plusieurs... tout ça reste encore à peaufiner.


    Canal 3.9µm en nuance de gris ; pixel rouge = feu vérifié ; pixel orange = feu potentiel

    Pytroll au coin du feu (1)

    Pour du m'pop corn grillé

    Suite à la présentation de Firemsg à la conférence Eumetsat 2017, on va ici rentrer dans le détail de l'application. Cela va nous permettre d'aborder l'utilisation de pytroll, une librairie bien utile pour la télédétection à partir d'images satellites. Comme l'explique l'équipe Pytroll, leur objectif est de fournir une librairie Python pour traiter les données issues de capteurs satellites :

    " The objective of Pytroll is to provide an easy to use, modular, free and open source python framework for the processing of earth observation satellite data. The provided python packages are designed to be used both in R&D environments and in 24/7 operational production. "

    Et il faut dire que c'est réussi, une bonne liste de capteurs opérationnels, une documentation accessible et des développeurs toujours prêts à dépanner. Alors on choisit son satellite préféré et on commence !

    Personnellement, j'ai utilisé Pytroll pour le développement d'un outil de surveillance d'incendie à partir d'images Meteosat. Je vous invite à y voir le résultat par ici. Obtenir des données Meteosat est simple comme bonjour, il suffit de s'enregistrer sur le Data Center d'Eumetsat et de demander un accès aux données. L'accès par FTP permet notamment de réaliser des séries temporelles à haute répétitivité temporelle. Il suffit d'automatiser le téléchargement périodiquement (table de cron sur Linux).


    Les noms des fichiers binaires ressemblent à cela : L-000-MSG3__-MSG3________-IR_039___-000001___-201709230245-C_.
    Le nom est composé :
  • du type de fichier binaire (L pour LRIT, H pour HRIT)
  • du nom du satellite (MSG3)
  • de la bande (IR_039)
  • du segment (entre 1 et 8, il faut 8 segments pour obtenir l'image complète du globe)
  • de la date (YYYY/MM/DD/HHMM)
  • de l'état de compression (C_ = compressé, __ = décompressé)

  • Mais des fichiers binaires, ça n'est pas tout à fait simple à traiter... sauf avec Pytroll ! Il suffit juste de suivre l'installation de Pytroll et de compiler le Public Wavelet Transform Decompression Library Software. On peut alors utiliser pytroll/mpop avec python.

    Mpop fonctionne à partir de scènes que l'on définit en fonction d'une date d'acquisition et de l'aire souhaitée.
    La classe d'objet à utiliser pour cela est GeostationaryFactory (pour les satellites géostationnaires). Alors on l'importe, tout comme la fonction get_area_def pour sélectionner notre aire de découpe, et la fonction debug_on, qui permet d'obtenir les logs détaillés de mpop.

    from mpop.satellites import GeostationaryFactory
    from mpop.projector import get_area_def
    from mpop.utils import debug_on

    La fonction debug_on( ) se lance le plus simplement du monde :
    debug_on()

    La définition de la scène se fait par la fonction create_scene( ) :
    #Scene Configuration
    time_slot = datetime.datetime(YYYY, MM, DD, hh, mm)
    try:
    	global_data = GeostationaryFactory.create_scene("meteosat", "10", "seviri", time_slot)
    except:
    	print "\nSATELLITE DEFINITION LOAD FAILED, CHECK THAT meteosat10.cfg EXISTS IN THE MPOP FOLDER\
    		 OR CHANGE ARGUMENT IF YOU USE ANOTHER SATELLITE DEFINITION."
    Les deux premiers arguments permettent de retrouver le nom du fichier de configuration qui a été placé dans le dossier d'installation de mpop. Le troisième correspond au capteur utilisé (cf. fichier de config meteosat10.cfg), afin de charger les paramètres du capteur. Le quatrième argument correspond à l'heure d'acquisition de l'image (issue du nom de fichier binaire).

    L'aire de définition est définie par get_area_def( ). L'argument fait référence nom des régions définit dans le fichier areas.def.
    try:
    	aire = get_area_def("AfSubSahara")
    except:
    	print "\nAREA DEFINITION LOAD FAILED, CHECK THAT areas.def"
    On définit pour finir les bandes à charger (variable 'IRchannelList' et 'VISchannelList') et utilise la fonction load( ) pour procéder au chargement des bandes.
    #Data load
    try:
    	if MSG_FILE_TYPE == 'L' :
    		IRchannelList = ['IR_039', 'IR_108']
    		global_data.load(IRchannelList, area_extent=aire.area_extent, calibrate=1)
    		print "\nDATA LOAD : OK"
    	else :
    		IRchannelList = ['IR_039','IR_108','IR_120']
    		VISchannelList = ['VIS006','VIS008']
    		global_data.load(IRchannelList, area_extent=aire.area_extent, calibrate=1)
    		global_data.load(VISchannelList, area_extent=aire.area_extent, calibrate=2)
    		print "\nDATA LOAD : OK"
    except:
    	print"\nDATA LOAD FAILED, CHECK DATA AND TIME SLOT."
    
    L'argument calibrate permet d'obtenir les valeurs radiométriques en température de surface (calibrate = 1) ou en radiance (calibrate = 2). Ici je charge les bandes 3.9 µm et 10.8 µm pour les besoins de la télédétection d'anomalies thermiques, mais on peut charger toutes les bandes. Attention toutes les bandes ne sont pas accessibles en LRIT !

    A partir de là, les bandes sont accessibles en précisant la bande entre crochet. Celles-ci peuvent être exportées grâce à gdal :
    		driver=gdal.GetDriverByName("GTiff")
    		
    		for channel in IRchannelList:
    			#Save channel into tiff
    			file_name=file_nameBT+channel+'.tiff'
    			driver.CreateCopy(file_name,gdal_array.OpenArray(data[channel].data,None))

    A titre d'exemple voilà le type de résultat que l'on peut obtenir grâce à Pytroll !