Archives mensuelles : février 2017

Traiter des images MODIS avec R pour calculer un NDVI

La Nasa propose depuis 1999 des images MODIS (Moderate-Resolution Imaging Spectroradiometer) à différentes résolutions spatiales et temporelles. Elles sont produites par deux satellites en orbite, Terra (1999) et Aqua  (2002), qui embarquent des capteurs pour le programme Earth Observing system.

Les instruments permettent d’enregistrer 36 bandes spectrales avec une fréquence de passage de 1 à 2 jours. La résolution des images diffère selon les bandes enregistrées et varie entre 0.25 et 1 km. Nous parlions il y a quelques jours du site reverb|echo de la Nasa, il se trouve qu’un certain nombre d’images sont disponibles sur le site.

Quelles données ?

Mais pour ça il faut savoir ce qu’on cherche, et la nomenclature des produits n’est pas vraiment transparente [1].  À titre d’exemple :

  • MOD9Q1 – Surface Reflectance 8-Day L3 Global 250m
  • MOD13Q1 – Vegetation Indices 16-Day L3 Global 250m
  • MOD11A2 – Land Surface Temperature/Emissivity 8-Day L3 Global 1km

Une fois qu’on en sait un peu plus sur le type de produit que l’on veut utiliser, il va falloir le télécharger et le traiter. Dans un post précédent, j’expliquais comment configurer sa machine pour utiliser wget. Si vous optez pour cette solution, il vous faudra utiliser le MRT pour reprojeter et assembler vos images. L’interface est assez intuitive et on peut faire du traitement de masse.

Mais je suis tombé aussi ce matin sur les logs de release d’un package R que je me suis empressé de tester : MODIStsp. Si vous avez déjà regardé MRT, MODISStsp n’est qu’un moyen de rester dans R pour télécharger les images. On pourra aussi regarder du côté du package MODIS qui lui n’intègre pas d’interface graphique.

Bon mais on parlait de NDVI dans le titre !

Traiter les images NDVI : MOD13Q1

En cherchant comment traiter les images MODIS, je suis tombé sur des ressources intéressantes. Notamment une page qui vous prend par la main pour traiter les données MODIS (MOD13Q1) dans R[2]. Je vous propose donc de faire un petit point pas à pas !

J’ai téléchargé et projeté les données grâce au package MODISStsp. J’ai donc un dossier qui contient les fichiers hrf, mais aussi un autre qui contient les fichiers tiff produits par le package.

Or quand on ouvre les données (dans Qgis), les gammes de valeurs ne correspondent pas à celles attendues pour du NDVI! On devrait avoir des valeurs entre -1 et 1 !

 

Freq. pour la tuille 16 – 17/01/2016

Il va falloir retravailler le raster pour obtenir les bonnes valeurs . Heureusement « the office of outer space affaire » nous donne la marche à suivre dans R[3]. Je vous propose donc ici un script commenté et customisé pour l’occasion.

library(raster)  
library(sp)
library(doParallel)
library(rgdal)
library(magrittr)

rm(list = ls()) ##clear env.

#create stack from tif NDVI files
path.v = "~/Documents/futurSahel/MOD13Q1_TS/VI_16Days_250m_v6/NDVI/"
l.files = list.files(path.v, pattern = ".tif")


my.stack = stack(paste0(path.v,l.files)) ##create raster stack of files in a directory
# load borders
border = shapefile("~/Documents/futurSahel/Senegal_gadm.org/SEN_adm0.shp") #ToDo: insert link to the shapefile with the country borders


## We use doParallel and magrittr packages to pipe different actions (crop and mask)
registerDoParallel(6) #we will use 6 parrallel thread
result = foreach(i = 1:dim(my.stack)[3],.packages='raster',.inorder=T) %dopar% {
  my.stack[[i]] %>% 
    crop(border) %>%
    mask(border)
}
endCluster()
ndvi.stack = stack(result)

ndvi.stack = ndvi.stack*0.0001 #rescaling of MODIS data
ndvi.stack[ndvi.stack ==-0.3]=NA #Fill value(-0,3) in NA
ndvi.stack[ndvi.stack<(-0.2)]=NA # as valid range is -0.2 -1 , all values smaller than -0,2 are masked out
names(ndvi.stack) = seq.POSIXt(from = ISOdate(2016,1,17), by = "16 day", length.out = 23) # atribut name as date for each layer

my_palette = colorRampPalette(c("red", "yellow", "lightgreen")) #Create a color palette for our values
## Plot X maps in the same layout
spplot(ndvi.stack, layout=c(4, 6),
       col.regions = my_palette(16))

Ce qui nous donne pour les 23 pas de temps étudiés l’image suivante. On voit bien reverdir la partie supérieure de la carte entre juillet et octobre en 2016.

Indice NDVI à partir des données MODIS MOD13Q1 sur le Sénégal pour 2016

Pour avoir un ordre d’idée, traiter 23 images et produire le graphe, en parallélisant une partie, me prend 11 min .

Des ressources :

[1]Type d’images disponible : https://modis.gsfc.nasa.gov/data/dataprod/

[2]http://www.un-spider.org/advisory-support/recommended-practices/recommended-practice-drought-monitoring/step-by-step/R

[3]http://www.un-spider.org/sites/default/files/Skript%20VCI%20final.R

Earthdata des données SIG pour l’environnement

Voilà plusieurs mois que je n’avais pas pris le temps de faire un post sur mon blog !
On trouve toujours beaucoup d’excuses, des papiers sur le feu, un déménagement, un nouveau post-doctorat… Alors pour recommencer doucement, une petite note sur un service de la NASA que j’ai découvert aujourd’hui : Reverb | echo[1]

logo earthdata
Un peu de contexte à cette découverte, mon nouveau post-doc m’a amené à sortir de ma zone de confort, on s’éloigne de la viticulture et de l’eau (encore que), pour travailler sur la zone sahélienne du Sénégal. Je cherchais des indicateurs de désertification… Et je suis tombé sur un blog génial de Martin Brandt[2] (chercheur au Department of Geosciences and Natural Resource Management à l’université de Copenhagen (Denmark)). En plus de pointer sur des papiers scientifiques, il y a aussi quelques posts sur des outils de télédétection (GRASS-GIS) et des liens vers des pourvoyeurs de données. C’est là que j’ai découvert LE truc.

Reverb est un silo de données environnementales propulsées par la NASA qui contient des choses incroyables et un puissant moteur de recherche. Je vous engage donc à regarder cette petite vidéo d’explication[3] (en anglais).

Je ne reviendrai pas sur la manière de trouver ses données, la vidéo est là, par contre une chose géniale qui pourrait inspirer d’autres sites : une fois que j’ai fait mon petit marché, le service me propose de télécharger un fichier TXT contenant les URLs des ressources demandées. L’idée est génial. Il y a quelques semaines, j’avais du faire une centaine de copier-coller depuis le site opendata.gouv.fr pour pouvoir automatiser un traitement sur le RGP de 2012[4] (j’en parlerai peut-être dans quelques mois quand l’article sera publié).

Donc une grande satisfaction, je disais, pour cette bonne pratique… je voulais donc me lancer dans un script bash et du wget… mais les ennuis arrivent …

wget http://e4ftl01.cr.usgs.gov//MODV6_Cmp_B/MOLT/MOD09Q1.006/2016.09.05/MOD09Q1.A2016249.h16v07.006.2016258071050.hdf

Connexion à urs.earthdata.nasa.gov (urs.earthdata.nasa.gov)|198.118.243.33|:443… connecté.
requête HTTP transmise, en attente de la réponse… 401 Unauthorized
Échec d’authentification par identifiant et mot de passe.
URL transformed to HTTPS due to an HSTS policy

Mince un login ! Nooon ! Pourtant Martin Brandt ne semble pas avoir ce problème. La solution est en deux étapes. Dans un premier temps, il faudra se créer un compte « earthdata » en cliquant sur login en haut à gauche. Ensuite, et la solution est sur le wiki [5], il faudra faire un petit fichier de config « curl ».

cd ~
touch .netrc
echo "machine urs.earthdata.nasa.gov login <log>; password <password>" > .netrc
chmod 0600 .netrc

Vous l’aurez compris, il ne reste qu’à remplacer par votre login et par votre mot de passe. Et cette fois-ci les requêtes passeront !

wget http://e4ftl01.cr.usgs.gov//MODV6_Cmp_B/MOLT/MOD09Q1.006/2016.09.05/MOD09Q1.A2016249.h16v07.006.2016258071050.hdf

Sauvegarde en : « MOD09Q1.A2016249.h16v07.006.2016258071050.hdf »
MOD09Q1.A2016249.h1 100%[===================&gt;]  64,76M   344KB/s    in 1m 51s

Merci pour les données, et merci pour cette facilité d’utilisation !!!

Données MODIS issu de earthdata
Les liens

[1]https://reverb.echo.nasa.gov
[2]https://matinbrandt.wordpress.com
[3]https://www.youtube.com/watch?v=iBXzzUv3b4w
[4]https://www.data.gouv.fr/fr/datasets/registre-parcellaire-graphique-2012-contours-des-ilots-culturaux-et-leur-groupe-de-cultures-majorita/
[5]https://wiki.earthdata.nasa.gov/display/EL/How+To+Access+Data+With+cURL+And+Wget