GeoExMachina

Blog géomatique

PisteCreator sur Qgis

Mieux que le p'tit poucet en forêt

Un billet pour prévenir de la sortie du plugin PisteCreator pour Qgis (version 2.14 minimum) ! PisteCreator a été accepté par la communauté Qgis et figure à présent dans la bibliothèque des extensions Qgis. Il est donc disponible directement dans le panneau d'extension de votre Qgis ! Attention le plugin est encore expérimental, pensez à activer l'affichage des extensions expérimentales dans les paramètres.


Désormais, plus besoin d'aller chercher le plugin sur Github. Gardez-y toutefois un oeil, ça reste encore le meilleur moyen de suivre les futurs développements !
Mais désormais il suffit juste de cocher PisteCreator dans le panneau d'extension !


Pour le tuto d'utilisation, j'écrirai un autre billet d'ici peu !

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