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.
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:
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(terra)
C2001 <- rast("Module5/Module5_donnees/Couvert_2001.tif")
C2016 <- rast("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
SpatExtent : 304441.92102709, 355642.02102709, 4449593.57677233, 4515068.77677233 (xmin, xmax, ymin, ymax)[1] 2212 1701 1[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
[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
[1] TRUE[1] TRUE TRUE[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 <- crop(C2001, ext(xFromCol(C2001, 700), xFromCol(C2001, 1199),
yFromRow(C2001, 1499), yFromRow(C2001, 1000)))
R2016 <- crop(C2016, ext(xFromCol(C2016, 700), xFromCol(C2016, 1199),
yFromRow(C2016, 1499), yFromRow(C2016, 1000)))b) Utiliser la fonction values() pour trouver les valeurs comprises dans le raster R2001, et déterminer la classe de ces valeurs.
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().
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()
[[1]]
ID Couvert_2001
1 11 11
2 21 21
3 22 22
4 23 23
5 24 24
6 41 41
7 42 42
8 43 43
9 52 52
10 71 71
11 81 81
12 82 82
13 90 90
14 95 95Leur nombre est de:
[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().
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.# 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:
[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)[[1]]
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.
Nous ajoutons l’attribut Class, spécifiant les types de couvert, au data.frame rat
ID Couvert_2001 classe
1 11 11 Open Water
2 21 21 Developed, Open Space
3 22 22 Developed, Low Intensity
4 23 23 Developed, Medium Intensity
5 24 24 Developed High Intensity
6 41 41 Deciduous Forest
7 42 42 Evergreen Forest
8 43 43 Mixed Forest
9 52 52 Shrub/Scrub
10 71 71 Grassland/Herbaceous
11 81 81 Pasture/Hay
12 82 82 Cultivated Crops
13 90 90 Woody Wetlands
14 95 95 Emergent Herbaceous WetlandsNous attachons maintenant la table d’attributs rat au raster R2001f en utilisant donc la fonction levels():
class : SpatRaster
size : 500, 498, 1 (nrow, ncol, nlyr)
resolution : 30.1, 29.6 (x, y)
extent : 325511.9, 340501.7, 4470698, 4485498 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 13N (EPSG:32613)
source(s) : memory
varname : Couvert_2001
categories : Couvert_2001, classe
name : Couvert_2001
min value : 11
max value : 95a) Refaite les manipulations précédentes pour le raster de la couverture terrestre de 2016.
Réponse
legende <- pal_nlcd()
type2016 <- unique(R2016f)[[1]]
sommaire <- legende[legende$ID %in% type2016, -4]
rat <- levels(R2016f)[[1]]
rat$classe <- sommaire$Class
levels(R2016f) <- rat
R2016fclass : SpatRaster
size : 500, 498, 1 (nrow, ncol, nlyr)
resolution : 30.1, 29.6 (x, y)
extent : 325511.9, 340501.7, 4470698, 4485498 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 13N (EPSG:32613)
source(s) : memory
varname : Couvert_2016
categories : Couvert_2016, classe
name : Couvert_2016
min value : 11
max value : 95Question 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
Warning: Found more colors (15) than zcol values (14)!
Trimming colors to match number of zcol values.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 freq() de la bibliothèque terra. Celle-ci compte le nombre de cellules de chaque catégorie et retourne un data.frame avec les colonnes layer, value et count. Par exemple:
layer value count
1 1 11 924
2 1 21 2419
3 1 22 652
4 1 23 43
5 1 24 3
6 1 41 67628
7 1 42 2740
8 1 43 1059
9 1 52 149398
10 1 71 3576
11 1 81 15472
12 1 82 1512
13 1 90 2558
14 1 95 1016Ajoutons ce calcul du nombre de cellules de chaque type de couvert au data.frame sommaire:
# A tibble: 14 × 4
ID Class Color Nbr2001
<dbl> <chr> <chr> <dbl>
1 11 Open Water #5475A8 924
2 21 Developed, Open Space #E8D1D1 2419
3 22 Developed, Low Intensity #E29E8C 652
4 23 Developed, Medium Intensity #ff0000 43
5 24 Developed High Intensity #B50000 3
6 41 Deciduous Forest #85C77E 67628
7 42 Evergreen Forest #38814E 2740
8 43 Mixed Forest #D4E7B0 1059
9 52 Shrub/Scrub #DCCA8F 149398
10 71 Grassland/Herbaceous #FDE9AA 3576
11 81 Pasture/Hay #FBF65D 15472
12 82 Cultivated Crops #CA9146 1512
13 90 Woody Wetlands #C8E6F8 2558
14 95 Emergent Herbaceous Wetlands #64B3D5 1016a) Répéter les manipulations précédentes pour l’année 2016.
b) 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… 924 1212 0.116
2 21 Developed… #E8D… 2419 2419 0
3 22 Developed… #E29… 652 653 0.000402
4 23 Developed… #ff0… 43 44 0.000402
5 24 Developed… #B50… 3 3 0
6 41 Deciduous… #85C… 67628 67616 -0.00482
7 42 Evergreen… #388… 2740 2759 0.00763
8 43 Mixed For… #D4E… 1059 1091 0.0129
9 52 Shrub/Scr… #DCC… 149398 147328 -0.831
10 71 Grassland… #FDE… 3576 3632 0.0225
11 81 Pasture/H… #FBF… 15472 14830 -0.258
12 82 Cultivate… #CA9… 1512 3879 0.951
13 90 Woody Wet… #C8E… 2558 2599 0.0165
14 95 Emergent … #64B… 1016 935 -0.0325 Nous 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.
[1] "Cultivated Crops"[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).