Portail Forêt Ouverte

Accéder au portail

Premièrement, il faut accéder au portail Forêt Ouverte. Ce portail web donne accès à un large éventail de données ouvertes géospatiales en lien avec la forêt québécoise, notamment plusieurs produits LiDAR. Il est maintenu par le Ministère des Forêts, de la Faune et des Parcs (MFPP).

Afficher les couches LiDAR

Dans le menu en haut à gauche, accéder au Cartes prédéfinies et choisir la catégorie LiDAR. Afficher la couche Téléchargement - LiDAR (pleine résolution).

Des informations sur les différents produits LiDAR du MFFP sont disponibles ici. Un guide d’utilisation plus détaillé est aussi disponible ici.

Centrer sur la région d’intérêt

Maintenant, centrons la carte sur la zone d’intérêt, dans notre cas le lac Croche de la Station de biologie des Laurentides (SBL). En cliquant sur la carte, cela fera apparaître un rectangle bleu correspondant à un feuillet 1/20,000 en particulier pour accéder aux données LiDAR. Il est à noter que la SBL, malheureusement, se trouve à l’intersection de quatre feuillets différents. Cela veut dire que pour obtenir les données LiDAR couvrant l’étendue complète de la SBL, il faudra télécharger les données de quatre feuillets (31G16NE, 31H13NO, 31I04SO et 31J01SE). Cela complique beaucoup l’acquisition des données, mais nous permettra d’apprendre des fonctions utiles, notamment comment réaliser une mosaïque d’images séparées.

Avertissement: Les données LiDAR sont généralement très volumineuses car la résolution spatiale des produits est très fine (e.g. 1-2 m).

Il est aussi utile de connaître l’année d’acquisition des données LiDAR sur le terrain. Pour ce faire, il faut afficher la couche Années d’acquisition du LiDAR sur la carte. Cela permet de constater que les données LiDAR que nous allons télécharger pour la SBL datent de 2014 (partie sud-est) et 2015 (reste de la SBL). Bien que le modèle numérique de terrain (MNT), lequel représente la surface du sol, est peu variable dans le temps, le modèle de hauteur de canopée (MHC) représente la hauteur des arbres dominants, laquelle varie selon la croissance des arbres.

Accéder au site FTP

Nous pouvons accéder au site FTP en cliquant sur un lien pointant vers l’un des produits LiDAR. Cela nous guidera vers la source des données LiDAR pour le feuillet selectionné sur la carte web.

Il est évidemment possible de télécharger manuellement les différents produits LiDAR pour chacun des feuillets 1/20,000 séparément à partir du site FTP dans le fureteur web. Cependant, il est fastidieux de procéder de la sorte, car dans notre cas il y a quatre feuillets (i.e. zones) différents, contenant chacun plusieurs produits LiDAR d’intérêt: modèles de hauteurs de canopée (MHC), modèles numériques de terrain (MNT) et pentes.

Nous procéderons donc par un téléchargement semi-automatique des différents produits LiDAR d’intérêt pour les quatre feuillets couvrant la SBL à travers R. Cela présente deux avantages principaux: (1) ça assure la reproductibilité de notre analyse (car tout le monde a accès aux données directement à travers le script) (2) facilite le téléchargement de fichiers lorsque ceux-ci sont nombreux.

Téléchargement des données avec R

Pour faciliter le téléchargement à travers R, nous avons besoin des liens pointant vers les fichiers. Premièrement, nous définissons le dossier FTP parent. Il s’agit simplement de copier-coller le lien URL du site FTP.

lien <- 'ftp://transfert.mffp.gouv.qc.ca/Public/Diffusion/DonneeGratuite/Foret/IMAGERIE/Produits_derives_LiDAR/31G/31G16NE/'

Deuxièmement, nous devons générer une liste tous les fichiers présents dans ce dossier. Pour ce faire nous utilisons la fonction getURL du package RCurl. Le code est ici un peu rébarbatif mais l’objectif est simplement d’obtenir un vecteur contenant le nom de tous les fichiers. Nous aurons besoin de ces noms de fichiers pour le téléchargement.

require(RCurl)
fichiers <- getURL(lien, ftp.use.epsv = TRUE, dirlistonly = TRUE)
fichiers <- unlist(strsplit(fichiers[[1]], "\n"))
fichiers
## [1] "MHC_31G16NE.lyr"       "MHC_31G16NE.tif"       "MNT_31G16NE.lyr"      
## [4] "MNT_31G16NE.tif"       "MNT_Ombre_31G16NE.lyr" "MNT_Ombre_31G16NE.tif"
## [7] "Pentes_31G16NE.lyr"    "Pentes_31G16NE.tif"    "Thumbs.db"

Troisièment, nous créons un dossier local (sur notre ordinateur) pour entreposer les données et puis nous procédons au téléchargement de tous les fichiers. Pour ce faire, nous utilisons une boucle for avec la fonction download.file(). Nous sauvegardons le nom du dossier dans un objet mon_dossier pour légèrement simplifier le code dans les prochaines étapes.

Avertissement: les fichiers .tif peuvent être très volumineux (plusieurs centaines de MB!)

mon_dossier <- 'donnees/'
# dir.create(mon_dossier) # Si vous voulez le créer à travers le script
for (i in 1:length(fichiers)) {
  download.file(paste0(lien, fichiers[i]),
         destfile = paste0(mon_dossier, fichiers[i]))
  }

Afin de télécharger toutes les données LiDAR couvrant le territoire de la SBL, il faut ensuite répéter ces étapes pour les feuillets 31G16NE, 31H13NO, 31I04SO et 31J01SE. Le code n’est pas répété ici pour éviter la redondance, mais les fichiers sont téléchargés vers le même dossier local sur notre ordinateur.

Inspecter les produits LiDAR

Une fois téléchargé, c’est une bonne idée d’inspecter les couches. Pour ce faire, un système d’information géographique (SIG), tel que le logiciel libre QGIS peut être utilisé. Ces logiciels sont optimisés pour la visualisation et manipulation des données géospatiales. Cependant, R est maintenant un SIG à part entière grâce à des nombreux packages spécialisés. Nous allons donc continuer avec R, car il est plus agréable de travailler dans un seul et même système.

Le nouveau package stars est un nouveau package conçu spécifiquement pour les données matricielles (raster) telles que produits LiDAR. Le package raster, plus ancien mais bien établi, pourrait aussi être utilisé mais à mon avis stars possède une syntaxe plus intuitive et interagit mieux avec les packages de l’univers tidyverse. Le package stars dépend du package sf, lequel traite les données vectorielles.

Un autre avantage du package stars est la capacité de travailler avec les données d’imagerie sans avoir à charger toute l’image au complet en mémoire vive (RAM). Cela est important pour R car les objets temporaires dans une session R sont stockées dans le RAM, lequel est généralement limité. Le package a créé un nouveau type d’objet spatial, les objets stars_proxy, pour répondre à ce besoin. Ces objets sont beaucoup plus légérs à manipuler et à visualiser dans R.

Pour charger une image sous la forme d’un objet stars_proxy, il faut utiliser la fonction générale read_stars(), mais avec l’argument proxy = TRUE. Pour illustrer ceci, chargeons en premier le MHC pour le feuillet 31G16NE. Nous utilisons le préfixe mon_dossier pour avoir le lien complet pointant vers le fichier.

require(stars)
MHC_31G16NE <- read_stars(paste0(mon_dossier, 'MHC_31G16NE.tif'),
                          proxy = TRUE)
MHC_31G16NE
## stars_proxy object with 1 attribute in file:
## $MHC_31G16NE.tif
## [1] "donnees/MHC_31G16NE.tif"
## 
## dimension(s):
##   from    to  offset delta                       refsys point values    
## x    1 21000  246000     1 PROJCS["unnamed",GEOGCS["... FALSE   NULL [x]
## y    1 15000 5096000    -1 PROJCS["unnamed",GEOGCS["... FALSE   NULL [y]

Taper le nom MHC_31G16NE de l’objet dans la console R donne quelques informations générales. En première ligne, on voit que l’objet est bien de classe stars_proxy et qu’il contient un seul attribut ou variable. C’est normal car ce produit ne doit que contenir uniquement la hauteur de la canopée (en mètres) pour chacun des pixels de l’image. Par la suite, le nom et la source du fichier sont donnés.

Le système de référence spatiale (SRS)

La section dimension(s) contient des informations importante de type géographique sur l’objet. On voit que l’image est une image 2D (i.e. elle contient 2 dimensions x et y). L’information la plus importante concerne le système de référence spatiale (SRS), dans la colonne refsys. Nous savons selon la section Informations complémentaires des métadonnées des produits LiDAR du MFFP que le SRS devrait être le NAD83(SCRS) comme projection cartographique la Mercator Transverse Modifiée (MTM).

La projection MTM comporte plusieurs zones, chacune correspondant à une bande de 3 degrés de longitude. Il est possible de télécharger des fichier .kml (qui peuvent être directement ouverts dans Google Earth) contenant les différentes zones MTM pour le Canada ici. On peut donc voir que le SBL correspond à la zone MTM 8 selon l’image ci-dessous.

Les SRS et projections les plus fréquemment utilisés ont un code EPSG afin de faciliter l’échange d’informations de SRS. Les codes EPSG des SRS et projections les plus fréquemment utilisés au Québec sont listés dans ce fichier PDF. Le code EPSG correspondant aux produits LiDAR pour la zone de la SBL est donc le EPSG:2950, puisque le SRS est le NAD83(SCRS) et la projection MTM zone 8. Il est possible de trouver plus d’informations sur des codes EPSG particuliers sur le site Spatial Reference. La base de données officielle de tous les codes EPSG, plus complexe à naviguer mais complète, est disponible ici.

Pour obtenir plus d’informations sur le SRS de notre objet MHC_31G16NE, on peut utiliser la fonction st_crs() (CRS = Coordinate Reference System, un synonyme de SRS) du package sf. Il est à noter que la majorité des fonctions de ce package débutent par st_, une abbréviation de space-time.

st_crs(MHC_31G16NE)
## Coordinate Reference System:
##   User input: PROJCS["unnamed",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-73.5],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",304800],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]] 
##   wkt:
## PROJCS["unnamed",
##     GEOGCS["GRS 1980(IUGG, 1980)",
##         DATUM["unknown",
##             SPHEROID["GRS80",6378137,298.257222101],
##             TOWGS84[0,0,0,0,0,0,0]],
##         PRIMEM["Greenwich",0],
##         UNIT["degree",0.0174532925199433]],
##     PROJECTION["Transverse_Mercator"],
##     PARAMETER["latitude_of_origin",0],
##     PARAMETER["central_meridian",-73.5],
##     PARAMETER["scale_factor",0.9999],
##     PARAMETER["false_easting",304800],
##     PARAMETER["false_northing",0],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]]]

Tous les détails du SRS et de la projection sont listés dans le format WKT. Ces informations ont été prélevées par la fonction read_stars() lors de l’importation car le fichier .tif est en fait un GeoTIFF, soit une image géoréférencée. Cependant, le code EPSG:2950 n’a pas été reconnu par stars et le nom de la projection dans l’objet est inconnue (unnamed). Toutefois, les paramètres du SRS et de la projection sont exacts. Cependant, pour éviter toute confusion, il est préférable d’associer explicitement le code EPSG:2950 au SRS de l’objet MHC_31G16NE; il faudra le faire pour toutes les images par la suite. On le fait à l’aide de st_set_crs().

MHC_31G16NE <- st_set_crs(MHC_31G16NE, 2950)
st_crs(MHC_31G16NE)
## Coordinate Reference System:
##   User input: EPSG:2950 
##   wkt:
## PROJCS["NAD83(CSRS) / MTM zone 8",
##     GEOGCS["NAD83(CSRS)",
##         DATUM["NAD83_Canadian_Spatial_Reference_System",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6140"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4617"]],
##     PROJECTION["Transverse_Mercator"],
##     PARAMETER["latitude_of_origin",0],
##     PARAMETER["central_meridian",-73.5],
##     PARAMETER["scale_factor",0.9999],
##     PARAMETER["false_easting",304800],
##     PARAMETER["false_northing",0],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]],
##     AXIS["E(X)",EAST],
##     AXIS["N(Y)",NORTH],
##     AUTHORITY["EPSG","2950"]]

Il y a une grande quantité d’information listée sur le SRS du code EPSG:2950. Il n’est pas important de comprendre tous les paramètres, mais on voit notamment que le bon nom de code EPSG est maintenant listé. On voit aussi que le méridien central de la zone MTM 8 correspond à la longitude -73.5\(^\circ\) (négatif car à l’ouest du méridien d’origine), ce qui est exact. De plus, on peut aussi voir que les unités des dimensions XY de la carte sont en m.

On comprend maintenant mieux les autres informations stockées dans dimension(s), qu’on peut obtenir directement avec st_dimensions(). Les détails derrière la structure des objets stars peuvent être consultés ici.

st_dimensions(MHC_31G16NE)
##   from    to  offset delta    refsys point values    
## x    1 21000  246000     1 EPSG:2950 FALSE   NULL [x]
## y    1 15000 5096000    -1 EPSG:2950 FALSE   NULL [y]

En gros, puisqu’il s’agit d’une projection MTM, la dimension x correspond au Easting (axe Est-Ouest) tandis que la dimension y correspond au Northing (axe Nord-Sud). Les colonnes from et to indiquent que l’image est 21 000 pixels de large et 15 000 pixels de haut. Les offset indiquent la position du pixel d’origine de l’image, laquelle correspond pour un GeoTIFF au coin supérieur gauche de l’image (nord-ouest). Les valeurs delta indiquent la résolution spatiale de chaque pixel en x et en y. Ici, chaque pixel correspond à un carré de 1 m par 1 m (parce que les unités de la projection sont des m). La raison pour laquelle la valeur delta et y est négative (-1) est parce que les indices de pixels vont du nord vers le sud. Pour l’instant, nous pouvons ignorer les colonnes point et values.

Visualisation de l’image

Après ce déluge d’informations, nous voulons maintenant voir l’image! La fonction générique plot() fonctionne pour les objets stars et stars_proxy peut être utilisée car c’est la plus simple. Il existe plusieurs autres options de cartographie plus avancées.

plot(MHC_31G16NE)

Par défaut, une image avec une seule bande est reproduite avec des nuances de gris. On voit que les hauteurs de canopée vont d’environ zéro (e.g. pour les cours d’eau et lacs) jusqu’à environ 25 m (pour les arbres les plus haut). Il est à noter que les MHC par LiDAR ont tendance à légèrement sous-estimer la hauteur réelle des arbres, comme il est expliqué dans ce document.

Bref, c’est réaliste, mais ce n’est pas très joli. Pour améliorer l’aspect visuel, nous pourrions utiliser la fonction topo.colors() et sélectionner une étendue de 20 couleurs allant de bleu (faible hauteur) à rouge (canopée très haute), en passant par le vert (hauteur intermédiaire). Le résultat ci-dessous est beaucoup plus intéressant.

plot(MHC_31G16NE, col = topo.colors(20))

Définir la zone d’intérêt: la SBL

Le fichier du MHC du feuillet 31G16NE présenté ci-dessus fait à lui seul 1.27 GB. C’est très lourd! La raison pour laquelle R arrive à facilement et rapidement à cartographier l’objet sans surchager notre RAM tient au fait que nous avons chargé l’image en tant qu’objet stars_proxy. En effet, pour ces objets le package stars charge en mémore seulement les métadonnées de l’image, pas les données de l’image elle-même. Lorsque nous demandons à R de faire une carte avec plot() à partir d’un objet stars_proxy, R réduit la résolution de l’image avant de charger l’image elle-même en mémoire de façon à avoir une résolution qui correspond à celle de notre moniteur (écran). En effet, il ne sert à rien de charger une image ayant une résolution beaucoup plus élevée que celle de notre moniteur.

Nous savons aussi que la SBL ne couvre qu’une toute petite partie de cette image (la partie nord-est de l’image, où on peut voir le Lac Croche), ainsi que des images des trois autres feuillets. Bref, une grande partie du 1.27 GB de cette image nous n’est d’aucune utilité. Pour alléger le fichier en ne conservant que les données pertinentes à la SBL, il serait judicieux d’extraire seulement la zone d’intérêt de chaque feuillet. Par la suite, nous pourrons fusionner les images extraites en une seule mosaïque.

Nous devons connaître les limites administratives de la SBL afin de procéder à l’extraction. C’est la prochaine étape.

Importation des limites de la SBL

Un service hebergé des limites administratives de la SBL est disponible dans ArcGIS Online. L’avantage de télécharger les données de cette source centralisée est au niveau de la reproducibilité: tout le monde peut accéder directement aux données via une connection internet.

Le code EPSG du service est le 3857, une projection souvent appellée Web Mercator. Cette projection est une modification de la projection Mercator standard (UTM). Popularisée par Google Maps, elle est devenue le standard du web. Cependant, il faut noter que le SRS géographique de base du système Web Mercator est généralement le EPSG:4326. C’est ce code qui en fait est utilisé lors de l’importation. Ces deux SRS sont intimement liés: en effet, le datum WGS84 est celui utilisé par le système GPS et aussi par Google Maps et autres. C’est le SRS géographique (i.e. latitude et longitude) de loin le plus répandu sur le web.

Bien qu’il soit possible d’exporter le service en format .shp à partir d’ArcGIS Online, il est plus simple d’importer directement le service dans R. Pour ce faire, nous utilisons le package esri2sf. Comme ce package n’est pas disponible sur CRAN, il faut l’installer directement de Github.

require(devtools)
install_github("yonghah/esri2sf")

Nous utilisons la fonction esri2sf() du package du même nom pour l’importation et la transformation à un objet de classe sf. Nous devons copier-coller le lien du service. Ce lien se trouve au page de la page web du service.

Avertissement: il faut ajouter /0 à la fin du lien pour que l’importation fonctionne.

Nous importons seulement le premier attribut (FID) dans la commande ci-dessous avec l’argument outFields.

require(esri2sf)
service <- 'https://services5.arcgis.com/ncjvTc8Zi7Btecnp/arcgis/rest/services/Limites_SBL/FeatureServer/0' # Ajouter /0 à la fin!!!
sbl <- esri2sf(service, outFields = 'FID')
## [1] "Feature Layer"
## [1] "esriGeometryPolygon"

Nous devons maintenant transformer l’objet sbl en EPSG:2950 afin qu’il soit compatible avec l’imagerie, car le EPSG:4326 est un SRS géographique alors que le EPSG:2950 est un SRS projeté. Notez aussi que le EPSG:2950 utilise le datum NAD93(SCRS), alors que le datum du EPSG:4326 est le WGS84.

La fonction st_transform() nous permet d’effectuer cette transformation de datum et de projection.

sbl_2950 <- st_transform(sbl, 2950)

Rajoutons maintenant les limites de la SBL en rouge sur la carte du MHC. Cela nous permet de vérifier que la SBL ne couvre qu’une toute petite partie (au nord-est) de l’image. Pour ajouter une deuxième couche sur une carte avec plot(), il faut utiliser reset = FALSE dans la première commande et add = TRUE dans la deuxième. Les autres arguments sont purement esthétiques.

plot(MHC_31G16NE, col = topo.colors(20), reset = FALSE)
plot(sbl_2950, col = NA, border = 10, add = TRUE, lwd = 3)

Définition d’une zone tampon

Afin d’extraire les données LiDAR, nous définissons une zone tampon de 250 m autour des limites de la SBL. Cela permet de sauvegarder les données en bordure de la SBL, qui pourraient être d’intérêt. Cela évite aussi de devoir refaire l’analyse si jamais les limites de la SBL n’étaient pas tout à fait précises.

Une zone tampon se définit à l’aide de st_buffer(). La largeur de la zone tampon est dans les unités de la projection, ici en m. Dans le graphique ci-dessous les limites de la SBL sont en rouge et la zone tampon en noir.

sbl_tampon <- st_buffer(sbl_2950, 250) # en m, car ce sont les unités du EPSG:2950
plot(sbl_2950, col = NA, border = 10, reset = FALSE, lwd = 3, main = NULL)
plot(sbl_tampon, col = NA, border = 1, add = TRUE, lwd = 3)

Extraction des données

L’approche choisie est d’extraire les données de la SBL pour chaque feuillet séparément avant de joindre ces sous-régions en une seule mosaïque. Cette approche est mois coûteuse en temps de calcul que le contraire. La fonction st_crop() nous permet d’extraire les données un utilisant une couche vectorielle (ici, sbl_tampon) comme bordure.

MHC_31G16NE_sbl <- st_crop(MHC_31G16NE, sbl_tampon)
plot(MHC_31G16NE_sbl, col = topo.colors(20), reset = FALSE)
plot(sbl_2950, col = NA, border = 10, add = TRUE, lwd = 3)

Il faut ensuite répéter ces mêmes opérations pour les trois autres feuillets. Cette fois, nous utilisons l’opérateur %>% (qu’on appelle un pipe). Cet opérateur indique une séquence d’opérations devant être enchaînées les unes après les autres. Il se lit comme ensuite. Cela simplifie le code et évite d’avoir à créer des objets temporaires pour sauvegarder les résultats intermédiaires. Le résultat de l’opération d’avant devient le premier argument de la fonction à l’opération suivante. La possibilité d’utiliser %>% est l’un des avantages des packages sf et stars, par rapport à sp et raster, respectivement.

MHC_31H13NO_sbl <- read_stars(paste0(mon_dossier, 'MHC_31H13NO.tif'),
                              proxy = TRUE) %>% 
   st_set_crs(2950) %>% 
   st_crop(sbl_tampon)

MHC_31I04SO_sbl <- read_stars(paste0(mon_dossier, 'MHC_31I04SO.tif'),
                              proxy = TRUE) %>% 
   st_set_crs(2950) %>% 
   st_crop(sbl_tampon)

MHC_31J01SE_sbl <- read_stars(paste0(mon_dossier, 'MHC_31J01SE.tif'),
                              proxy = TRUE) %>% 
   st_set_crs(2950) %>% 
   st_crop(sbl_tampon)

Regardons les zones extraites des trois autres feuillets.

Création d’une mosaïque de hauteur de canopée

Nous pouvons maintenant fusionner les différents jeux de données provenant des quatre feuillets en une seule mosaïque à l’aide de st_mosaic(). Il suffit de lister les différents objets stars_proxy un après l’autre. La création de mosaïques peut être une opération assez lourde pour l’ordinateur (quelques minutes).

Avertissement: au moment d’écrire ces lignes une limitation du package stars fait en sorte qu’il est absolument nécessaire d’ajouter l’option options = c("-srcnodata", "nan") afin d’obtenir une mosaïque complète, sans données manquantes. Le problème est documenté ici. Il est possible que ce comportement soit modifié dans une version future.

MHC_sbl <- st_mosaic(MHC_31G16NE_sbl,
                     MHC_31H13NO_sbl,
                     MHC_31I04SO_sbl,
                     MHC_31J01SE_sbl,
                     options = c("-srcnodata", "nan"))

Il est important de sauvegarder le résultat final dans un fichier, ici un GeoTIFF. Pour ce faire, il faut en premier transformer l’objet stars_proxy (lequel, souvenez-vous, ne contient que les métadonnées) en objet stars (lequel contient les données) avec st_as_stars(). Pour créer le GeoTiff, on utilise write_stars() avec l’extension .tif.

MHC_sbl %>% 
   st_as_stars() %>% 
   write_stars('resultats/MHC_sbl.tif')

Finalement, on peut visualiser le résultat. Joli!

Bravo! Cette mosaïque peut être maintenant utilisée pour des analyses futures, par exemple sur la distribution de la hauteur de la canopée à la SBL en fonction de la topographie.

Mosaïques des autres produits LiDAR

Les autres produits LiDAR d’intérêts du portail Forêt Ouverte sont les:

La procédure pour extraire les données pertinentes pour la SBL et les fusionner en mosaïque est exactement la même que celle décrite précédemment pour le MHC. Nous ne répéterons donc pas le code, mais présentons seulement les résultats finaux. Des fichiers GeoTiff sont sauvegardés pour les analyses futures.

```