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    1
res(C2001)
[1] 30.1 29.6


c) 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] TRUE
res(C2001) == res(C2016)
[1] TRUE TRUE
st_crs(C2001) == st_crs(C2016)
[1] TRUE


Question 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 95


Leur nombre est de:

dim(levels(R2001f)[[1]])[1]
[1] 14


Question 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] 20


Nous 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 Wetlands


Nous attachons maintenant la table d’attributs rat au raster R2001f en utilisant donc la fonction levels():

levels(R2001f) <- rat

R2001f
class      : 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 Wetlands


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

R2016f
class      : 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 Wetlands


Question 4

a) Visualiser le raster catégorique de l’année 2001 avec la fonction mapview(). En particulier:

  • Utiliser la colonne Color du data.frame sommaire pour définir l’argument col.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   1016


Ajoutons 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    1016


a) 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]]$COUNT


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…     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.0324


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.

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




  1. Les données peuvent être téléchargées sur le site du Multi-Resolution Land Characteristics Consortium.↩︎