Manipuler les rasters avec R et le package Terra

GEO UNIV’R Tunisie 2024

Note

Ce support est inspiré les informations de ces différentes sources : Géomatique avec R (Giraud et Pecout 2024), Analyse d’image raster (et télédétection) (Madelin 2021) et Image processing and all things raster (Nowosad 2021). Les exemples proposés sont adaptés au contexte tunisien.

Packages utilisés

• sf : importer, manipuler et exporter des données géographiques vectorielles.
• terra : importer, manipuler et exporter des fichiers rasters.

Si vous n’avez pas ces packages, vous pouvez les installer en exécutant la ligne suivante dans la console.

install.packages(c( 'sf', 'terra'))

Importer des données rasters avec terra

Activer le package terra

library("terra")

Importer des objets vecteur

La ligne qui servira au transect est un objet vectoriel. On peut la charger en tant que tel avec la fonction vect() de terra. L’objet chargé est de classe SpatVector.

Il est aussi possible de la charger avec le package sf conçu pour manipuler les objets vectoriels en priorité. Pour cela, voir le module de l’école thématique portant sur les vecteurs : manipuler des objets vectoriels avec sf

Ici nous chargeons un objet vectoriel au format .shp

# Import du fichier de ligne
Kairouan_Mahdia <- vect("Data/Formes/Kairouan_Mahdia.shp") 

class(Kairouan_Mahdia)

[1] "SpatVector"
attr(,"package")
[1] "terra"

Il est aussi possible de charger des objets au format .gpkg. Nous chargeons le masque qui servira aux prochaines opérations.

# Charger un gpkg
zone_etude <- vect("Data/Formes/Sebkha.gpkg")

Importer des objets raster

La fonction pour importer des objets raster est rast()quelque soit le type d’image.

Importer des MNT

# MNT
MNT_ouest <- rast("Data/DEM/n35_e010_1arc_v3.tif")
MNT_est <- rast("Data/DEM/n35_e011_1arc_v3.tif")

Importer des images satellites en bandes séparées

# Importer les différentes bandes des images sentinels

#La bande bleue
bleu <- rast("Data/Sentinel2/oliviers_2024-02-04_S2_B02.tif")

#La bande verte
vert <- rast("Data/Sentinel2/oliviers_2024-02-04_S2_B03.tif")

#La bande rouge
rouge <- rast("Data/Sentinel2/oliviers_2024-02-04_S2_B04.tif")

#La bande proche infrarouge
proche_infrarouge <- rast("Data/Sentinel2/oliviers_2024-02-04_S2_B08.tif")

Les données sont chargées. À partir de ces bandes, nous pouvons créer une dalle ou une tuile dans le visible.

#On peut créer un objet qui contient déjà les trois bandes dans le visible
image_RGB = c(bleu, vert, rouge)

#Affichage des métadonnées
image_RGB

class       : SpatRaster 
dimensions  : 6988, 8771, 3  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 601840, 689550, 3914010, 3983890  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632) 
sources     : oliviers_2024-02-04_S2_B02.tif  
              oliviers_2024-02-04_S2_B03.tif  
              oliviers_2024-02-04_S2_B04.tif  
names       :   B02,   B03,   B04 
min values  :     1,     1,     1 
max values  : 17368, 16408, 15912 

Importation d’un raster avec plusieurs bandes incluses

RGB_PIR <-  rast("Data/Sentinel2/RGB_PIR.tif")

#Affichage des métadonnées
RGB_PIR

class       : SpatRaster 
dimensions  : 6988, 8771, 4  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 601840, 689550, 3914010, 3983890  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632) 
source      : RGB_PIR.tif 
names       :   B02,   B03,   B04,   B08 
min values  :     1,     1,     1,     1 
max values  : 17368, 16408, 15912, 15448 

Visualiser des images avec terra

Visualiser un raster

On peut faire varier l’affichage en paramétrant l’argument type =de la fonction plot(). Chaque type fait appel à d’autres arguments qui lui sont propres.

Contininous est le mode par défaut, range permet de définir une plage d’affichage, ici, entre les valeurs minimales et maximales Comme dans les logiciels de SIG, ‘minmax()’ doit calculer les statistiques pour certaines taches, il suffit de rajouter le paramètre ‘minmax(image, compute = TRUE)’,

# Avec le type continuous
plot(MNT_ouest, type = "continuous", range = c(-10, max(minmax(MNT_ouest, compute = TRUE)))) 

Interval est pratique pour faire des classes rapides, breaks est le nombre de classes voulues

# avec le type interval
plot(MNT_ouest, type = "interval", breaks = 10)

classes permet d’afficher une classification. Si le raster a déjà été classifier alors il suffit d’indiquer type = 'classes'. Sinon il faut construire cette classification.

Dans notre exemple nous construisons la classification à la volée avec la fonction classify(). Ici, la classification se fait de dix en dix en commençant de zéro à la valeur maximale.

# Avec le type classes et la fonction classify

plot(x = classify(MNT_ouest, seq(0, max(minmax(MNT_ouest)), by = 10)), type = "classes")

Visualiser des images satellites

Visualiser une bande seule

On peut visualiser les bandes une par une avec la fonction plot() en r-base.

Depuis l’objet bleu contenant la bande bleue

# La bande unique contenant le bleu
plot(bleu)

Depuis l’objet RGB_PIR contenant toutes les bandes. Ici on affiche la bande verte en la sélectionnant avec son nom B03

# La bande nommée B03 dans l'image à quatre canaux
plot(RGB_PIR$B03)

Depuis l’objet RGB_PIR on affiche la bande proche infrarouge avec son index

# Le canal numéro 4 de l'image à 4 bande, soit le proche infrarouge ('col = grey(0:100/100)' permet de mettre la bande en noir et blanc)
plot(RGB_PIR[[4]], col = grey(0:100/100))

Visualiser plusieurs bandes

Pour visualiser plusieurs bandes en même temps, on utilise la fonction plotRGB() du package terra.

#On essaye d'afficher l'image

plotRGB(image_RGB)

Error in grDevices::rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = out$alpha, : color intensity 1.29679, not in [0,1]

La visualisation d’un raster multibande nécessite de connaître nos images.

Ce message d’erreur est dû au fait que les valeurs des pixels ne sont pas dans une plage de 0 à 255.

Connaître les valeurs min et max avec la fonction minmac() permet de renseigner la plage d’affichage des pixels.

valeurs_min_et_max <- minmax(image_RGB)
max(valeurs_min_et_max[2,])

[1] 17368

On peut renseigner cette valeur dans l’argument scale = de la fonction plotRGB()

plotRGB(image_RGB, scale = max(valeurs_min_et_max[2,]), stretch="lin")

Il peut être nécessaire de préciser à quelles couleurs correspondent les bandes. Ici il s’agit de l’index des couches dans l’objet ìmage_RGB, mais on peut utiliser leur nom.

plotRGB(image_RGB, r = 3, g = 2, b = 1,  scale = max(valeurs_min_et_max[2,]), stretch="lin")

Visualiser une image multibande depuis des objets distincts

# Pour afficher une image en créant un raster multibande pour la commande même
plotRGB(x = c(rouge, vert, bleu),  scale = max(minmax(c(bleu, vert, rouge))[2,]), stretch = "lin")

Visualiser une image multibande depuis un objet contenant toutes les bandes.

Ici on précise la couleur de chaque bande pour obtenir une image en vraies couleurs.

# Pour afficher le spectre du visible à partir de l'image multicanal
plotRGB(RGB_PIR, r = "B04", g="B03", b="B02",  scale= max(minmax(RGB_PIR)[2,]), stretch="lin")

Ici en fausse couleur

# En fausse couleur
plotRGB(RGB_PIR, r = "B08", g="B04", b="B03",  scale = max(minmax(RGB_PIR)[2,]), stretch="lin")

Fusionner, reprojeter et découper des dalles

L’objectif de cet exercice est de parvenir à la réalisation d’un profil topographique entre les villes de Kairouan à Mahdia à partir des deux tuiles MNT.

Fusionner deux tuiles

Nous allons fusionner les deux tuiles chargées dans la Section 1.3.1

par(mfrow = c(1, 2)) # permet de jouxter deux figures (une ligne, deux colonnes). 
  
plot(MNT_ouest)
plot(MNT_est)

Dans un script r il faut utiliser la fonction dev.off() pour remettre les paramètres d’affichage à zéro. Ici ceci permet de revenir à l’affichage d’une seule figure à la fois.

Pour fusionner deux tuiles, on utilise la fonction merge()

# La fusion est possible avec la fonction merge.
MNT <- merge(MNT_ouest, MNT_est)

# Affichage des métadonnées
MNT

class       : SpatRaster 
dimensions  : 3601, 7201, 1  (nrow, ncol, nlyr)
resolution  : 0.0002777778, 0.0002777778  (x, y)
extent      : 9.999861, 12.00014, 34.99986, 36.00014  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
varname     : n35_e010_1arc_v3 
name        : n35_e010_1arc_v3 
min value   :              -50 
max value   :              235 

Regardons le résultat

# Regardons le résultat
plot(MNT)

Agréger, ou comment baisser la résolution d’un raster

Les images peuvent être lourdes pour nos machines, pour de nombreuses raisons, il peut être utile de baisser la résolution de l’image. Pour cela, la fonction ‘aggregate()’ va nous aider. Le paramètre ‘fact = n’ renseigne le nombre de pixels que l’on veut en horizontal/vertival. ‘fun =’ permet de préciser la logique par lequel vont être agrégées les valeurs des pixels (mean, max, min, median, sum, etc.)

Illustration des fonctions ‘aggregate()’ et ‘disagg()’

#Changeons la résolution du MNT avec un groupement de 100 x 100 par pixel
MNT_allege <- aggregate(MNT, fact = 100, fun = 'mean')
plot(MNT_allege)

Désagréger ou augmenter la résolution d’un raster

Inversement, lorsque des traitements impliquent plusieurs rasters, ils peuvent nécessiter que les rasters soient à la même résolution, il est possible d’augmenter la résolution avec la fonction ‘disagg()’. Comme la fonction ‘aggregate()’, il faut renseigner le facteur, le nombre de pixels voulus par pixel en horizontal/vertival. ‘method’ renseigne la méthode de désagrégation (near ou bilinear)

#Changeons la résolution pour 2x2 pixels pour chaque pixel
MNT_enrichi <- disagg(MNT, fact = 2, method = 'bilinear')

|---------|---------|---------|---------|
=========================================
                                          

Pour connaître la résolution des pixels, il est possible de regarder les métadonnées ou d’utiliser directement la fonction ‘res()’. L’unité est celle utilisée dans le système de coordonnées.

#Résolution des MNT
res(MNT)

[1] 0.0002777778 0.0002777778

res(MNT_allege)

[1] 0.02777778 0.02777778

res(MNT_enrichi)

[1] 0.0001388889 0.0001388889

#Résolution de la bande verte de sentinel
res(bleu)

[1] 10 10

Reprojeter

Visualisons les deux objets qui serviront au transect, le MNT et la ligne vectorielle chargée plus haut.

#Affichons le MNT avec la ligne cette ligne 
plot(MNT)
plot(Kairouan_Mahdia, add = T) # add = T permet d'ajouter la ligne à un plot déjà ouvert

La ligne ne s’affiche pas : vérifions si les projections sont identiques avec la fonction same.crs()

same.crs(MNT, Kairouan_Mahdia)

[1] FALSE

Les objets n’ont pas la même projection, pour les connaître regardons leurs métadonnées

MNT

class       : SpatRaster 
dimensions  : 3601, 7201, 1  (nrow, ncol, nlyr)
resolution  : 0.0002777778, 0.0002777778  (x, y)
extent      : 9.999861, 12.00014, 34.99986, 36.00014  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
varname     : n35_e010_1arc_v3 
name        : n35_e010_1arc_v3 
min value   :              -50 
max value   :              235 

Kairouan_Mahdia

 class       : SpatVector 
 geometry    : lines 
 dimensions  : 1, 1  (geometries, attributes)
 extent      : 599872.4, 689145.4, 3930644, 3949245  (xmin, xmax, ymin, ymax)
 source      : Kairouan_Mahdia.shp
 coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632) 
 names       :    id
 type        : <int>
 values      :    NA

Le raster est en WGS 84 EPSG:4326, et la ligne en WGS 84 EPSG:32632.

Reprojetons le raster pour qu’il corresponde à la ligne avec la fonction project() et la fonction crs() permettant de récupérer le crs d’un objet.

MNT <- project(x = MNT, crs(Kairouan_Mahdia))
  
# Réessayons d'afficher la ligne sur le MNT
plot(MNT)
plot(Kairouan_Mahdia, add = T)

Découper un raster

Il est possible de découper un raster de deux façons, soit avec la fonction crop() qui prend en compte l’emprise (la bbox) d’un objet, soit avec la fonction mask() qui découpe précisément selon la forme de l’objet.

Pour notre exemple nous allons utiliser la fonction crop() pour découper le MNT à l’emprise de la ligne.

MNT_Kairouan_Mahdia <- crop(MNT, Kairouan_Mahdia)

plot(MNT_Kairouan_Mahdia)
plot(Kairouan_Mahdia, add = T)

On peut aussi découper un raster selon un masque chargé dans la Section 1.2. Pour cela on utilise la fonction mask()

MNT_mask <- mask(MNT, zone_etude)

plot(MNT_mask)
plot(zone_etude, add = TRUE, border = "red", lwd = 4)

Le zone_etude ne réduit pas l’emprise de la zone d’étude, mais remplace les valeurs en dehors du masque pas NA. Pour obtenir une zone d’étude à la bonne emprise et découper, on combine les deux fonctions crop() et mask().

MNT_mask <- crop(MNT_mask, zone_etude)

plot(MNT_mask)
plot(zone_etude, add = TRUE, border = "red", lwd = 4)

Intersecter un vecteur et un raster

On peut extraire les valeurs d’un raster en fonction de la forme d’un vecteur. Ici on extrait la valeur des cellules du raster se trouvant sur la ligne.

Kairouan_Mahdia_points <- extract(x = MNT_Kairouan_Mahdia, y = Kairouan_Mahdia)

Pour savoir le type d’objet que nous avons, ‘class()’ permet de nous aider.

# Interrogation du type d'objet
class(Kairouan_Mahdia_points)

[1] "data.frame"

L’objet récupéré est de type data.frame. ‘head()’ permet d’afficher les premières lignes d’une data.frame.

# Visualisation des premières lignes de la df
head(Kairouan_Mahdia_points)

  ID n35_e010_1arc_v3
1  1         64.84071
2  1         64.19274
3  1         62.45517
4  1         60.80166
5  1         60.09102
6  1         60.59118

Pour connaître le nombre de lignes d’une df, la fonction ‘rnow()’ est faite pour nous aider.

nrow(Kairouan_Mahdia_points)

[1] 4080

Note utile : mesurer la longueur des différents objets

Suivant le type d’objet (connaissable avec la fonction ‘classe()’), voici comment les mesurer :

‘ncol()’ : permet de connaître le nombre de colonnes dans une table

# Nombre de colonnes 
ncol(MNT)

[1] 6952

ncol(Kairouan_Mahdia_points)

[1] 2

‘nrow()’ : permet de connaître le nombre de lignes dans une table

#Nombre de ligne
nrow(MNT)

[1] 4339

nrow(Kairouan_Mahdia_points)

[1] 4080

‘lenght()’ : permet de connaître le nombre d’objets dans un vecteur (c())

# Nombre d'objets dans un vecteur
length(Kairouan_Mahdia_points$n35_e010_1arc_v3)

[1] 4080

length(c(4, 8, 9, 4, 2, 3))

[1] 6

‘nlyr()’ : permet de connaître le nombre de couches que contient un SpatRaster

# Nombre de couches dans un objet SpatRaster
nlyr(RGB_PIR)

[1] 4

‘ncell()’ : permet de connaître le nombre de pixels dans un raster

# Nombre de pixels
ncell(MNT)

[1] 30164728

On a intersecté 4080 cellules. Visualisons notre coupe topographique

# Affichons le MNT et le profil sur un même plot
par(mfrow = c(2, 1)) # permet de jouxter deux figures (deux lignes, une colonne).

plot(MNT_Kairouan_Mahdia)
plot(Kairouan_Mahdia, add = T)
title("Zone de la coupe")

plot(x = 1:nrow(Kairouan_Mahdia_points), y = Kairouan_Mahdia_points$n35_e010_1arc_v3, type = 'line')
title("Profil topographique")

Réaliser la coupe avec SF

Réaliser la coupe avec sf permet de définir un nombre de points à intersecter.

# La méthode pour le profil consiste à créer des points réguliers le long de la ligne, d'extraire les valeurs de pixel du raster pour chaque point et créer un graphique. Nous allons utiliser le package sf pour créer les points et donc, de convertir l'objet SpatVector en objet sf.
library(sf)

# st_as_sf() pour la conversion vers un objet sf
Kairouan_Mahdia_points <- st_as_sf(st_sample(x = st_as_sf(Kairouan_Mahdia), type = "regular", size = 1000), crs = 32632)

# La fonction extract permet d'extraire les valeurs du raster
Kairouan_Mahdia_points <- extract(x = MNT_Kairouan_Mahdia, y = Kairouan_Mahdia_points)

# Voyons à quoi ressemblent les données avec les premières lignes avec la fonction head()
head(Kairouan_Mahdia_points)

# Affichons le profil avec plot
plot(x = Kairouan_Mahdia_points$n35_e010_1arc_v3, y = 1:1000) 

# Il est nécessaire de préciser les deux axes et il est nécessaire que les deux séries aient le même nombre de valeurs. La fonction nrow() permet de connaître le nombre de lignes dans un tableau. Changeons les points par des lignes
plot(x = 1:nrow(Kairouan_Mahdia_points), y = Kairouan_Mahdia_points$n35_e010_1arc_v3, type = 'line')

  
# Affichons le MNT et le profil sur un même plot
par(mfrow = c(2, 1)) # permet de jouxter deux figures (deux lignes, une colonne).
  
plot(MNT_Kairouan_Mahdia)
plot(Kairouan_Mahdia, add = T)
plot(x = 1:nrow(Kairouan_Mahdia_points), y = Kairouan_Mahdia_points$n35_e010_1arc_v3, type = 'line')

Autres opérations sur la zone d’étude

Agrégation desagragatio

segragate

Algebre spatial

On peut manipuler les valeurs des rasters selon quatre grands types d’opérations d’algèbre spatiale : - Locale : chaque pixel indépendamment, par exemple calculer un NDVI - Focale : la valeur du pixel dépend de la valeur des pixels l’entourant - Globale et Zonale : la valeur du pixel dépend de la superficie de son espace de référence. Lorsque l’opération est globale, il s’agit de l’emprise de l’image, lorsque l’opération est zonale, il s’agit d’une partie de l’emprise de l’image.

Opération globale

Les opérations globales vont interroger l’ensemble des valeurs d’un raster pour obtenir un résultat, il s’agit souvent d’opération statistique.

‘terra’ regroupe quelque opération dans la fonction ‘global()’

Essayons de mieux connaître nos images

Essayons de connaître l’altitude moyenne

global(MNT, fun = 'mean', na.rm = TRUE)

                     mean
n35_e010_1arc_v3 32.32533

global(MNT, fun = 'max', na.rm = TRUE)

                      max
n35_e010_1arc_v3 234.9012

global(MNT, fun = 'min', na.rm = TRUE)

                       min
n35_e010_1arc_v3 -31.23336

global(MNT, fun = 'range', na.rm = TRUE)

                       min      max
n35_e010_1arc_v3 -31.23336 234.9012

Essayons d’afficher un histogramme de notre MNT

hist(MNT)

Essayons d’afficher un histogramme de notre image satellite avec les quatre bandes

density(bleu, col = "blue")

density(vert, col = "green", add = TRUE)

density(rouge, col = "red")

density(proche_infrarouge, col = "purple")

#plot(x = MNT, lines(v) )
#plot(y = c(freq(round(MNT, 0))$count), x = c(freq(round(MNT, 0))$value), type = "line") 
#pixel_max <- max(minmax(bleu),minmax(vert),minmax(rouge),minmax(proche_infrarouge))
plot(y = c(freq(round(bleu, -2))$count), x = c(freq(round(bleu, -2))$value), type = "line", col = "blue", xlim = c(0,6000), ylab = "nombre de pixel", xlab = "Valeur des pixels") 
lines(y = c(freq(round(vert, -2))$count), x = c(freq(round(vert, -2))$value), col = "green") 
lines(y = c(freq(round(rouge, -2))$count), x = c(freq(round(rouge, -2))$value), col = "red") 
lines(y = c(freq(round(proche_infrarouge, -2))$count), x = c(freq(round(proche_infrarouge, -2))$value), col = "purple") 

Opérations locales

Isolons les surfaces en eau du Sebkha Sidi El Hani grâce au NDWI.

Voici notre zone d’étude

plot(MNT_mask)

La première étape est de réduire l’emprise des bandes pour qu’elles correspondent à l’emprise de la zone d’études. Pour cela nous utilisons le crop.

# Nous aurons besoin de la bande 3 (vert), de la bande 8 (Proche infrarouge) et d'un masque pour alléger le calcul
Sebkha_vert <- crop(vert, zone_etude)
Sebkha_Proche_infrarouge <- crop(proche_infrarouge, zone_etude)


par(mfrow = c(1,2))
plot(Sebkha_vert, main = "Vert") 
plot(Sebkha_Proche_infrarouge, main = "Proche infrarouge")

dev.off()

null device 
          1 

Puis on calcule le NDWI

:::callout-tip title=“Pour rappel” Pour Sentinel : NDWI = (Band 3 – Band 8)/(Band 3 + Band 8) Pour Sentinel : NDWI = (vert – proche infrarouge)/(vert + proche infrarouge) :::

# Calculons le NDVI sous le même principe de la calculatrice raster d'Arcmap ou Qgis

NDWI <- (Sebkha_vert - Sebkha_Proche_infrarouge)/(Sebkha_vert + Sebkha_Proche_infrarouge)

plot(NDWI)

Allons plus loin dans l’opération locale, créons un raster binaire entre surfaces humides/en eau et les autres

On crée un histogramme pour visualiser la répartition du calcul, éventuellement faire des seuillages.

# La fonction hist permet de créer un histogramme, ce qui nous permet de définir un seuil
hist(NDWI)
abline(v = 0)

On peut choisir de visualiser le résultat du NDWI de façon binaire : eau / non-eau

# Création d'un raster True / False : toutes les valeurs supérieures à zéro sont TRUE les autres FALSE
NDWI_binaire <- NDWI > 0

head(NDWI_binaire)

    B03
1 FALSE
2 FALSE
3 FALSE
4 FALSE
5 FALSE
6 FALSE

plot(NDWI_binaire )

Mais on peut aussi faire une classification sur la base des plages standards

• 0,2 – 1 : Surface de l’eau
• 0.0 – 0,2 : Haute humidité
• -0,3 – 0.0 : Sécheresse modérée, surfaces non aqueuses
• -1 – -0.3 : Sécheresse, surfaces non aqueuses

NDWI_class <- classify(NDWI, c(minmax(NDWI)[1,], -0.3, 0,0.2, minmax(NDWI)[2,]))

head(NDWI_class)

       B03
1 (-0.3–0]
2 (-0.3–0]
3 (-0.3–0]
4 (-0.3–0]
5 (-0.3–0]
6 (-0.3–0]

plot(NDWI_class, type = "classes", levels = c("Sécheresse, surfaces non aqueuses",
                                              "Sécheresse modérée, surfaces non aqueuses",
                                              "Haute humidité",
                                              "Surface de l’eau"))

Opération focale

Les opérations focales sont des traitements qui, pour chaque pixel, utiliseront les valeurs des pixels l’entourant. La création d’une carte de pente, les calculs de TPI (terrain position index), et tous les outils relatifs à l’hydrologie sont des opérations focales.

Certaines opérations classiques existent déjà dans le package terra et sont accessibles avec la fonction terrain(), c’est le cas du calcul de pente, du TPI ou de l’ombrage.

On indique l’objet sur lequel effectuer l’opération focale, le type d’opération, l’unité de mesure et surtout le nombre de voisins directs à comparer.

Quelques exemples d’autres calculs possibles

par(mfrow= c(2,2))

# Pente
pente <- terrain(x = MNT_mask, v = "slope", unit = "radians", neighbors = 8)
plot(pente)
title("pente")

# Aspect
aspect <- terrain(MNT_mask, "aspect", unit="radians", neighbors=8)
plot(aspect)
title("aspect")

# TPI
plot(terrain(MNT_mask, "TPI", neighbors=8))
title("TPI")

Avec la fonction focal() nous pouvons paramétrer précisément notre opération. L’argument w= spécifie le nombre de pixel distants du pixel local, l’argument fun= prend la valeur d’une fonction de calcul statistique.

Par exemple il est possible d’adoucir le profil d’élévation calculé dans la Section 3.6. D’abord on lisse les valeurs du MNT

#  moyenne focale 
MNT_Kairouan_Mahdia_focal_9 <- focal(x= MNT_Kairouan_Mahdia, w = 9, fun = "mean")
MNT_Kairouan_Mahdia_focal_19 <- focal(x= MNT_Kairouan_Mahdia, w = 19, fun = "mean")

Puis on intersecte le MNT avec la ligne vectorielle

# La fonction extract permet d'extraire les valeurs du raster

Kairouan_Mahdia_points_focal_9 <- extract(x = MNT_Kairouan_Mahdia_focal_9, y = Kairouan_Mahdia)
Kairouan_Mahdia_points_focal_19 <- extract(x = MNT_Kairouan_Mahdia_focal_19, y = Kairouan_Mahdia)

Puis on compare le résultat

par(mfrow = c(3, 1)) # permet de jouxter deux figures (deux lignes, une colonne).

# On reprend le profil original
plot(x = 1:nrow(Kairouan_Mahdia_points), y = Kairouan_Mahdia_points$n35_e010_1arc_v3, type = 'line')
title("Coupe initiale")

#Pour le plot, supprimons les valeurs na, les opérations focales ne peuvent calculer les pixels en bordure d'image
plot(x = 1:nrow(na.omit(Kairouan_Mahdia_points_focal_9)), y = na.omit(Kairouan_Mahdia_points_focal_9$focal_mean), type = 'line')
title("Coupe adoucie : focal = 9")

plot(x = 1:nrow(na.omit(Kairouan_Mahdia_points_focal_19)), y = na.omit(Kairouan_Mahdia_points_focal_19$focal_mean), type = 'line')
title("Coupe adoucie : focal = 19")

dev.off()

null device 
          1 

Représenter le terrain grâce aux opérations focales et locales

On peut combiner les objets construits grâce aux opérations focales et locales précédentes. Ici on construit une représentation du terrain mêlant l’ombrage et le NDWI.

Pour construire l’ombrage on mobilise nos objets pente et aspect dans la fonction shade() à laquelle on rajoute quelques paramètres.

ombrage <- shade(pente, aspect, angle = 40, direction = 300, normalize = FALSE)

# MNT 
plot(MNT_mask)

# ombrage
plot(ombrage, col = grey(0:100/100), alpha = 0.5, add = T)

# NDWI
plot(NDWI_binaire, type = "classe", col = c(rgb(255, 255, 255, alpha = 0, maxColorValue = 255), "skyblue2"), add = T)

En extra : Testons avec plusieurs images

Cet exemple permet de regarder l’évolution de la surface en eau à partir d’images prises à plusieurs dates. Bien que le code fasse appel à une boucle, il n’est pas bien compliqué puisqu’il reprend ce que l’on a vu jusqu’ici.

D’abord on charge les données automatiquement. Ici on fait la liste de tous les fichiers contenant les caractères “B03” dans le dossier Data/Sentinel2 avec la fonction list.files(), on charge les objets de cette liste avec la fonction rast() puis on les découpe avec crop() avec le masque fait plus haut.

On utilise la syntaxe avec les ‘pipe’ |> pour indiquer une suite d’opération, mais on peut aussi les imbriquer.

On répète ensuite l’opération avec “B08”.

# Chargement et préparation de la bande du proche infrarouge (B08) suivant la logique des 'pipe'
B08 <-  list.files("Data/Sentinel2", pattern = "B08", full.names = T) |> 
        rast() |> 
        crop(y = zone_etude)

# Chargement et préparation de la bande du proche infrarouge (B08) suivant la logique des fonctions imbriquées
B03 <-  crop(rast(list.files("Data/Sentinel2", pattern = "B03", full.names = T)), zone_etude)

Puis on réalise l’affichage

# Prépare une figure avec 5 lignes et 3 colonnes
par(mfrow = c(3,3)) 

for(i in 1:nlyr(B03)){ #1:nlyr(B03) permet de créer une suite de nombre commençant de 1 jusqu'au nombre de couches que contient l'objet B03. L'objet `i` aura pour valeur chaque chiffre de la séquence 1:nlyr(B03)
  # Calcul le NDWI pour la couche i
  NDWI_i <- (B03[[i]] - B08[[i]]) / (B03[[i]] + B08[[i]])  
  # Affiche le NDWI pour i
  plot(NDWI_i)
  # Dessine l'histogramme pour i
  hist(NDWI_i)
  # Affiche le raster binaire eau oui / non
  plot(classify(NDWI_i, c(minmax(NDWI_i)[1,], 0, minmax(NDWI_i)[2,])))
  # efface l'objet 
  rm("NDWI_i")
}

Zonale

Exporter un fichier

Vers gpkg, vers tif, vers format sf

Exporter un raster avec ‘SpatRaster()’

writeRaster(MNT_Kairouan_Mahdia_focal_19, "Data/Export/MNT_lisse.tif", overwrite = TRUE)

Exporter un vecteur avec ‘writeVector()’

writeVector(Kairouan_Mahdia, "Data/Export/Ligne.shp", overwrite = TRUE)
writeVector(Kairouan_Mahdia, "Data/Export/Ligne.gpkg", overwrite = TRUE)

Convertir un objet SpatVector vers un objet sf

#Vérifions que notre fichier est bien un objet SpatVector avec la fonction class()
class(zone_etude)

[1] "SpatVector"
attr(,"package")
[1] "terra"

#Activons le package sf

library(sf)

Utilisons la fonction ‘st_as_sf()’ pour convertir le type d’objet et vérifions avec la fonction ‘class’

zone_etude_SF <- st_as_sf(zone_etude)
class(zone_etude_SF)

[1] "sf"         "data.frame"

Regardons si on plot notre MNT avec nos deux objets vecteurs, l’un sf, l’autre SpatVector, s’ils s’affichent et se superposent bien.

#Affichage du MNT en nuances de gris
plot(MNT, col = grey(0:100/100))
#Affichage de notre objet sf en rouge
plot(zone_etude_SF, add = TRUE, border = "red", lwd = 10)
#Affichons notre SpatVector en blanc
plot(zone_etude, add = TRUE, lwd = 4, border = "white")

---

A propos de ce document

Ce support a été créé pour la semaine de formation franco-tunisienne GEO UNIV’R Tunisie 2024 - “Enseigner la statistique, la cartographie et l’analyse spatiale avec R qui se tient à Sousse en mai 2024.

Références

sessionInfo()

R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS 15.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Paris
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] sf_1.0-19    terra_1.7-78

loaded via a namespace (and not attached):
 [1] digest_0.6.37      codetools_0.2-20   fastmap_1.2.0      xfun_0.49         
 [5] magrittr_2.0.3     e1071_1.7-16       KernSmooth_2.23-24 knitr_1.49        
 [9] htmltools_0.5.8.1  rmarkdown_2.29     classInt_0.4-10    cli_3.6.3         
[13] grid_4.4.1         DBI_1.2.3          proxy_0.4-27       class_7.3-22      
[17] compiler_4.4.1     rstudioapi_0.17.1  tools_4.4.1        evaluate_1.0.1    
[21] Rcpp_1.0.13-1      yaml_2.3.10        rlang_1.1.5        jsonlite_1.8.9    
[25] htmlwidgets_1.6.4  units_0.8-5       

Giraud, Timothée, et Hugues Pecout. 2024. Géomatique avec R. https://doi.org/https://doi.org/10.5281/zenodo.5906212.

Madelin, Malika. 2021. « Analyse d’images raster (et télédétection) ». https://mmadelin.github.io/sigr2021/SIGR2021_raster_MM.html.

Nowosad, Jakub. 2021. « Image processing and all things raster ». https://nowosad.github.io/SIGR2021/workshop2/workshop2.html.