DATA\WAX

Author: Louise Leroux

Corrélation de deux rasters multibandes avec R

Toujours dans le cadre de mes travaux de recherches, je devais comparer deux façons de calculer un NDVI cumulé annuel, 1) en intégrant toutes mes valeurs bi-mensuelles sur toute l'année, 2) en intégrant mes valeurs correspondant uniquement à la saison des pluies, c'est à dire environ de Mai à Octobre (contexte africain, au Niger).

Mes données d'entrée sont donc deux rasters calculés:

  • raster1 = la manière 1 et comptant autant de bandes que d'années (ici 11)
  • raster2 = la manière 2 et comptant autant de bandes que d'années (ici 11)

​Le résultat que je souhaitais obtenir en sortie était quelque chose comme ceci :

scatter

mais avec un graphique par année (ou par bande).

J'ai donc décidé de faire cela sur …. R (^^).

Voici donc le code que j'ai mis en place. Dans mon cas j'avais besoin de le faire tourner sur plusieurs zones donc sur plusieurs images, il est donc fait pour tourner sur plusieurs paires de fichiers, mais en supprimer la première boucle vous pouver le faire tourner sur uniquement deux fichiers.

######### Comparaison de rasterS multibandes #######

# input : 2 Rasters multibandes format ENVI (dans ce cas si, mais fonctionne avec autres formats raster)
# output : graphique au format pdf de la "corrélation" de la bande 1…n d'un raster1 avec la bande 1…n d'un raster2
# Leroux Louise, CIRAD, UMR TETIS 2014

######## Initialisation #########
## Chargement de packages nécéssaires

library(raster)
library(sp)
library(Hmisc)

## Initialisation de l'espace de travail
setwd("D:/MyFolder")

##Importation des données raster##
# Pour le premier raster

liste.file1<-list.files(path=".", pattern="motif1")
liste.file1.hdr<-liste.file1[-grep(« .hdr », liste.file1, fixed=T)] #Récupération du nom sans l'extension
 

# Pour le second raster
liste.file2<-list.files(path=".", pattern="motif2")
liste.file2.hdr<-liste.file2[-grep(« .hdr », liste.file2, fixed=T)]  #Récupération du nom sans l'extension

## Boucle sur les couples de fichiers ##
for (j in 1:length(liste.file1.hdr)){
  name1<-liste.file1.hdr[j]
  name2<-liste.file2.hdr[j]
  nband<-nbands(raster(name1)) #variable stockant le nombre de bandes dans l'image
  cell<-ncell(raster(name1)) #variable stockant le nombre de pixel d'une bande
  mat1<-matrix(0, ncol=nband, nrow=cell)
  mat2<-matrix(0, ncol=nband, nrow=cell)
  
  # Stockage des bandes du raster 1 dans un matrice
  for (i in 1:nband){
    mat1[,i]<-getValues(raster(name1, band=i))
  }
  
  # Stockage des bandes du raster 2 dans un matrice
  for (l in 1:nband){
    mat2[,l]<-getValues(raster(name2, band=l))
  }
  
  ## Création du graphique ##
  out_name<-paste(out_name, ".pdf",sep="")
  pdf(out_name, height=10,width=10)
  ye<-seq(2000,2010,1) # Séquence pour nommer les graphiques, ici, 1 bande= 1 année
  par(mfrow=c(4,3)) # Préparation de la fenêtre graphique : 4 lignes, 3 colonnes pour mettre un graphique par bande dans chaque case

  # Boucle sur les bandes
  for (k in 1:nband){
    df <- data.frame(mat1[,i],mat2[,k])
    x <- densCols(mat1[,k],mat2[,k], colramp=colorRampPalette(c("black", "white")))
    # Création d'une palette pour représenter la concentration des individus
    df$dens <- col2rgb(x)[1,] + 1L
    cols <-  colorRampPalette(c("#000099", "#00FEFF", "#45FE4F","#FCFF00", "#FF9400", "#FF3100"))(256)
    df$col <- cols[df$dens]
    plot(mat2[,k]~mat1[,k], data=df, pch=20, col=col, xlab=name_xlab, ylab=name_ylab, main=ye[k], cex.lab=0.75)
    reg<-lm(mat2[,k]~mat1[,k])
    abline(coef=coef(reg), lwd=1)
  }
  dev.off()
}

Et voici le résultat en sortie pour l'un de mes fichiers :

Result

Carte de Flux sous R

Dans le cadre de mes travaux actuels, j'avais besoin d'avoir une idée des relation de l'Afrique de l'Ouest avec le reste du monde en termes d'échanges alimentaires. Plutôt que de faire un simple tableau avec les données d'export/import, je me suis dit "faisons une carte des flux avec des liaisons de taille proportionnelle à l'intensité de la variables". Voilà Voilà je me suis donc lancée là dessus …. sur R!

Les données que proviennent de la World Integrated Trade Solution , j'ai récupéré les valeurs des échanges alimentaires (Food Products) pour les imports/exports pour les 10 pays ayant le plus de poids.

Le code a été adapté de http://flowingdata.com/2011/05/11/how-to-map-connections-with-great-circles/

Voici donc mon code R:

# Création d'une carte de flux sous R

# Chargement des Packages
library(maps)
library(geosphere)

#Initialisation de l'espace de travail
setwd("D:/MyFolder")

#Importation des données
##un tableau comportant au moins 4 colonnes et n lignes correspondant aux n connections
##latitude du point d'origine
##longitude du point d'origine
##latitude du point de destination
##longitude du point de destination
## + une colonne supplémentaire dans le cas où 
##l'on souhaite des traits de taille proportionelle
##à la variable représentée

data<-read.csv("Export.csv", sep=";", dec=".", header=TRUE)

#Importation du fond de carte et mise en forme
map("world", fill=FALSE, boundary=TRUE,bg="white", lwd=1.5, col="lightslategrey")
map.axes()
map.scale(120,-70,ratio=FALSE, relwidth=0.10, cex.main=2)

#Ajout des lignes
##Création d'une palette de couleur

pal <- colorRampPalette(c("royalblue4", "royalblue"))
colors <- pal(10)

##ajout des lignes/connections
for (i in 1:10){
  lat_o<-data[i,1]
  long_o<-data[i,2]
  lat_d<-data[i,3]
  long_d<-data[i,4]
  val<-data[i,5]
  inter <- gcIntermediate(c(long_o, lat_o), c(long_d, lat_d), n=1000, addStartEnd=TRUE)
  lines(inter, col=colors[i],lwd=val)
}

## Finalisation de la mise en page
title("Exports Food Products for West Africa in 2011",font=3, family="serif")
legend("bottomleft",lty=1, col="royalblue",lwd=1,legend="60760 US$ Thousand",cex=1.25)

Et voici mon résultat pour les exports:

Rplot_Trade_Export

Raster et « Jolie » matrice de corrélation sur R et ggplot, petit tuto

Tout est partis de la rédaction d'un article pour une revue de télédétection dans laquelle j'avais une matrice de corrélation à partir d'un raster multibandes à présenter! Bon une matrice de corrélation c'est pas franchement "sexy" à présenter donc je me suis dit que je pouvais grâce à R et toutes les possibilités qu'il offre réussir à faire quelque chose qui change de l'ordinaire!

Dans mon cas, il s'agissait de vérifier la non-stabilité (c'est  à dire la potentielle dynamique) du domaine cultivé en Afrique de l'Ouest à partir du produit MODIS land Cover! Initialement j'ai donc une image avec 11 bandes binaires (culture / non culture, entre 2001 et 2011), présentées ci dessous:

11 bandes (2001-2011) MODIS affichées dans R

Voici donc la démarche que j'ai suivi :

########################################
# Initialisation de l'espace de travail

setwd("D:/MyFolder")

#Packages nécéssaires
library(raster)
library(sp)
library(Hmisc)
library(plyr)
library(reshape2)
library(ggplot2)

#############Importation des données raster et stockage des bandes dans une matrice#######
filename<-'mon_image.tif'# variable stockant le nom du fichier
nband<-nbands(raster(filename)) #variable stockant le nombre de bandes dans l'image
cell<-ncell(raster(filename)) #variable stockant le nombre de pixel d'une bande
mat<-matrix(0, ncol=nband, nrow=cell) #initialisation de la matrice reçevant les valeurs des bandes de l'image

for (i in 1:nband){
  mat[,i]<-getValues(raster(filename, band=i))
}

#############Création de la matrice de corrélation#######
cor<-rcorr(mat,type="pearson") # création de la matrice de corrélation
coef.corr<-cor$r  # récupère la matrice de coefficients de corrélation

#############Mise en forme de la matrice de corrélation pour la création du graphique#######
coef.corr[coef.corr==1]<-0 #valeurs de la diagonale mises à 0

# Boucle garder que la première moitié du tableau de corrélation (éviter la duplication dans le graphique)
for (i in 1:nband){
  for (j in 1:nband){
    if (coef.corr[i,j]==0){
      coef.corr[i:nband,j]<- 0
    }
  }
}

coef.corr[coef.corr==0]<-NA
mat.m<-melt(coef.corr, na.rm=FALSE)
mat.n<-ddply(mat.m, .(Var1), transform)
r<-round(mat.n$value,2) #arrondissement des coefficients de corrélation à deux chiffres après la virgule
mat.r<-cbind(mat.n[,1:2],r)

#############Diagramme à deux dimensions avec ggplot2#######
labelsx<-seq(2001,2011,1) #noms de l'axe des x
labelsy<-seq(2001,2011,1) #noms de l'axe des y

#Graphique

ggplot(mat.r, aes(factor(Var1), factor(Var2), label=r))+ geom_tile(aes(fill=r ),colour="white")+scale_fill_gradient(low="seagreen4", high="violetred",       na.value="white")+geom_text(aes(fontface=2))+ylab("Years")+xlab("Years")+labs(fill="Correlation \nCoefficient (R)")+labs(title="Heatmap of Correlation Coefficient")+theme(axis.title.x=element_text(size=20), axis.title.y=element_text(size=20),axis.text.x=element_text(size=10), axis.text.y=element_text(size=10),plot.title=element_text(size=20))+theme_bw(base_size=20, base_family=2)+scale_x_discrete(labels=labelsx)+scale_y_discrete(labels=labelsy)+coord_flip()

########################################

Et voilà, vous obtiendrez un chouette graphique que celui-ci:

Diagramme deux dimensions R

ChangeMatters : Interface de Visualisation des dynamiques de la végétation à partir de l’imagerie landsat

Voici une application découverte très récemment et qui s'avère être vraimenent intéréssante. Elle permet de visualiser la végétation à deux dates et de donner une carte de dynamique entre ces deux dates, le tout basé sur de l'imagerie Landsat!

Pour cela, il suffit juste :

  • d'entrer le lieu désiré
  • de spécifier le type de visualisation voulu (infra rouge, indice de végétation, couleur naturelle…)
  • et la période souhaitée

L'interface de visualisation est ensuite très bien faite :

ChangeMatter

Copyright © 2019 DATA\WAX

Theme by Anders NorenUp ↑