Agriculture de précision

Remarque : on peut avoir un message Warning à l'ouverture d'un raster parlant de rectify(). Cela implique de faire une correction du raster :

r0 <- terra::rectify(r0) # Faire une rotation du raster si message d'alterne rectify à l'ouverture

Très important ! Appliquez à chaque raster la fonction names() pour confirmer les longueurs d'onde qu'il présente :

names(r_NIR)

Affichage en couleurs naturelles

# AFFICHAGE EN COULEURS NATURELLES

terra::plotRGB(r_rgb,1,2,3) # 1, 2, 3 ordre des couches rgb - parfois à réordonner selon les drones. On peut avoir ainsi le rouge en deuxième place et le vert en première, il faudrait alors mettre 2,1,3

# Si l'image ne sort pas car valeurs très en-dessous de 255, on peut les multiplier par un nombre empirique :

terra::plotRGB(r_NIR*40000,3,2,1) # ne marche pas car les couches ne sont pas codés de 0 à 255, il faut donc multiplier par 255/max(value) - la valeur de 1810 est donc à adapter au cas par cas selon les drones...

Optionnellement, je peux aussi forcer le raster à raisonner en longitude-latitude pour pouvoir ajouter des points de repères où le placer sur une carte.

# Définir la nouvelle projection

new_crs <- "+proj=longlat +datum=WGS84"

r1_long_lat<- project(r1, new_crs) # terra

Pour se simplifier la vie, on ne va pas le faire ici.

Quelques fonctions  seront ici bien utiles :

• res(), permet de connaître la distance qui correspond à un seul pixel.

• dim(), permet de connaître les dimensions d'une image et le nombre de canaux par pixel

• aggregate(), va permettre d'agréger les valeurs des pixels selon la moyenne ou, mieux, la médiane, pour la compresser par une baisse de la résolution.

• resample(),permet d'aligner la résolution d'images capturées avec différentes caméra (drone, satelitte)

Remarque importante : la fonction res() permet d'avoir la dimension en mètres des pixels. 

res(r1) #  0.01679497 0.01679497 Mes pixels ont une précision de 1,6 cm. Ca va loin...

# Si, par exemple, cela ne sert à rien d'avoir une résolution au-delà de 5 mètres si l'on souhaite faire un traitement personnalisé avec un gros tracteur.

5/res(r0_ortho) # Cela donne le facteur de réduction que je peux appliquer à mon raster (entre 200 et 300 ici) 

# Attention : 5 m n'est qu'un exemple, et plutôt sévère...

Aligner la résolution de toutes les images

Dans cet exemple, l'image satellite a une résolution différente des images r_rgb_comp et r_r_NIR_comp...

dim(r_sat) # 64 79 11 - mesures à 11 lambdas, image de 64 sur 79 pixels.

De la même façon, r_rgb_comp et r_r_NIR_comp n'ont pas la même résolution exactement, il faut donc les aligner en résolution :

# Resample le raster NIR en le forçant à la résolution du raster rgb

r_NIR_comp <- resample(r_NIR_comp, r_rgb_comp, method = "bilinear")

Il ne reste plus qu'à visualiser ou calculer les indices...

Si vos images Pix4Dfields n'ont pas des noms correspondant à ceux indiqués, voici une légende des 11 bandes présentes dans l'ordre : 

"Blue"      "Green"     "Red"       "Red edge"  "Red edge"  "Red edge"  "NIR"       "NarrowNIR" "SWIR"      "SWIR"      "Alpha"  

# VARI

vari<- function(r,g,b) {

  if (is(r)[1] == "SpatRaster") {

vari <- clamp(((g - r) / (g+ r - b)),-1,1)# Version pix4Dfields normalisée entre -1 et 1

              } else {

        vari <- (g - r) / (g+ r - b)  

  }        

return(vari)

}

# TGI

tgi <- function(r,g,b) {

rr <- as.data.frame(r) ; gg <- as.data.frame(g) ; bb <- as.data.frame(b)

# Normaliser le pixel pour normaliser le TGI

tgi <- (g-0.39*r-0.61*b)/max(apply(cbind(rr,gg,bb),2,max))

#bbb <- -0.5*(190*(r-v) - 120*(r-b))#/max(c(r,v,b))

#print(cor.test(as.data.frame(tgi)[,1],as.data.frame(bbb)[,1]))

#plot(as.data.frame(tgi)[,1],as.data.frame(bbb)[,1])

return(tgi)

}

# NDWI

ndwi <- function(g, nir) {

ndwi<- (g - nir)/(g + nir)

return(ndwi)

}

# On notera que certains suggèrent de faire un NDWI/NDMI avec cette formule  (NIR - SWIR) / (NIR + SWIR)  

Remarque importante : il existe d'autres indices peut être nettement plus pertinents, pour en citer certains :

Normalized Difference Vegetation Index (NDVI) : Utilisé pour estimer la biomasse, l'indice foliaire (LAI) et l'activité photosynthétique dans les prairies mixtes et les chaparrals en Californie du Sud, avec une taille de pixel optimale de 6 m ou moins pour des études fonctionnelles (Rahman et al., 2003)​​.

Renormalized Difference Vegetation Index (RDVI) : Comparaison avec d'autres indices pour l'estimation de l'indice foliaire (LAI) dans des écosystèmes de prairies mixtes, montrant une légère amélioration par rapport aux autres indices (He et al., 2006)​​.

Adjusted Transformed Soil-Adjusted Vegetation Index (ATSAVI) : Utilisé pour améliorer l'estimation du LAI dans les prairies, particulièrement lorsqu'il est combiné avec l'indice d'absorption de cellulose (CAI) pour un nouvel indice (L-ATSAVI), augmentant la capacité d'estimation de l'AGB d'environ 10% (He et al., 2006)​​.

Water Band Index (WBI) : employé comme indicateur de contenu en eau des plantes dans les prairies et les chaparrals, avec des analyses de mise à l'échelle spatiale pour optimiser les procédures d'échantillonnage (Rahman et al., 2003)​​.

Photochemical Reflectance Index (PRI) : Utilisé pour les études de la biomasse, du contenu en eau des plantes et de l'activité photosynthétique (Rahman et al., 2003)​​.

Modified Chlorophyll Absorption Ratio Index 2 (MCARI2) : Performant pour estimer le LAI dans les prairies mixtes (He et al., 2006)​​.

Pure Vegetation Index (PVI) : Un index extrait du VI normal pour minimiser les effets du sol, améliorant ainsi les estimations de l'AGB des prairies (Li et al., 2016)​​.

Cellulose Absorption Index (CAI) : Utilisé en combinaison avec ATSAVI pour améliorer les estimations de LAI en tenant compte de la matière morte au sol (He et al., 2006)​​.

Soil Adjusted Total Vegetation Index (SATVI) : Développé pour estimer la couverture végétale totale, la hauteur et la biomasse dans les prairies arides (Marsett et al., 2006)​​.

Et d'autres ?

Convertissons maintenant toutes les données des différents rasters en un tableau data.frame (ici df). Il sera possible sur ce tableau de faire des ACP et autres calculs, mais c'est peut-être plus pertinent de faire directement un calculs des indices.

# Attention : vérifier pour chaque raster qu'ils ont bien la même résolution avant de passer à une fusion

dim(r_rgb)

dim(r_NIR)

dim(r_vari)

# (...)

# Attention : si les données n'ont pas les mêmes dim() ==> resample()

# Mise au format tableau et fusion

# L'option x et y permet d'ajouter au data.frame la position de chaque pixel en ajoutant des colonnes x ety

df1 <- as.data.frame(r_rgb,xy=TRUE)

df2 <- as.data.frame(r_NIR,xy=TRUE)

# Fusionner les dataframes en utilisant les colonnes x et y comme clés de fusion

df <- merge(df1, df2, by = c("x", "y"))

# Si je veux ajouter un indice au tableau :

# 1) Je convertis le raster indice en tableau

df_vari <- as.data.frame(r_vari,xy=TRUE) 

# 2) Je renomme la colonne qui présente l'indice vari (ici la 3e)

colnames(df_vari)[3] <- "VARI"

# 3) J'utilise merge() pour remplacer df par df et la nouvelle colonne

df <- merge(df1, df_vari, by = c("x", "y"))

colnames(df) # pour vérifier que tout est OK

# Et ainsi de suite...

# A la fin, je dégage les données manquantes

df<- na.omit(df)

# Produire df_indices qui ne contient que les variables que l'on souhaite garder .

df_indices <- df[,c("VARI","Red",(...)]

myacp <- prcomp(na.omit(df_indices),center=T,scale=T)

library(KefiR)  # Pour installer KefiR

pareto(myacp$sdev)

biplt(myacp,labels="",cex=0.2)

# Visulasation en 3D

library(rgl)

plot3d(myacp$x[,1:3])

Rien n'empêche de faire de la sélection de variables sur l'ACP. Si j'ai fait un clustering (partie d'après) sur quelques indices, je peux forcer à ce que mon ACP soit faite uniquement sur ces indices.

Visualiser les résultats de l'ACP en fausses couleurs

On peut aussi transformer les résultats de l'ACP pour utiliser les 3 principales pour coloriser le rouge, vert et bleu pour reproduire une image combinant ces informations.

Il faut regrouper les valeurs de l'ACP dans un seul tableau :

# install.packages("tidyterra")

library(tidyterra)

# Récupérer les positons x et y et les agglomérer aux composantes 1 à 3

tablo <- cbind(df[,1:2],myacp$x[,1:3])

Il faut d'abord mettre les résultats à l'échelle car les composantes qui peuvent avoir des valeurs négatives doivent avoir des valeurs entre 0 et 255 pour être reconnues en tant que couleurs.

# Mettre la valeur min de chaque composante à 0 et le max à 255

myraster <- as_spatraster(tablo)

<- apply(tablo[,3:5] ,2,function(x) {(x-min(x))/(max(x)-min(x))*255})

# Voir si la distribution des composantes est bonne

boxplot(table[,3:5])

# Si une composante semble avoir une médiane beaucoup plus haute que les autres, on peut les forcer à avoir la même moyenne

# Exemple : forçons la Composante 1 (colonne 3) à avoir la même moyenne que la composante 2 (colonne 4) sans changer le min et le max.

tablo[,3]<- KefiR::exp_dyn(tablo[,3],mean(tablo[,4],255)

Convertir ensuite le tableau en raster :

myraster <- as_spatraster(tablo)

Revenons plus tard sur le biplt(), on pourra colorer tous les points en fonctions des zones identifiées par clustering.

Faire du regroupement, c'est identifier les zones sans avoir les réponses à l'avance. Il existe de nombreuses techniques :

• k-means, pas toujours performant mais rapide.

• CAH, très lent, assez bon, qui a l'avantage de permettre d'aider à trouver le nombre de types de zones.

• Mélange de gaussienne, très performant, mais qui exige que l'on fixe soi-même le nombre de types de zone.

Attention, le clustering nécessite de centrer-réduire avec scale() si vos indicateurs n'ont pas la même échelle. Pensez à utliser scale si nécessaire.

Pour df_indice (df regroupant que les variables utiles, cf. partie ACP).

Avez-vous pensé à utiliser scale() ? Avez-vous bien sélectionné uniquement les (bons) indicateurs et sorti les autres colonnes telles x et y ?

Le SVM peut s'avérer  trop lourd, ou trop simpliste. Une autre approche, très efficace est de faire glisser une fenêtre de w sur w pixels et de lui attribuer à chaque fois le groupe majoritaire.

r_smooth <- focal(raster_catego, w = 9, fun = modal, na.policy = "omit", pad = TRUE)

Avec w : la taille de la fenêtre (idéalement un nombre impaire).

On voit ainsi que la catégorie 3 présente un indice NDWI (celui déterminé à partir des données du drone) inférieur à -0.66. Cela implique que cette zone présente un fort stress hydrique qui justifiera qu'on lui donne plus d'eau (ou qu'on change d'autres modalités de culture). Au passage, on notera que les zones considérées en catégorie 3 par cet arbre ne le seront effectivement qu'à 89% (11% seront en réalité de la catégorie 1). On peut toujours affiner un tel arbre en faisant des forêts aléatoires qui permettront de ne conserver que les indices les plus influents avant de recommencer le cycle regroupement(catégorisation) des étapes ci-avant.

Pour séparer la catégorie 1 de la catégorie 2, seul l'indice VARI (celui déterminé par drone) a été retenu. Il permet de décrire la zone 1 comme présentant une faible vigueur de végétation. Sans surprise, la catégorie 2 représente de la végétation plus verte et plus vigoureuse.

Il faudrait donc, en théorie, faire maintenant appel à un.e expert.e pour proposer des traitements différenciés.

Les données compilées en raster peuvent être exportées en tiff, recompilées avec pix4Dfield pour permettre d'automatiser le travail d'un tracteur.

Idéalement, on peut aller beaucoup plus loin qu'un arbre de décision simpliste :

en faisant des forêts aléatoires... Cela va donner l'importance  relative de chaque indicateur, c'est donc beaucoup mieux.

On peut aussi directement traiter les données statistiquement et graphiquement sous R

Si on fait voler le drone plusieurs fois dans l'année, il devient possible de suivre les indicateurs dans le temps.

1. Charger tous les TIFF

# Initialiser des vecteurs pour stocker les mois et les années

dates <- c()

months <- numeric(length(tif_files))

years <- numeric(length(tif_files))

days <- numeric(length(tif_files))

# Pour chaque fichier .tif

for (i in seq_along(tif_files)) {

  # Extraire la date du nom du fichier

  file_name <- basename(tif_files[i])

  date_str <- substr(file_name, 1, 10)

  # Convertir la date en objet Date

  date <- as.Date(date_str)

  # Extraire le mois et l'année de la date

  dates <- c(dates,as.character(format(date, "%Y-%m-%d")))

  days[i] <- as.numeric(format(date, "%d"))

  months[i] <- as.numeric(format(date, "%m"))

  years[i] <- as.numeric(format(date, "%Y"))

}

2. Ne conserver que les TIFF de l'année écoulée

# Nettoyer la liste de TIFF du dossier pour ne garder que ceux de l'année en cours

dates <- strptime(dates,format="%Y-%m-%d")

indices <- which(as.Date(dates) >= as.Date("2022-05-31")) # Ici on ne garde que les TIFF ayant dépassé mai 2022

tif_files <- tif_files[indices]

3. Calculer les indices pour chaque TIFF puis faire la médiane pour chaque TIFF

# Initier les listes des valeurs extraites de chaque TIFF

dates <- c()

months <- numeric(length(tif_files))

years <- numeric(length(tif_files))

days <- numeric(length(tif_files))

indice <- numeric(length(tif_files))

# Pour chaque tif de l'année en cours

for (i in seq_along(tif_files)) {

  # Ouvrir le fichier .tif en tant qu'objet raster

  r0 <- terra::rast(tif_files[i])

  r1 <- rectify(r0)

  # Calcule de l'indice du TIFF (ici NDWI)

  indice.temp <- ndwi(r1[["Green"]],r1[["NIR"]]) 

  # Extraire la date du nom du fichier

  file_name <- basename(tif_files[i])

  date_str <- substr(file_name, 1, 10)

  # Convertir la date en objet Date

  date <- as.Date(date_str) 

  # Extraire le mois et l'année de la date

  dates <- c(dates,as.character(format(date, "%Y-%m-%d")))

  days[i] <- as.numeric(format(date, "%d"))

  months[i] <- as.numeric(format(date, "%m"))

  years[i] <- as.numeric(format(date, "%Y"))

  # Calcule de l'indice médian de l'image

  indice[i] <- median(as.data.frame(indice.temp[[1]])[,1],na.rm=T)

}

4. Tracer le graphique

dates <- strptime(dates,format="%Y-%m-%d")

plot(dates,indice,type="o",axes=F,lwd=2) ; grid()

# Abscisses exprimant dans ses graduations uniquement le nom des mois

axis.POSIXct(1, at=seq(from=strptime("01/06/2022","%d/%m/%Y" ), to=strptime("31/05/2023","%d/%m/%Y" ), by="month"), format="%b", las=2) # force indiquant la date de chaque jour

# Ordonnées

axis(2,seq(from = -1, to=1,by=0.1))

Suivre l'évolution climatique implique d'avoir des pratiques culturales constantes sur la région sur au moins cinq ans. Inutile donc de le faire si le terrain étudié change de récoltes.

library(terra)

setwd("G:/Mon Drive/IUT/statistiques/SAE_AGRO/Drone/Export/La Forge/NDVI/La Forge")

tif_files <- list.files(pattern = "Délimitation 2.data.tif$", full.names = TRUE)

# Initialiser des vecteurs pour stocker les mois et les années

months <- numeric(length(tif_files))

years <- numeric(length(tif_files))

days <- numeric(length(tif_files))

indice <- numeric(length(tif_files))

dates <- c()

# Pour chaque fichier .tif

for (i in seq_along(tif_files)) {

  # Ouvrir le fichier .tif en tant qu'objet raster

  r0 <- terra::rast(tif_files[i])

  r1 <- rectify(r0)

  indice.temp <- vari(r1[["Red"]],r1[["Green"]],r1[["Blue"]])

  #indice.temp <- ndwi(r1[["Green"]],r1[["NIR"]])

  # Extraire la date du nom du fichier

  file_name <- basename(tif_files[i])

  date_str <- substr(file_name, 1, 10)

  # Convertir la date en objet Date

  date <- as.Date(date_str)

  # Extraire le mois et l'année de la date

  dates <- c(dates,as.character(format(date, "%Y-%m-%d")))

  days[i] <- as.numeric(format(date, "%d"))

  months[i] <- as.numeric(format(date, "%m"))

  years[i] <- as.numeric(format(date, "%Y"))

  indice[i] <- median(as.data.frame(indice.temp)[,1],na.rm=T)

}

Tracer le graphique et analyser la régression

dates <- strptime(dates,format="%Y-%m-%d") 

synth <- data.frame("dates"=as.numeric(dates),"indice"=indice)

plot(dates,indice,axes=F,data=synth)

reg <- lm(indice~as.numeric(dates))

summary(reg)

abline(reg,col="red",lwd=3)

axis.POSIXct(1, at=seq(from=strptime(dates[1],"%Y-%m-%d" ), to=strptime(dates[length(dates)],"%Y-%m-%d" ), by="year"), format="%Y", las=2) # force indiquant la date de chaque jour

# Ordonnées

axis(2,seq(0,1,by=0.1))

Call:

lm(formula = indice ~ as.numeric(dates))

Residuals:

     Min       1Q   Median       3Q      Max 

-0.26641 -0.11099 -0.04687  0.11155  0.45676 

Coefficients:

                    Estimate Std. Error t value Pr(>|t|)    

(Intercept)       -2.562e+00  7.002e-01  -3.660 0.000599 ***

as.numeric(dates)  1.565e-09  4.312e-10   3.629 0.000659 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1601 on 51 degrees of freedom

Multiple R-squared:  0.2052,    Adjusted R-squared:  0.1896 

F-statistic: 13.17 on 1 and 51 DF,  p-value: 0.0006588

Permettons nous une prédiction. Comment sera l'indice VARI en 2030 ?

Actuellement, lors de cet étude, il est de 0.21 (mai 2023), mais il devrait passer à 0.42 en mai 2030.

Si j'applique le même raisonnement en changeant d'indice, je trouve une évolution moins significative du NDWI qui devrait le faire passer de -0.62 à -0.67. Cela permet de diagnostiquer que, sans surprise, le stress hydrique semble à la hausse, mais que, pour l'instant, la végétation profite de la hausse des températures* (*Attention diagnostic sur seulement cinq ans).

median(df_indices_categories$col_vari) # Mai 2023, mois de l'étude

Prédire sa valeur pour une date.

date_to_predict <- "2030-05-30"

predict(reg,data.frame("dates"=as.numeric(strptime(date_to_predict,format="%Y-%m-%d"))))