library(terra)
# import de l'occupation du sol
occ_sol <- rast('./Data/ESA_WorldCover_2021_fr.tif')
# affichage
plot(occ_sol)
Géotraitements et rasters 1
Objectifs du deuxième volet
Le but de cette deuxième partie est de croiser plusieurs couches vecteurs à l’aide de géotraitements afin d’en tirer des informations supplémentaires. Nous en profiterons également pour faire quelques manipulations basiques sur des fichiers raster. À la fin de cette partie, vous serez capables de faire les manipulations suivantes avec R :
- importer un raster dans R via terra
- exporter un raster
- calculer des statistiques zonales sur un raster par comparaison à une couche vecteur
- reprojeter un raster
- découper et masquer un raster selon une couche de masque
- vectoriser un raster
- automatiser des tâches
Le fil rouge de cette partie sera de déterminer, pour un département donné, la commune ayant le plus de superficie boisée proportionnellement à sa surface.
Données fournies
Pour cet exercice, nous utiliserons les données fournies suivantes :
- la couche des communes dans laquelle nous avons joint les populations (communes_France_pop.gpkg cf partie précédente)
- la couche des départements de France métropolitaine (dep_France.gpkg cf partie précédente)
- la couche raster de l’occupation des sols de 2021 pour une portion de la France (ESA_WorldCover_2021_fr.tif)
La couche d’occupation du sol ESA WorldCover est un raster présentant l’occupation du sol pour le monde entier à une résolution spatiale de 10 m. Pour les besoins de l’exemple, seule un petit extrait pour une partie centre orientale de la France est fourni et cet extrait a été rééchantillonné à 50 m de résolution afin d’alléger le poids du fichier. La nomenclature est la suivante (Tableau 1) :
| ID classe | Classe | Couleur |
|---|---|---|
| 10 | Espace boisé | #006400 |
| 20 | Broussailles | #ffbb22 |
| 30 | Prairie | #ffff4c |
| 40 | Culture | #f096ff |
| 50 | Bâti | #fa0000 |
| 60 | Sol nu / Végétation éparse | #b4b4b4 |
| 70 | Neige et glace | #f0f0f0 |
| 80 | Eau permanente | #0064c8 |
| 90 | Zone humide herbeuse | #0096a0 |
| 95 | Mangrove | #00cf75 |
| 100 | Mousse et lichen | #fae6a0 |
Tableau 1 : Nomenclature ESA World Cover
Importer, explorer et reprojeter un raster
Nous allons commencer par importer notre raster d’occupation du sol dans R via terra. Nous obtiendrons ainsi un objet spatRaster.
Nous pouvons explorer certaines propriétés utiles de cette couche, comme la résolution spatiale et le système de coordonnées de référence (SCR) utilisé.
# extraction de la résolution spatiale du raster (en X et en Y)
res_X <- res(occ_sol)[1]
res_Y <- res(occ_sol)[2]Quelle est la forme des pixels sur notre raster ? Quelles sont les dimensions de cette forme ?
# extraction du SCR avec toutes les informations
scr_rast <- crs(occ_sol, describe=T)
# extraction du code EPSG seulement
epsg_rast <- crs(occ_sol, describe=T)$codePouvez-vous directement croiser cette couche avec votre couche des communes ?
Réponse
Notre raster d’occupation du sol est en WGS84/UTM Zone 31 Nord et notre vecteur des communes en Lambert 93. Nous devons donc au préalable reprojeter notre raster en Lambert 93 (EPSG 2154). Nous spécifions ici une méthode de rééchentillonnage au plus proche voisin (method = ‘near’) car nous avons affaire à un raster discret (catégoriel).
# reprojection du raster
occ_sol <- project(occ_sol, 'EPSG:2154', method = 'near')Découper et masquer un raster
Nous souhaitons maintenant caractériser l’occupation du sol à l’échelle de la Nièvre, dont nous avons extrait les contours dans la partie précédente. Commençons par (ré)importer le fichier vecteur de ce département.
# import du département
dep_58 <- vect('./Data_processed/dep_58.gpkg')Nous allons commencer par découper et masquer notre raster d’occupation du sol selon les contours du département.
Découper et masquer
Découper (clip) et masquer (mask) sont deux opérations différentes et complémentaires. La première permet de découper un raster selon l’étendue d’un polygone. Le raster est donc réduit mais les valeurs des pixels externes au polygone sont encore présentes. La seconde permet de transformer ces pixels externes en no data. D’un point de vue informatique, les pixels sont toujours là, mais ils n’ont plus de valeurs intelligibles.
L’ordre des opérations importe peu. Nous commencerons par découper puis nous masquerons. Afficher les deux rasters pour vous rendre compte.
# découpe du raster d'occupation du sol
occ_sol_58 <- crop(occ_sol, dep_58)
# masque du raster selon les limites du département
occ_sol_58 <- mask(occ_sol_58, dep_58)
plot(occ_sol_58)
Afin de croiser notre occupation du sol avec les autres couches, il est bien de vectoriser ce raster.
Vectoriser un raster
Dans le monde du SIG nous avons l’habitude de travailler principalement avec des données vecteurs. Nous allons ainsi vectoriser notre raster d’occupation du sol. Chaque groupe de pixels homogènes sera intégré au sein d’un polygone dont un attribut prendra la valeur des pixels du groupe homogène, dans notre cas il s’agira de l’identifiant de classe.
# vectorisation du raster d'occupation du sol
occ_sol_dep <- as.polygons(occ_sol_58)
occ_sol_dep class : SpatVector
geometry : polygons
dimensions : 8, 1 (geometries, attributes)
extent : 688357.1, 793143.2, 6616936, 6720873 (xmin, xmax, ymin, ymax)
coord. ref. : RGF93 v1 / Lambert-93 (EPSG:2154)
names : ESA_WorldCover_2021_fr
type : <int>
values : 10
20
30Nous pouvons vérifier que notre nouvelle couche contient un attribut dans lequel nous trouvons l’identifiant des classes d’occupation du sol. Quel est le nom de cet identifiant ? En vous inspirant de la première partie, renommez cet attribut avec un nom pertinent (id_classe par exemple).
Réponse
Pour renommer cet attribut, nous suivons la syntaxe suivante.
names(occ_sol_dep)[names(occ_sol_dep) == 'ESA_WorldCover_2021_fr'] <- 'id_classe'En suivant une procédure vue précédemment, extrayez les polygones de forêt (classe 10) et stockez les dans une variable nommée foret_dep. Quelle syntaxe allez-vous utiliser ? Combien y a-t-il d’entités dans la couche résultat ?
Réponse
L’extraction des polygones de forêt va se faire selon une sélection attributaire sur le champ id_classe.
foret_dep <- occ_sol_dep[occ_sol_dep$id_classe == 10, ]
foret_dep class : SpatVector
geometry : polygons
dimensions : 1, 1 (geometries, attributes)
extent : 688457, 793143.2, 6617036, 6720873 (xmin, xmax, ymin, ymax)
coord. ref. : RGF93 v1 / Lambert-93 (EPSG:2154)
names : id_classe
type : <int>
values : 10Bien qu’il y ait de très nombreux patchs forestiers, la couche ne contient qu’une seule entité. Il s’agit en fait d’un polygone multi-parties.
Quelle est la superficie du département couvert par la forêt ? Quel pourcentage cela représente-t-il ?
Réponse
Pour répondre à cette question, il suffit d’ajouter un attribut de superficie et de le convertir en kilomètres carrés.
# ajout d'un attribut de superficie en km2
foret_dep$sup_km2 <- expanse(foret_dep) / 1000000
# on récupère la valeur dans une variable
sup_f <- foret_dep$sup_km2
# on calcule le pourcentage en comparant à la superficie totale du 58
pc_f <- sup_f / (expanse(dep_58) / 1000000) * 100Au total, sur notre département, la forêt recouvre 2876.84 km^2, ce qui représente 41.85 % du territoire départemental.
Nous avons maintenant une vision de la forêt à l’échelle d’un département via une procédure automatique. Il est possible de répéter la manipulation pour un autre département.
À vous de jouer
Quel est le pourcentage de prairies sur le département de l’Allier ?
Quelques géotraitements
Une des raisons d’être des SIG est de croiser des couches différentes afin de créer de nouvelles données répondant à une question. Les géotraitements permettent de croiser différentes couches entre elles afin de créer de nouvelles données. Nous allons ici nous intéresser, d’une part aux portions de communes sous coucvert forestier, afin de connaître la commune la plus boisée du département, et d’autre part nous allons extraire les portions de communes ne se trouvant pas sous forêt.
Quel géotraitement va-t-il vous permettre de quantifier la part de forêt dans chaque commune ?
Réponse
Le géotraitement permettant d’extraire l’espace occupé par deux entités à la fois (ici les communes et la forêt) est l’intersection.
La mise en application de ce géotraitement se fait de la façon qui suit. Au préalable il est nécessaire d’extraire les communes du département dans une nouvelle variable. Avant de faire l’intersection proprement dite, nous allons ajouter un attribut de superficie aux communes. Cela nous permettra par la suite de calculer une proportion de forêt par commune.
# si nécessaire, on importe la couche des communes
comm <- vect('./Data/communes_France.gpkg')
# extraction des communes du département
comm_dep <- comm[comm$CODE_DEPT == '58', ]
# ajout de la superficie communale en km2
comm_dep$sup_comm_km2 <- expanse(comm_dep) / 1000000
# intersection de ces communes et des forêts du département
comm_dep_foret <- intersect(comm_dep, foret_dep)
# calcul de la superficie de chaque patch forestier
comm_dep_foret$sup_km2 <- expanse(comm_dep_foret) / 1000000
# pourcentage de forêt dans chaque commune
comm_dep_foret$pc_foret <- (comm_dep_foret$sup_km2 / comm_dep_foret$sup_comm_km2) * 100
#comm_dep_foretNous pouvons ainsi identifier les communes les plus boisées et les moins boisées (en pourcentage).
# la commune au pourcentage de forêt max
comm_max <- comm_dep_foret[comm_dep_foret$pc_foret == max(comm_dep_foret$pc_foret)]$NOM_COM
pc_max <- comm_dep_foret[comm_dep_foret$pc_foret == max(comm_dep_foret$pc_foret)]$pc_foret
# la commune au pourcentage de forêt min
comm_min <- comm_dep_foret[comm_dep_foret$pc_foret == min(comm_dep_foret$pc_foret)]$NOM_COM
pc_min <- comm_dep_foret[comm_dep_foret$pc_foret == min(comm_dep_foret$pc_foret)]$pc_foretLa commune proportionnellement la plus boisée est LAVAULT-DE-FRETOY avec un pourcentage de forêt de 86.87 et la moins boisée est POUGNY avec un pourcentage de forêt de 3.24
À vous de jouer
Faites maintenant la manipulation inverse et extrayez les portions de communes qui ne sont pas boisées. Vous devrez utiliser le géotraitement différence.