Portail Forêt Ouverte

Le portail Forêt Ouverte a été décrit précédemment pour l’acquisition des données LiDAR. Ce portail contient plusieurs autres produits d’intérêt pour l’écologie végétale. L’un des produits les plus intéressants est la carte des peuplements écoforestiers.

Peuplements écoforestiers

La carte des peuplements écoforestiers est réalisée à partir d’interprétation de photographies aériennes ainsi que d’autres données pouvant aider à l’interprétation. Contrairement aux produits LiDAR, lesquels sont en format matriciel (i.e. images composées de pixels), les peuplements forestiers sont des données vectorielles (polygones ou points) et tabulaires (tables d’attributs).

La description des différentes composantes de la base de données est présentée dans ce document.

Accéder aux données

La façon la plus simple d’accéder aux données des peuplements écoforestiers est de charger la couche Catalogue > Données écoforestières (DIF) > Téléchargement des produits DIF dans Forêt Ouverte.

Comme nous avons découvert dans l’atelier sur les données LiDAR, la Station de biologie des Laurentides (SBL) se situe par malchance à la jonction de non pas deux, mais bien quatre feuillets DIF! C’est vraiment un malheureux hasard lorsqu’on considère que le feuillet 31H couvre à lui seul la grande région de Montréal, toute la Montérégie et une bonne partie de l’Estrie. Cela veut dire que pour obtenir les données couvrant l’étendue complète de la SBL, nous devons télécharger les données de ces quatre feuillets, bien que nous n’ayons en réalité besoin que d’une infime partie (un maigre 16 km\(^2\) pour être exact). Le bon côté est que cela nous forcera à apprendre de nouveaux outils pour fusionner des données vectorielles et tabulaires.

Le site FTP du MFFP

Sélectionner Accéder au répertoire de données Résultats d’inventaire et carte écoforestière nous mène au site FTP du MFFP.

Si une fenêtre de connexion apparaît, on peut se connecter en tant qu’invité.

Téléchargement des fichiers

Vous devriez alors voir deux fichiers .zip. Ce sont les mêmes données mais dans deux formats différents spécifiques à des versions différences de ArcMap (SIG propriétaire utilisé par le MFFP). L’un est pour ArcMap 9.3 (suffixe _93.zip)., l’autre pour ArcMap 10.0 (suffixe _10.zip). Nous avons besoin uniquement de celui terminant par _10.zip.

Il serait possible de télécharger les fichiers manuellement. Cependant, cela pourrait devenir fastidieux si on avait beaucoup de fichiers à télécharger. Nous allons donc utiliser une méthode de téléchargement directement via R. Notre code sera ainsi plus facile à reproduire. Pour ce faire nous définissons le lien FTP parent (lien ci-dessous) et utilisons la fonction download.file().

Avertissement: chacun des fichiers .zip fait au-dessus de 400 MB!

# Lien vers dossier FTP parent
lien <- 'ftp://transfert.mffp.gouv.qc.ca/Public/Diffusion/DonneeGratuite/Foret/DONNEES_FOR_ECO_SUD/Resultats_inventaire_et_carte_ecofor/'

# Téléchargher les fichiers .zip des quatre zones (31G, 31H, 31I, 31J)
download.file(paste0(lien, "31G/PRODUITS_IEQM_31G_10.zip"),
                     destfile = "donnees/PRODUITS_IEQM_31G_10.zip")
download.file(paste0(lien, "31H/PRODUITS_IEQM_31H_10.zip"),
                     destfile = "donnees/PRODUITS_IEQM_31H_10.zip")
download.file(paste0(lien, "31I/PRODUITS_IEQM_31I_10.zip"),
                     destfile = "donnees/PRODUITS_IEQM_31I_10.zip")
download.file(paste0(lien, "31J/PRODUITS_IEQM_31J_10.zip"),
                     destfile = "donnees/PRODUITS_IEQM_31J_10.zip")

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

Importation des données du service

Nous aurons besoin des limites administratives de la SBL afin d’extraire uniquement les données pertinentes à la station. Nous utiliserons la stratégie discutée précemment pour l’atelier sur l’acquisition des données LiDAR, soit importer directement le service hégergé sur ArcGIS Online.

Pour ce faire, 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. De plus, nous transformons les données (avec st_transform()) dans le système de référence spatiale (SRS) EPSG:2950. Il s’agit du datum NAD83(SCRS) avec projection MTM zone 8 (voir atelier LiDAR). Nous utilisons ce SRS car c’est une projection métrique (les unités sont en m) qui est la même que celle de nos données LiDAR.

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

Création de la zone tampon

Pour la zone d’extraction des données nous définissons une zone tampon de 250 m autour des limites de la SBL.

sbl_tampon <- st_buffer(sbl, 250)
plot(sbl, col = NA, border = 10, reset = FALSE, lwd = 3, main = NULL)
plot(sbl_tampon, col = NA, border = 1, add = TRUE, lwd = 3)

Transformation au format WKT

Par la suite, nous transformons cet zone tampon en format Well-known-text (WKT). Cette couche WKT sera par la suite utilisée comme filtre spatial lors de l’importation des données dans R, nous permettant ainsi de charger dans R uniquement les données d’intérêt. Cette stratégie évite de surcharger la mémoire vive de R pour rien.

Étant donné que les données de peuplements écoforestiers sont dans le SRS projeté EPSG:32198, soit la projection Conique conforme de Lambert (voir codes EPGS les plus utilisés au Québec), notre filtre WKT doit aussi être dans cette projection. C’est pour cela que nous utilisons st_transform(32198). Par la suite, nous laissons tomber les attributs de sbl_tampon en ne conservant que la géométrie (même s’il n’y avait qu’un seul attribut, soit FID), puis nous transformons le tout en texte avec st_as_text().

sbl_tampon_wkt <- sbl_tampon %>% 
  st_transform(32198) %>%
  st_geometry() %>% 
  st_as_text()

Importer les géodatabases ESRI

Nous importons les données pour le feuillet 31G à l’aide de read_sf(). Cette fonction reconnait l’extension .gdb (ESRI geodatabase), un format propriétaire d’ArcGIS. La fonction st_read() nous permet d’utiliser un filtre spatial en format WKT à l’aide de l’argument wkt_filter afin d’importer seulement les entités couvrant la SBL.

La géodatabase contient plusieurs classes d’entités (points, polygones) ainsi que des tableaux de données. Chacune des ces entités ou tableaux sont des couches. La couche PEE_ORI_31G est la couche centrale à importer en premier. PEE est l’abbréviation de PEuplements Écoforestiers, ORI est pour ORIginale (carte originale) C’est une couche de polygones avec plusieurs attributs importants. Chacun des polygone correspond à un peuplement écoforestier en particulier.

Nous transformons immédiatement les données de leur SRS d’origine, le EPSG:32198, au SRS EPSG:2950, pour que toutes nos données (incluant les données LiDAR) sont dans le même SRS.

pee_31G <- read_sf('donnees/PRODUITS_IEQM_31G_10/PRODUITS_IEQM_31G_10.gdb',
                   layer = 'PEE_ORI_31G',
                   wkt_filter = sbl_tampon_wkt) %>% 
  st_transform(2950)

Visualiser les données

Il est intéressant à ce point-ci de visualiser les données importées. Afin de pouvoir visualiser les polygones ainsi que les attributs associés, une carte interactive est utile. Une façon rapide de créer une carte interactive est à travers le package mapview. On peut visualiser une variable en particulier à l’aide de l’argument zcol. Dans ce cas, nous sélectionnons l’attribut TYPE_ECO (type écologique) qui classifie les différents peuplements.

require(mapview)
mapview(pee_31G, zcol = 'TYPE_ECO')

Les définitions des différents champs, incluant les types écologiques, sont détaillés dans ce document Excel. Par exemple, le type écologique FE32 est une Érablière à bouleau jaune sur dépôt minéral de mince à épais, de texture moyenne, de drainage mésique. Il s’agit du type écologique attendu selon le domaine bioclimatique de la SBL sur site mésique. Mais est-ce que cela correspond vraiment à la réalité sur le terrain? Est-ce que la SBL est couverte d’érablières à bouleau jaune, tel que la carte suppose? Nos inventaires de terrain nous permetterons de le vérifier.

Nous pouvons ignorer les autres champs pour l’instant.

Importer les données des autres feuillets

Comme la SBL couvre quatre feuillets, nous devons importer les données des trois autres feuillets (31H, 31I, 31J). La procédure est la même que celle utilisée pour le feuillet 31G donc nous montrons seulement le résultat, sans répéter le code R.

Nous utilisons la fonction générique plot() pour des cartographies statiques. On y ajoute les limites de la SBL et la zone tampon.

Fusionner les couches

Il faut maintenant fusionner les couches des différents feuillets en une seule. Comme le package sf tente de s’intégrer le mieux possible à l’univers tidyverse, plusieurs fonctions similaires conservent le même nom. Dans ce cas, on utilise la fonction générique rbind() pour row-bind étant donné que fusionner des couches sfconsiste à fusionner des lignes. En effet, selon le modèle de données chaque ligne correspond à un polygone dans un objets sf.

Si on observe bien, on peut constater dans les graphiques ci-dessus que certains polygones sont présents dans plus d’une couche. Pour éviter la duplication de ces polygones, on utilise distinct() avec l’argument .keep_all = TRUE pour conserver tous les attributs des polygones (et non seulement leurs géométries).

Finalement, nous utilisons mapview() pour visualiser de façon interactive le résultat.

pee_sbl <- rbind(pee_31G,
                pee_31H,
                pee_31I,
                pee_31J) %>% 
  distinct(.keep_all = TRUE)
mapview(pee_sbl, zcol = 'TYPE_ECO')

Importer les données tabulaires

Les géodatabases des peuplements écoforestiers contiennent plusieurs autres entités et tables, dont deux sont d’intérêt particuliers pour le cours. Les détails de toutes les tables et champs sont disponibles ici.

Les deux tables d’intérêt pour le cours sont les suivantes:

  1. ETAGE_ORI_PROV: Variables dendrométriques photointerprétées détaillées par étage des peuplements écoforestiers originaux
  2. ESSENCE_ORI_PROV: Composition en essences détaillée des peuplements écoforestiers originaux

Ces variables ne sont pas disponibles pour tous les peuplements mais seulement un sous-ensemble de ces derniers.

Variables dendrométriques par étage

Nous importons les données par étage pour chaque feuillet à l’aide de read_sf(). Nous sélectionnons seulement les données de la SBL en filtrant à l’aide de semi_join(). Cette fonction permet de sélectionner les entrées dans le premier objet sf(dans ce cas-ci, pee_31G_etage, qui contient les données de tout le feuillet 31G au complet) qui sont présents dans le deuxième objet (ici, pee_sbl). La variable permettant de joindre les deux jeux de données ensemble est GEOCODE. Comme elle est la seule variable présente dans les deux objets il n’est pas nécéssaire de la spéficier.

# 31G
pee_31G_etage <- read_sf('donnees/PRODUITS_IEQM_31G_10/PRODUITS_IEQM_31G_10.gdb',
                   layer = 'ETAGE_ORI_31G') %>% 
  semi_join(pee_sbl)
## Warning: no simple feature geometries present: returning a data.frame or
## tbl_df
## Joining, by = "GEOCODE"

La fonction nous avertit qu’aucune géométrie n’est présente dans le jeu de données. C’est normal car il ne s’agit que d’un tableau de données.

Nous faisons la même chose pour les trois autres feuillets (code R omis). Cela génère trois autres objets ayant la même structure que pee_31G_etage:

  • pee_31H_etage
  • pee_31I_etage
  • pee_31J_etage

Fusionner les variables par étage

La fusion des données des variables par étage se fait encore une fois avec rbind(). Nous utilisons distinct() afin de supprimer les entrées qui sont présentes dans plus d’un jeu de données et seraient donc dupliquées.

pee_etage_sbl <- rbind(pee_31G_etage,
                pee_31H_etage,
                pee_31I_etage,
                pee_31J_etage) %>% 
  distinct()
pee_etage_sbl
## # A tibble: 67 x 7
##    GEOCODE       ETAGE TY_COUV_ET DENSITE HAUTEUR CL_AGE_ET ETA_ESS_PC     
##    <chr>         <chr> <chr>        <int>   <int> <chr>     <chr>          
##  1 -421964,37+2… SUP   F               95      11 30        EO30BP20ES20FN…
##  2 -422020,29+2… SUP   F               95      20 VIN       ES50EO20HG20BP…
##  3 -422141,63+2… SUP   F               95      11 30        EO30BP20ES20PE…
##  4 -422197,79+2… SUP   M               75      12 30        BP20PE20EO20FN…
##  5 -422317,43+2… SUP   F               85      18 70        BP70EO10BJ10RX…
##  6 -422346,15+2… SUP   M               65      16 JIN       BJ20EO20BP10ES…
##  7 -422365,46+2… SUP   F               65      17 70        BP30EO20ES10BJ…
##  8 -422464,74+2… SUP   M               75      20 70        BP50EO10BJ10SB…
##  9 -422494,32+2… SUP   F               85      19 JIN       ES50EO20BP10BJ…
## 10 -422513,35+2… SUP   F               85      19 JIN       ES50EO20BP10BJ…
## # … with 57 more rows

Le résultat indique que seulement 67 entités sur les 322 présentes dans pee_sbl ont fait l’objet d’une photo-interprétation détaillée des variables dendrométriques.

Visualiser les peuplements avec données par étage

Afin de visualiser où se situent les 67 entités pour lesquelles des variables dendrométriques détaillées sont disponibles, il est possible de filtrer la couche pee_sbl à l’aide de la couche pee_etage_sbl et de fusionner leurs attributs à l’aide de right_join().

pee_etage_sbl_geom <- pee_sbl %>% 
  right_join(pee_etage_sbl)
## Joining, by = "GEOCODE"
mapview(pee_etage_sbl_geom, zcol = 'TYPE_ECO')

Comme nous pouvons les constater, les entités autour du Lac Croche n’ont pas de mesures détaillées de la composition en espèces et autres variables dendrométriques.

Composition en essences détaillée

L’importation des données de composition en essences détaillées est essentiellement réalisée de la même façon qu’avec les variables dendrométriques. Nous montrons ci-dessous le code R seulement pour la zone 31G.

# 31G
pee_31G_essence <- read_sf('donnees/PRODUITS_IEQM_31G_10/PRODUITS_IEQM_31G_10.gdb',
                   layer = 'ESSENCE_ORI_31G') %>% 
  semi_join(pee_sbl)
## Warning: no simple feature geometries present: returning a data.frame or
## tbl_df
## Joining, by = "GEOCODE"

Fusionner les données d’essences

La même méthode décrite pour les variables dendrométriques s’applique pour fusionner les données des quatre feuillets.

pee_essence_sbl <- rbind(pee_31G_essence,
                pee_31H_essence,
                pee_31I_essence,
                pee_31J_essence) %>% 
  distinct()
pee_essence_sbl
## # A tibble: 302 x 4
##    GEOCODE              ETAGE ESSENCE ST_ESS_PC
##    <chr>                <chr> <chr>       <int>
##  1 -421964,37+237575,08 SUP   EO             30
##  2 -421964,37+237575,08 SUP   FN             10
##  3 -421964,37+237575,08 SUP   PE             10
##  4 -421964,37+237575,08 SUP   ES             20
##  5 -421964,37+237575,08 SUP   BP             20
##  6 -421964,37+237575,08 SUP   SB             10
##  7 -422020,29+238516,08 SUP   HG             20
##  8 -422020,29+238516,08 SUP   BP             10
##  9 -422020,29+238516,08 SUP   EO             20
## 10 -422020,29+238516,08 SUP   ES             50
## # … with 292 more rows

Sauvegarder les données

Nous allons sauvegarder les données dans un GeoPackage. Contrairement aux format Geodatabase (.gdb) lequel est un format propriétaire défini par ESRI (propriétaire d’ArcGIS), le GeoPackage est un format libre défini par le Open Geospatial Consortium. Tout comme une géodatabase ESRI, un GeoPackage peut contenir plusieurs couches en un seul fichier.

Pour ce faire nous utilisons st_write() et sauvegardons la couche des peuplements écoforestiers (les 322 polygones).

# dir.create('resultats') # Pour créer un dossier pour stocker les résultats
st_write(pee_sbl,
         dsn = 'resultats/pee_sbl.gpkg')

Malheureusement st_write() ne permet pas (du moins pour l’instant) d’ajouter des couches sans entités géométriques (i.e. des tableaux de données) dans un .gpkg. Mais puisque qu’un GeoPackage est essentiellement une base de données SQLite, il est possible d’inclure les tables pee_etage_sbl et pee_essence_sbl à l’aide d’outils R pour ces bases de données.

require(RSQLite)

# Se connecter au GeoPackage
con <- dbConnect(SQLite(), dbname = 'resultats/pee_sbl.gpkg')

# Ajouter les tables dans la base de données SQLite (le .gpkg)
dbWriteTable(con, "pee_etage_sbl", pee_etage_sbl)
dbWriteTable(con, "pee_essence_sbl", pee_essence_sbl)

Et voilà! Nous sommes maintenant prêts à faire des requêtes et analyses sur les peuplements écoforestiers de la SBL.