Agriculture de précision

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.

Remarque importante : la fonction res() permet d'avoir la dimension en mètres des pixels. 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.

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

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

Convertissons maintenant toutes les donnée des deux rasters en un tableau data.frame. 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.

# Mise au format tableau et fusion

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

df2 <- as.data.frame(r1_ortho,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"))

df <- na.omit(df) # suppression des lignes sans données

colnames(df)

df[-ncol(df)]-> df # étape optionnelle : suppression de la dernière colonne (inutile : masque de la zone)

colnames(df)[c(14,15,16)] <- c("Red_ortho","Green_ortho","Blue_ortho")

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)  

Le plus simple est toujours de travailler avec les données en tableurs.

# Reprenons la data.frame df ci-dessus.

colnames(df)

col_vari <- vari(df$Red,df$Green,df$Blue)

col_vari2 <- vari(df$Red_ortho,df$Green_ortho,df$Blue_ortho)

col_tgi <- tgi(df$Red,df$Green,df$Blue)

col_tgi2 <- tgi(df$Red_ortho,df$Green_ortho,df$Blue_ortho)

col_ndwi <- ndwi(df$Green,df$NIR)

col_ndwi2 <- ndwi(df$Green_ortho,df$NIR*255/max(df$NIR))

df_indices <- data.frame(col_vari,col_vari2,col_tgi,col_tgi2,col_ndwi,col_ndwi2)

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)​​. On peut aussi inclure le NDVI2 certainement plus pertinent.

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)​​.

apply(df_indices,2,function(x){any(is.na(x))}) # Vérifier qu'une donnée manquante ne s'est pas glissée dans le tableau

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

library(KefiR)  # Pour installer KefiR

pareto(myacp$sdev)

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

library(rgl)

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

On peut aussi transformer les résultats de l'ACP pour reproduire une image.

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.

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.

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

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))

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"))))