5.2 Exercices
Dans cette section, vous mettrez en pratique certains concepts vus dans la section leçon de ce module et apprendrez à manipuler des données matricielles catégoriques. Bien que la réponse à chaque question soit disponible, il est très important de tenter d’y répondre par vous-même!
Les questions qui suivent se rapportent à des données matricielles de couvertures terrestres d’une petite région du Colorado aux États-Unis. Ce sont les données utilisées dans le tutoriel de Nikki Inglis et qui proviennent du National Land Cover Database des États-Unis37.
Nous allons également utiliser la bibliothèque FedData qui permet d’extraire certaines données en provenance de sources fédérales américaines. Installer cette bibliothèque avant de débuter les questions:
install.packages("FedData")Question 1
a) Lire les rasters Couvert_2001.tif et Couvert_2016.tif correspondant aux couvertures terrestres lors des années 2001 et 2016 respectivement.
Réponse
library(raster)
C2001 <- raster("Module5/Module5_donnees/Couvert_2001.tif")
C2016 <- raster("Module5/Module5_donnees/Couvert_2016.tif")b) Trouver l’étendue, les dimensions, et la résolution du raster associé à l’année 2001.
Réponse
extent(C2001)class : Extent
xmin : 304442
xmax : 355642
ymin : 4449594
ymax : 4515069 dim(C2001)[1] 2212 1701 1res(C2001)[1] 30.1 29.6c) En utilisant la fonction st_crs() de la library sf, déterminer en format proj4string le SCR utilisé ainsi que l’unité de mesure de la projection du raster de la couverture terrestre de 2001.
Réponse
library(sf)
st_crs(C2001)$proj4string
st_crs(C2001)$units[1] "+proj=utm +zone=13 +datum=WGS84 +units=m +no_defs"[1] "m"d) Confirmer que l’étendue, la résolution et le système de coordonnées sont les mêmes pour les rasters des années 2001 et 2016.
Réponse
extent(C2001) == extent(C2016)[1] TRUEres(C2001) == res(C2016)[1] TRUE TRUEst_crs(C2001) == st_crs(C2016)[1] TRUEQuestion 2
a) Créer deux rasters, R2001 et R2016, en sélectionnant une région carrée de 500 x 500 cellules à l’intérieure des deux rasters de couverture terrestre initiaux.
Réponse
Rappelons que les dimensions des rasters initiaux, C2001 et C2016, sont de 2212 par 1701 cellules. Sélectionnons une région carrée qui soit relativement centrale pour éviter les bordures.
R2001 <- C2001[1000:1499, 700:1199, drop=FALSE]
R2016 <- C2016[1000:1499, 700:1199, drop=FALSE]b) Utiliser la fonction getValues() pour trouver les valeurs comprises dans le raster R2001, et déterminer la classe de ces valeurs.
Réponse
val <- getValues(R2001)
class(val)[1] "numeric"c) Puisque les rasters R2001 et R0016 représentent une couverture terrestre, les valeurs ne devraient pas être numériques mais plutôt catégoriques. Transformer les rasters en facteur en utilisant la fonction as.factor().
Réponse
R2001f <- as.factor(R2001)
R2016f <- as.factor(R2016)d) Déterminer l’identifiant des catégories de couverture terrestre. Combien de catégories différentes exite-t-il?
Réponse
Les niveaux sont donnés par la fonction levels()
levels(R2001f)[[1]]
ID
1 11
2 21
3 22
4 23
5 24
6 41
7 42
8 43
9 52
10 71
11 81
12 82
13 90
14 95Leur nombre est de:
dim(levels(R2001f)[[1]])[1][1] 14Question 3
Pour connaître à quel type de couverts ces catégories correspondent, allons charger la légende disponible dans la bibliothèque FedData, en utilisant la fonction pal_nlcd().
library(FedData)You have loaded FedData v4.
As of FedData v4 we have retired
dependencies on the `sp` and `raster` packages.
All functions in FedData v4 return `terra` (raster)
or `sf` (vector) objects by default, and there may be
other breaking changes.legende <- pal_nlcd()
legende# A tibble: 20 × 4
ID Class Color Description
<dbl> <chr> <chr> <chr>
1 11 Open Water #547… Areas of o…
2 12 Perennial Ice/Snow #FFF… Areas char…
3 21 Developed, Open Space #E8D… Areas with…
4 22 Developed, Low Intensity #E29… Areas with…
5 23 Developed, Medium Intensity #ff0… Areas with…
6 24 Developed High Intensity #B50… Highly dev…
7 31 Barren Land (Rock/Sand/Clay) #D2C… Areas of b…
8 41 Deciduous Forest #85C… Areas domi…
9 42 Evergreen Forest #388… Areas domi…
10 43 Mixed Forest #D4E… Areas domi…
11 51 Dwarf Scrub #AF9… Alaska onl…
12 52 Shrub/Scrub #DCC… Areas domi…
13 71 Grassland/Herbaceous #FDE… Areas domi…
14 72 Sedge/Herbaceous #D1D… Alaska onl…
15 73 Lichens #A3C… Alaska onl…
16 74 Moss #82B… Alaska onl…
17 81 Pasture/Hay #FBF… Areas of g…
18 82 Cultivated Crops #CA9… Areas used…
19 90 Woody Wetlands #C8E… Areas wher…
20 95 Emergent Herbaceous Wetlands #64B… Areas wher…Remarquer que les attributs de ce data.frame:
ID: correspond à l’identifiant de chaque type de couvert;Description: donne une description du type de couvert;Class: correspond à la classe du type de couvert. Par exemple, tous les couverts de forêt (Deciduous Forest, Evergreen Forest, Mixed Forest) sont de classe forest;Color: correspond aux couleurs (en format hexadécimal) à utiliser pour représenter les types de couvert.
La quantité de couverts différents répertoriés dans cette légende est:
dim(legende)[1][1] 20Nous remarquons que le nombre de catégories dans la légende est supérieur au nombre de catégories dans le raster catégorique R2001f.
Créons un tableau sommaire qui reprend les éléments de legende en conservant seulement les types de couverts présents dans le raster R2001f:
type2001 <- unique(R2001f)
sommaire <- legende[legende$ID %in% type2001, -4]
# oùlegende$ID %in% type2001 est un vecteur logique, TRUE lorsque ID est dans type2001 FALSE sinon
# où -4 retire la 4ième colonne du tableau legende contenant une description détaillée du couvert
sommaire# A tibble: 14 × 3
ID Class Color
<dbl> <chr> <chr>
1 11 Open Water #5475A8
2 21 Developed, Open Space #E8D1D1
3 22 Developed, Low Intensity #E29E8C
4 23 Developed, Medium Intensity #ff0000
5 24 Developed High Intensity #B50000
6 41 Deciduous Forest #85C77E
7 42 Evergreen Forest #38814E
8 43 Mixed Forest #D4E7B0
9 52 Shrub/Scrub #DCCA8F
10 71 Grassland/Herbaceous #FDE9AA
11 81 Pasture/Hay #FBF65D
12 82 Cultivated Crops #CA9146
13 90 Woody Wetlands #C8E6F8
14 95 Emergent Herbaceous Wetlands #64B3D5À partir du nouveau data.frame sommaire, nous allons maintenant créer un data.frame appelé rat qui permettra “d’attacher” ces attributs au raster R2001f.
La première colonne d’un data.frame rat doit toujours se nommer “ID” et lister les identifiants des différentes catégories présentes dans le raster.
# Premiere colonne de rat
rat <- levels(R2001f)[[1]]Nous ajoutons l’attribut Class, spécifiant les types de couvert, au data.frame rat
rat$classe <- sommaire$Class
rat ID classe
1 11 Open Water
2 21 Developed, Open Space
3 22 Developed, Low Intensity
4 23 Developed, Medium Intensity
5 24 Developed High Intensity
6 41 Deciduous Forest
7 42 Evergreen Forest
8 43 Mixed Forest
9 52 Shrub/Scrub
10 71 Grassland/Herbaceous
11 81 Pasture/Hay
12 82 Cultivated Crops
13 90 Woody Wetlands
14 95 Emergent Herbaceous WetlandsNous attachons maintenant la table d’attributs rat au raster R2001f en utilisant donc la fonction levels():
levels(R2001f) <- rat
R2001fclass : RasterLayer
dimensions : 500, 500, 250000 (nrow, ncol, ncell)
resolution : 30.1, 29.6 (x, y)
extent : 325482, 340532, 4470698, 4485498 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs
source : memory
names : layer
values : 11, 95 (min, max)
attributes :
ID classe
from: 11 Open Water
to : 95 Emergent Herbaceous Wetlandsa) Refaite les manipulations précédentes pour le raster de la couverture terrestre de 2016.
Réponse
legende <- pal_nlcd()
type2016 <- unique(R2016f)
sommaire <- legende[legende$ID %in% type2016, -4]
rat <- levels(R2016f)[[1]]
rat$classe <- sommaire$Class
levels(R2016f) <- rat
R2016fclass : RasterLayer
dimensions : 500, 500, 250000 (nrow, ncol, ncell)
resolution : 30.1, 29.6 (x, y)
extent : 325482, 340532, 4470698, 4485498 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs
source : memory
names : layer
values : 11, 95 (min, max)
attributes :
ID classe
from: 11 Open Water
to : 95 Emergent Herbaceous WetlandsQuestion 4
a) Visualiser le raster catégorique de l’année 2001 avec la fonction mapview(). En particulier:
- Utiliser la colonne
Colordu data.framesommairepour définir l’argumentcol.regions. - Utiliser l’argument
att =afin que la légende affiche les types de couvertures.
Réponse
mapview(R2001f, col.regions = sommaire$Color, att = "classe", legend = TRUE)b) Utiliser la fonction latticeview() de la bibliothèque leafsync pour afficher les cartes des deux années côte à côte. N’ajouter pas de légende.
Réponse
library(leafsync)
m2001 <- mapview(R2001f, col.regions = sommaire$Color, legend = FALSE)
m2016 <- mapview(R2016f, col.regions = sommaire$Color, legend = FALSE)
leafsync::latticeview(m2001, m2016, ncol = 2)Question 5
En comparant les deux cartes pour les années 2001 et 2016, nous pouvons observer visuellement des changements dans la couverture terrestre. Cependant, nous désirons faire une comparaison plus quantitative de ces changements.
Pour ce faire, nous pouvons utiliser la fonction ratify() de la bibliothèque raster. Celle-ci possède l’argument count qui permet de compter le nombre de cellules de chaque catégorie^(Noter que la fonction ratify()permet également de transformer un raster numérique en raster catégorique: voir la documentation sur cette fonction.). Par exemple:
R2001_nbr<-ratify(R2001f, count = TRUE)
levels(R2001_nbr)[[1]]
ID COUNT
1 11 925
2 21 2435
3 22 653
4 23 43
5 24 3
6 41 68050
7 42 2762
8 43 1083
9 52 149858
10 71 3586
11 81 15511
12 82 1512
13 90 2563
14 95 1016Ajoutons ce calcul du nombre de cellules de chaque type de couvert au data.frame sommaire:
sommaire$Nbr2001 <- levels(R2001_nbr)[[1]]$COUNT
sommaire# A tibble: 14 × 4
ID Class Color Nbr2001
<dbl> <chr> <chr> <dbl>
1 11 Open Water #5475A8 925
2 21 Developed, Open Space #E8D1D1 2435
3 22 Developed, Low Intensity #E29E8C 653
4 23 Developed, Medium Intensity #ff0000 43
5 24 Developed High Intensity #B50000 3
6 41 Deciduous Forest #85C77E 68050
7 42 Evergreen Forest #38814E 2762
8 43 Mixed Forest #D4E7B0 1083
9 52 Shrub/Scrub #DCCA8F 149858
10 71 Grassland/Herbaceous #FDE9AA 3586
11 81 Pasture/Hay #FBF65D 15511
12 82 Cultivated Crops #CA9146 1512
13 90 Woody Wetlands #C8E6F8 2563
14 95 Emergent Herbaceous Wetlands #64B3D5 1016a) Répéter les manipulations précédentes pour l’année 2016.
Réponse
R2016_nbr<-ratify(R2016f, count = TRUE)
sommaire$Nbr2016 <- levels(R2016_nbr)[[1]]$COUNTb) Déterminer le type de couvert qui a le plus augmenté et celui qui a le plus diminué entre l’année 2001 et l’année 2016.
Réponse
Pour chaque type de couvert, nous calculons la différence entre le nombre de cellules en 2016 et le nombre de cellules en 2001. Nous divisons ensuite par le nombre de cellules totales dans l’un ou l’autre des rasters (ils ont le même nombre), ce qui permet de comparer les catégories entre elles.
Diff <- 100*(sommaire$Nbr2016 - sommaire$Nbr2001)/ncell(R2016f)
sommaire$Diff2016_2001 <- Diff
sommaire# A tibble: 14 × 6
ID Class Color Nbr2001 Nbr2016 Diff2016_2001
<dbl> <chr> <chr> <dbl> <dbl> <dbl>
1 11 Open Water #547… 925 1213 0.115
2 21 Developed… #E8D… 2435 2435 0
3 22 Developed… #E29… 653 654 0.0004
4 23 Developed… #ff0… 43 44 0.0004
5 24 Developed… #B50… 3 3 0
6 41 Deciduous… #85C… 68050 68038 -0.0048
7 42 Evergreen… #388… 2762 2781 0.0076
8 43 Mixed For… #D4E… 1083 1115 0.0128
9 52 Shrub/Scr… #DCC… 149858 147788 -0.828
10 71 Grassland… #FDE… 3586 3642 0.0224
11 81 Pasture/H… #FBF… 15511 14869 -0.257
12 82 Cultivate… #CA9… 1512 3879 0.947
13 90 Woody Wet… #C8E… 2563 2604 0.0164
14 95 Emergent … #64B… 1016 935 -0.0324Nous pouvons maintenant utiliser les fonctions which.max() et which.min() qui identifient la position, dans un vecteur, des éléments de valeur maximale et minimale respectivement.
sommaire$Class[which.max(sommaire$Diff2016_2001)][1] "Cultivated Crops"sommaire$Class[which.min(sommaire$Diff2016_2001)][1] "Shrub/Scrub"Ce sont donc les champs cultivés qui ont le plus augmenté et les peuplements d’arbustes qui ont le plus diminué au cours de cette période. Notez que votre réponse dépendra de la sélection de la région carrée faite en a).
Les données peuvent être téléchargées sur le site du Multi-Resolution Land Characteristics Consortium.↩︎