Using rstac and gdalcubes packages

IO is a catalog of environmental raster data layers that can be used in biodiversity-related modeling contexts. This resource is in the form of a STAC catalog, which is accessible through various methods, including R packages rstac and gdalcubes.

The available layers in the catalog can be visualized here.

The data layers in IO are in COG (Cloud Optimized Geotiff) format, which allows optimal access with remote queries. For example, it is possible to extract only a small region from a global file and transform its resolution and reference coordinate system without ever having to download it entirely.

Here are some examples of queries to interact with the IO STAC catalog and COG files.

To start, install the following packages. Once the packages are installed, it is not necessary to re-install them at each R session.

library(gdalcubes)
library(rstac)
library(terra)

Connect to the STAC catalog

s_obj <- stac("https://io.biodiversite-quebec.ca/stac/")

List collections

collections <- s_obj |> collections() |> get_request()

See collections and their descriptions

library(knitr)
df<-data.frame(id=character(),title=character(),description=character())
for (c in collections[['collections']]){
  df<-rbind(df,data.frame(id=c$id,title=c$title,description=c$description))
}
kable(df)
id title description
chelsa-clim CHELSA Climatologies CHELSA Climatologies
gfw-lossyear Global Forest Watch - Loss year Global Forest Watch - Loss year
ghmts Global Human Modification of Terrestrial Systems The Global Human Modification of Terrestrial Systems data set provides a cumulative measure of the human modification of terrestrial lands across the globe at a 1-km resolution. It is a continuous 0-1 metric that reflects the proportion of a landscape modified, based on modeling the physical extents of 13 anthropogenic stressors and their estimated impacts using spatially-explicit global data sets with a median year of 2016.
gfw-treecover2000 Global Forest Watch - Tree cover 2000 Global Forest Watch - Tree cover 2000
gfw-gain Global Forest Watch - Gain Global Forest Watch - Gain
chelsa-monthly CHELSA monthly timeseries CHELSA monthly timeseries
esacci-lc ESA Land cover time series Global annual land cover map at 300 m based on the full archives of 300 m MERIS, 1 km SPOT-VEGETATION, 1 km PROBA-V and 1 km AVHRR surface reflectance 7-day composites
chelsa-clim-proj CHELSA Climatologies Projections CHELSA Climatologies Projections
stressors_qc Potential Stressors of Québec’s Biodiversity Indicators calculated across Québec that capture the various stressors of the biodiversity in the province.
silvis Silvis Dynamic Habitat Indices The DHIs are designed for biodiversity assessments and to describe habitats of different species.
earthenv_topography EarthEnv - Topography A fully standardized and global multivariate product of different terrain features to support many large-scale research applications. The product is based on the digital elevation model products of global 250 m GMTED2010. From https://www.earthenv.org/topography
earthenv_landcover EarthEnv - Consensus Land Cover - Full version Consensus land cover dataset containing 12 data layers, each of which provides consensus information on the prevalence of one land-cover class. All data layers contain unsigned 8-bit values and the valid values range from 0-100, representing the consensus prevalence in percentage. All data layers have a spatial extent from 90ºN - 56ºS and from 180ºW - 180ºE, and have a spatial resolution of 30 arc-second per pixel (~1 km per pixel at the equator).
earthenv_habitat_heterogeneity EarthEnv - Habitat heterogeneity The datasets contain 14 metrics quantifying spatial heterogeneity of global habitat at multiple resolutions based on the textural features of Enhanced Vegetation Index (EVI) imagery acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS).
accessibility_to_cities Accessibility - Time to access cities This global accessibility map enumerates land-based travel time to the nearest densely-populated area for all areas between 85 degrees north and 60 degrees south for a nominal year 2015. Densely-populated areas are defined as contiguous areas with 1,500 or more inhabitants per square kilometer or a majority of built-up land cover types coincident with a population centre of at least 50,000 inhabitants.
qc_pilot_env Environmental layers for the Quebec pilot BON optimization project. Environmental layers for the Quebec pilot BON optimization project.
global-mammals Global habitat availability for mammals from 2015-2055 Habitat availability for each mammal species (total 5090 species) in 0.25 degrees cells, from 2015 to 2055 in 5 year intervals, for the SSP/RCP scenarios. Developed for the BES-SIM project by the University of Sapienza team using the InSiGHTS (Integrated Scenarios of Global Habitat for Terrestrial Species) model
fragmentation-rmf Relative Magnitude of Fragmentation (RMF) The Relative Magnitude of Fragmentation (RMF) measures the fragmentation of specific land cover types using the entropy-based local indicator of spatial association (ELSA). The values of ELSA vary between 0 and 1 (multiplied by a scale factor of 10000) denoting lowest and highest fragmentation. We use an existing spatially contiguous, global remote-sensing data product (i.e. the 27-year annual ESA CCI land cover maps which can be categorized as an EBV ‘Ecosystem Distribution’) to derive an annual (27 year) time-series of the Relative Magnitude of Fragmentation (RMF) at a global scale and with a spatial resolution of 300m. From this derived EBV data product, we can calculate a RMF indicator of ecosystem degradation, i.e. the change, and rate of change, in fragmentation of ecosystems (e.g. forests) over the last 27 years. This can provide important information for measuring biodiversity change as it directly links to the draft monitoring framework of the zero draft of the post-2020 global biodiversity framework of the Convention on Biological Diversity (CBD), especially Draft 2050 Goal 1 and the related Draft 2030 Target 1 (see Annex of the zero draft). We define the ‘forest’ class by aggregating all 14 tree cover related land cover types from the ESA CCI product into one class. We further define eight non-forest classes (agriculture, grassland, wetland, settlement, sparse vegetation, bare area, water, permanent snow and ice) that we use as multinomial categorical data, or as binary categorical data (to define forest vs. non-forest). This classes follow the reclassification used by Mousivand & Arsanjani 2019 (Applied Geography 106: 82-92). For deriving the RMF, we either calculate ELSA using the binary categorical data (forest vs. non-forest) or the multinomial categorical data (forest vs. the eight non-forest classes).
earthenv_topography_derived EarthEnv - Derived topographic categorical variables A fully standardized and global multivariate product of different terrain features to support many large-scale research applications. The product is based on the digital elevation model products of global 250 m GMTED2010. From https://www.earthenv.org/topography
distance_to_roads Distance to roads Distance to all roads in the Global Roads Inventory Project (GRIP) dataset (see https://www.globio.info/download-grip-dataset). The resolution is 300 meters at the equator and the units are the number of 300m resolution pixels. To get a distance in meters at the equator, multiply the value by 300.
cec_land_cover CEC North American Land Cover Land cover map generated by the Commission for Environmental Cooperation (CEC) for North America at 30 meter resolution. The map contains 19 land cover categories. See http://www.cec.org/north-american-environmental-atlas/land-cover-30m-2020/.
gbif_heatmaps Occurrence density maps created from GBIF data Density maps of all occurrences in GBIF, organized by taxonomic group at an approximately 1 km resolution.
cec_land_cover_percentage CEC North American Land Cover Percentage at 300 m resolution Percentager of each land cover at 300 m. reslution from land cover map generated by the Commission for Environmental Cooperation (CEC) for North America at 30 meter resolution. The collection contains one map for each of the 19 land cover categories. See http://www.cec.org/north-american-environmental-atlas/land-cover-30m-2020/.
cec_derived_maps Rasters derived from CEC North American Land Cover map to serve as inputs for various pipelines. Rasters at 300 m. resolution derived from the land cover map generated by the Commission for Environmental Cooperation (CEC) for North America at 30 meter resolution. See http://www.cec.org/north-american-environmental-atlas/land-cover-30m-2020/. These maps serve as inputs for various pipelines for species distribution and EBV modeling.
soilgrids Soil Grids datasets SoilGrids aggregated datasets at 1km resolution, from https://www.isric.org/explore/soilgrids
ouranos_climate_projections Ouranos projections climatiques Variables bioclimatiques projetées dérivées des projections Ouranos de la réanalyse canadienne de surface de l’ensemble de simulations ESPO-G6-R2 v1.0. Les valeurs présentées sont les percentiles ou les moyennes des 15 modèles utilisés. ESPO-G6-R2 v1.0 est un ensemble de simulations climatiques multimodèles issues du Projet d’intercomparaison de modèles couplés (CMIP6). Un post-traitement par quantile avec conservation de la tendance est appliqué à cet ensemble climatique afin de réduire les biais et faire la mise à l’échelle statistique des sorties de modèles globaux. ESPO-G6-R2 v1.0 est composé de séries quotidiennes sur une échelle de 10X10 km pour trois variables : température maximale quotidienne, température minimale quotidienne et précipitation quotidienne. Les données des modèles climatiques globaux sont post-traitées en utilisant comme cible la Réanalyse Canadienne de Surface (RCaSv2.1) d’Environnement et Changement climatique Canada. L’ensemble ESPO-G6-R2 v1.0 est disponible de 1950 à 2100 sur l’Amérique du Nord [https://www.ouranos.ca/fr/donnees-climatiques-ouranos/espo-g6-r2]. [Lavoie, J., Bourgault, P., Smith, T.J. et al. An ensemble of bias-adjusted CMIP6 climate simulations based on a high-resolution North American reanalysis. Sci Data 11, 64 (2024). [https://doi.org/10.1038/s41597-023-02855-z]]
aafc_land_use Cartes d’utilisation des terres d’Agriculture et Agroalimentaire Canada L’utilisation des terres d’AAC (2020) représente une méta-analyse culminante et organisée de plusieurs ensembles de données spatiales de haute qualité produits entre 1990 et 2021 par diverses méthodes employées par des équipes de chercheurs au gré de l’évolution des techniques et des capacités. L’information des ensembles de données d’entrée a été regroupée et intégrée à chaque pixel de 30 x 30 m pour composer des historiques sur pixels et dégager des milliers de combinaisons uniques de données prêtes à être examinées avec soin. Au moyen de nombreuses sources de données probantes de grande qualité et d’observations visuelles d’imagerie Google Earth, nous appliquons une stratégie progressive pour parvenir à la meilleure compréhension cohérente et actuelle de ce qui s’est passé dans chaque pixel pendant la série chronologique. https://ouvert.canada.ca/data/fr/dataset/7a098ea9-cc31-4d79-b326-89f6cd1fbb7d
geologie_du_socle SIGÉOM Géologie du socle Affinitées géologiques sélectionnées à partir des données de géologie du socle de SIGÉOM (MINISTÈRE DES RESSOURCES NATURELLES ET DES FORÊTS. Géologie du socle, [Jeu de données], dans Données Québec, 2018, mis à jour le 13 février 2026. [https://www.donneesquebec.ca/recherche/dataset/geologie-du-socle], (consulté le 27 février 2026)).
hirst-pearson_2022_canadian_human_footprint The Canadian Human Footprint (Hirst-Pearson, 2022) Canada’s human footprint reveals large intact areas juxtaposed against areas under immense anthropogenic pressure https://doi.org/10.1139/facets-2021-0063
mhc Modèle de hauteur de canopée Modèle de hauteur de canopée dérivé du lidar. Les données originales ont une résolution de 1 m. Une hauteur de 0 a été attribuée aux pixels sans valeur. MINISTÈRE DES RESSOURCES NATURELLES ET DES FORÊTS. Lidar - Modèles numériques (terrain, canopée, pente, courbe de niveau), [Jeu de données], dans Données Québec, 2016, mis à jour le 28 janvier 2025. [https://www.donneesquebec.ca/recherche/dataset/produits-derives-de-base-du-lidar], (consulté le 12 mai 2025)
twi Indice d’humidité topographique Indice d’humidité topographique (Topographic Wetness Index) dérivé du lidar. Les données originales ont une résolution de 1 m. Une valeur de 24 a été attribuée aux pixels sans valeur. MINISTÈRE DES RESSOURCES NATURELLES ET DES FORÊTS. Indice d’humidité topographique issu du LiDAR, [Jeu de données], dans Données Québec, 2020, mis à jour le 28 janvier 2025. [https://www.donneesquebec.ca/recherche/dataset/indice-humidite-topographique-issu-du-lidar], (consulté le 12 mai 2025).
sigeom_zones_morphosedimentologiques_percentage SIGÉOM Zones morphosédimentologiques pourcentage Zones morphosédimentologiques de la géologie du Quaternaire. Provient d’une rasterization de fichiers vectoriels et d’une convesion en pourcentage. Voir https://sigeom.mines.gouv.qc.ca/signet/classes/I1102_aLaCarte#QUAT.
canadian_human_footprint The Canadian Human Footprint Canada’s human footprint reveals large intact areas juxtaposed against areas under immense anthropogenic pressure https://doi.org/10.1139/facets-2021-0063
pither_functional_connectivity Functional landscape connectivity In this study we tested the hypothesis that functional connectivity for multiple species can be estimated across Canada using a single, upstream connectivity model. We developed a movement cost layer with cost values assigned using expert opinion to anthropogenic land cover features and natural features based on their known and assumed effects on the movement of terrestrial, non-volant fauna. We used Circuitscape to conduct an omnidirectional connectivity analysis for terrestrial landscapes, in which the potential contribution of all landscape elements to connectivity were considered and where source and destination nodes were independent of land tenure. Our resulting map of mean current density provided a seamless estimate of movement probability at a 300 m resolution across Canada. https://doi.org/10.1371/journal.pone.0281980
sothe_2022_canadian_carbon_stocks Carbon storage and distribution in terrestrial ecosystems of Canada This collection contains the maps produced in the project Carbon storage and distribution in terrestrial ecosystems of Canada. We estimated carbon stocks in plants (AGB, BGB, dead plants) and soils (0-1m depth) of Canada at 250 m spatial resolution. We used multisource satellite, climate and topographic data and a machine-learning algorithm. https://doi.org/10.4121/c.5421810
ndvi_max_qc NDVImax sur la province de Québec Nous avons utilisé toutes les données de niveau 1 (Tier 1) disponibles pour la province du Québec entre 1985 et 2024 pour les mois de juin, juillet, août et septembre dans la collection Landsat Collection 2 fournie par l’U.S. Geological Survey (USGS) et hébergée sur Google Earth Engine (GEE). Les produits de données de niveau 1 analysés incluent la réflectance de surface des capteurs Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper + (ETM+), et Landsat 8 Operational Land Imager (OLI). Nous avons exclusivement sélectionné des images ayant une couverture nuageuse moyenne inférieure à 80 %, car des scènes à forte couverture nuageuse peuvent compromettre la précision de la calibration géométrique. La version C de Function of Mask (CFmask) a été appliquée pour classer chaque pixel comme clair (terre/eau), neige, nuage, adjacent à un nuage, ou ombre de nuage (Zhu et al., 2015; Zhu & Woodcock, 2014). Les pixels affectés par la neige, les nuages, les zones adjacentes aux nuages ou les ombres de nuages ont été exclus de l’analyse. Les effets angulaires résultant de la variation de la géométrie d’observation et solaire tout au long de la série temporelle Landsat ont été identifiés comme une source significative de variation dans la réflectance directionnelle obtenue, indépendamment des changements liés au sol (Nagol et al., 2015). La correction de ces effets est cruciale en raison de la variabilité de la géométrie capteur-surface-soleil, décrite comme la Fonction de Répartition Directionnelle de la Réflectance (BRDF), au moment de l’acquisition, qui varie à la fois spatialement et temporellement. Roy et al. (2016) ont introduit une méthode complète utilisant un ensemble fixe de paramètres basés sur le modèle BRDF RossThick-LiSparse (Schaaf et al., 2002). Cette méthode permet de normaliser toute la série temporelle Landsat au nadir (angle de zénith de visionnage à 0°) et à une SZA constante, réduisant les variations de réflectance liées à la BRDF. Nous avons appliqué cette correction directement dans GEE en utilisant le script développé par Erik Lindquist et Daniel Wiell, disponible dans les groupes des développeurs Earth Engine. Notre implémentation de la correction a impliqué l’utilisation de coefficients globaux et d’une SZA normalisée optimale et constante par emplacement (Zhang et al., 2016). Berner et al. (2020) ont identifié des divergences radiométriques systématiques entre les satellites Landsat 5 TM, 7 ETM+ et 8 OLI. Ils ont proposé une méthode basée sur la proximité phénologique pour effectuer une calibration croisée des satellites Landsat 5 vers Landsat 7, et Landsat 8 vers Landsat 7, en utilisant les années où les deux satellites étaient disponibles. Landsat 7 a servi de référence en raison de son chevauchement temporel avec les deux autres capteurs. Nous avons effectué une calibration croisée des capteurs pour corriger les divergences radiométriques systématiques dans les bandes spectrales des satellites Landsat 5 TM, 7 ETM+ et 8 OLI en utilisant les méthodes et coefficients pour l’Arctique de Berner et al. (2023). Nous avons calculé le NDVI en utilisant les bandes Rouge et NIR corrigées pour la BRDF et la calibration croisée des capteurs, puis calculé le maximum annuel pendant la saison de croissance. Berner et al. (2020) ont démontré que les estimations de NDVImax dépendent du nombre d’observations durant la saison de croissance. Bayle et al. (2024) ont montré que, comme les observations Landsat augmentent au fil de la série temporelle, le maximum de NDVI en début de série a tendance à être sous-estimé comparé aux années récentes. Berner et al. (2023) ont montré que ce biais d’échantillonnage peut être partiellement corrigé en modélisant la phénologie sur une base pixel par pixel et en ajustant les estimations de NDVImax en fonction de la phénologie moyenne. Par conséquent, nous avons appliqué une modélisation phénologique pour ajuster les valeurs de NDVI avant de calculer le NDVImax. Nous avons utilisé la méthode de reconstruction par Analyse Harmonique des Séries Temporelles (HANTS) comme utilisée par Choler et al. (2024). Cette procédure n’a été appliquée qu’aux observations de NDVI avec des valeurs supérieures à 0,15 pour éliminer les pixels non végétalisés, car la correction repose sur la modélisation phénologique. La procédure complète, reposant sur des méthodes à la pointe, nous a permis de dériver une série robuste de NDVImax annuel pour la province du Québec. Le jeu de données est constitué de 40 rasters sous le format GEOTIFF avec un raster par année contenant le maximum de NDVI. Pour certaines années où le nombre d’observation est limitée, des manques de données de NDVImax sont visibles et signalés par des NA. Par: Arthur Bayle. Code: https://github.com/arthurbayle/quebec_ndvimax
inat_canada_heatmaps iNaturalist Canada Occurrence density maps Density maps of all occurrences in iNaturalist Canada, organized by taxonomic group at an approximately 100m or 1km resolution.
mhp2023 Milieux humides potentiels 2023 La cartographie des milieux humides potentiels du Québec (CMHPQ) 2023. Voir https://www.donneesquebec.ca/recherche/dataset/milieux-humides-potentiels.
sigeom_zones_morphosedimentologiques SIGÉOM Zones morphosédimentologiques Zones morphosédimentologiques de la géologie du Quaternaire. Provient d’une rasterization de fichiers vectoriels. Voir https://sigeom.mines.gouv.qc.ca/signet/classes/I1102_aLaCarte#QUAT.
mhp2023_percentage Milieux humides potentiels 2023 pourcentage La cartographie des milieux humides potentiels du Québec (CMHPQ) 2023 transformée en pourcentage à 100 m. Voir https://www.donneesquebec.ca/recherche/dataset/milieux-humides-potentiels.
ouranos_past_climate_monthly Ouranos climat passé mensuel Variables climatiques mensuelles passée calculées à partir de xclim [https://xclim.readthedocs.io/en/stable/] et basées sur les données Ouranos [Lavoie, J., Bourgault, P., Smith, T.J. et al. An ensemble of bias-adjusted CMIP6 climate simulations based on a high-resolution North American reanalysis. Sci Data 11, 64 (2024). https://doi.org/10.1038/s41597-023-02855-z]
verdissement_qc Verdissement sur la province de Québec Cette carte représente le taux de verdissement calculé à partir des images Landsat de 1984 à 2024 et utilisées pour obtenir des cartes de NDVI max pour chaque année. Ensuite, un pente a été calculée pour chaque pixel et la valeur de la pente est représentée sur cette carte. Par: Arthur Bayle. Code: https://github.com/arthurbayle/quebec_ndvimax
ouranos_projected_climate_models Ouranos projections climatiques détaillées Variables bioclimatiques projetées dérivées des projections Ouranos de la réanalyse canadienne de surface de l’ensemble de simulations ESPO-G6-R2 v1.0. ESPO-G6-R2 v1.0 est un ensemble de simulations climatiques multimodèles issues du Projet d’intercomparaison de modèles couplés (CMIP6). Un post-traitement par quantile avec conservation de la tendance est appliqué à cet ensemble climatique afin de réduire les biais et faire la mise à l’échelle statistique des sorties de modèles globaux. ESPO-G6-R2 v1.0 est composé de séries quotidiennes sur une échelle de 10X10 km pour trois variables : température maximale quotidienne, température minimale quotidienne et précipitation quotidienne. Les données des modèles climatiques globaux sont post-traitées en utilisant comme cible la Réanalyse Canadienne de Surface (RCaSv2.1) d’Environnement et Changement climatique Canada. L’ensemble ESPO-G6-R2 v1.0 est disponible de 1950 à 2100 sur l’Amérique du Nord [https://www.ouranos.ca/fr/donnees-climatiques-ouranos/espo-g6-r2]. [Lavoie, J., Bourgault, P., Smith, T.J. et al. An ensemble of bias-adjusted CMIP6 climate simulations based on a high-resolution North American reanalysis. Sci Data 11, 64 (2024). [https://doi.org/10.1038/s41597-023-02855-z]]
carte_eco_code_terrain_depot Code terrain et dépôt de surface sélectionnés (Carte écoforestière) Distance par rapport à différentes entités sélectionnées à partir des codes de terrain et des dépôts de surface de la carte écoforestière provinciale à jour en format GPKG [https://diffusion.mffp.gouv.qc.ca/Diffusion/DonneeGratuite/Foret/DONNEES_FOR_ECO_SUD/Cartes_ecoforestieres_perturbations/02-Donnees/PROV/CARTE_ECO_MAJ_PROV_GPKG.zip]. Dernière mise à jour de la carte écoforestière le 2025-07-14. Accédé le 2026-02-25.
cop-dem-glo Elevation and derived variables (Copernicus DEM) The Copernicus DEM is a digital surface model (DSM), which represents the surface of the Earth including buildings, infrastructure, and vegetation. This DSM is based on radar satellite data acquired during the TanDEM-X Mission, which was funded by a public-private partnership between the German Aerospace Centre (DLR) and Airbus Defence and Space. Copernicus DEM is available at both 30-meter and 90-meter resolution; this dataset has a horizontal resolution of approximately 30 meters. See the dataset page on OpenTopography: [https://doi.org/10.5069/G9028PQB]. (European Space Agency (2024). Copernicus Global Digital Elevation Model. Distributed by OpenTopography. https://doi.org/10.5069/G9028PQB. Accessed 2026-02-21. The DEM source is © DLR e.V. 2010-2014 and © Airbus Defence and Space GmbH 2014-2018 provided under COPERNICUS by the European Union and ESA; all rights reserved. The derived variables were produced using Copernicus WorldDEM-30 © DLR e.V. 2010-2014 and © Airbus Defence and Space GmbH 2014-2018 provided under COPERNICUS by the European Union and ESA; all rights reserved)
geomorphons_percentages Geomorphons percentages derived from Copernicus DEM Geomorphons percentages derived from the 30m Copernicus Global Digital Elevation Model (DEM) (https://doi.org/10.5069/G9028PQB) using WhiteboxTools (https://www.whiteboxgeo.com/). Copernicus Global DEM acquired from the Planetary Computer. The derived variables were produced using Copernicus WorldDEM-30 © DLR e.V. 2010-2014 and © Airbus Defence and Space GmbH 2014-2018 provided under COPERNICUS by the European Union and ESA; all rights reserved.
sigeom_zones_morphosedimentologiques_combined SIGÉOM Zones morphosédimentologiques combinées Zones morphosédimentologiques de la géologie du Quaternaire. Provient d’une rasterization de fichiers vectoriels. Voir https://sigeom.mines.gouv.qc.ca/signet/classes/I1102_aLaCarte#QUAT.
geomorphons Geomorphons derived from Copernicus DEM Geomorphons derived from the 30m Copernicus Global Digital Elevation Model (DEM) (https://doi.org/10.5069/G9028PQB) using WhiteboxTools (https://www.whiteboxgeo.com/). Copernicus Global DEM acquired from the Planetary Computer. The derived variables were produced using Copernicus WorldDEM-30 © DLR e.V. 2010-2014 and © Airbus Defence and Space GmbH 2014-2018 provided under COPERNICUS by the European Union and ESA; all rights reserved.
sigeom_zones_morphosedimentologiques_percentage_combined SIGÉOM Zones morphosédimentologiques pourcentage combinées Zones morphosédimentologiques de la géologie du Quaternaire. Provient d’une rasterization de fichiers vectoriels et d’une convesion en pourcentage. Voir https://sigeom.mines.gouv.qc.ca/signet/classes/I1102_aLaCarte#QUAT.
GRHQ GRHQ Distances au réseau hydrographique du Québec Distances (m) au réseau hydrique pour différents éléments extraits et simplifiés de la Géobase du réseau hydrographique du Québec (GRHQ). Les lacs incluent les plans d’eau où il n’y a pas d’écoulement (lacs, réservoirs, étangs, mares, parfois constructions artificielles). Les rivières considérées sont les sections représentées par des polygones et sont donc des rivières d’une certaine taille. La distance aux ruisseaux concernent la distance au réseau linéaire ce qui inclus les ruisseaux, les rivières, les canaux, certains fossés et le Saint-Laurent. Les milieux humides ne concernent que les milieux humides du sud du Québec qui sont répertoriés dans la base de données. Finalement, la distance au St-Laurent prend également en compte certains affluents principaux (Outaouais, Richelieu). MINISTÈRE DES RESSOURCES NATURELLES ET DES FORÊTS. Géobase du réseau hydrographique du Québec (GRHQ), [Jeu de données], dans Données Québec, 2019, mis à jour le 23 octobre 2025. [https://www.donneesquebec.ca/recherche/dataset/grhq], (consulté le 27 novembre 2025).
ouranos_past_climate_period Ouranos climat passé par période Climat passé par période de 30 ans Ouranos concernant les 19 variables bioclimatiques. [Lavoie, J., Bourgault, P., Smith, T.J. et al. An ensemble of bias-adjusted CMIP6 climate simulations based on a high-resolution North American reanalysis. Sci Data 11, 64 (2024). https://doi.org/10.1038/s41597-023-02855-z]
ouranos_past_climate_yearly Ouranos climat passé annuel Climat annuel passé Ouranos concernant les 19 variables bioclimatiques. [Lavoie, J., Bourgault, P., Smith, T.J. et al. An ensemble of bias-adjusted CMIP6 climate simulations based on a high-resolution North American reanalysis. Sci Data 11, 64 (2024). https://doi.org/10.1038/s41597-023-02855-z]
test_collection This is a test collection This is a test collection description

Search for a specific collection (earthenv_landcover)

it_obj <- s_obj |>
  stac_search(collections = "earthenv_landcover") |>
  post_request() |> items_fetch()
it_obj
## ###Items
## - features (12 item(s)):
##   - class_9
##   - class_8
##   - class_7
##   - class_6
##   - class_5
##   - class_4
##   - class_3
##   - class_2
##   - class_12
##   - class_11
##   - class_10
##   - class_1
## - assets: data
## - item's fields: 
## assets, bbox, collection, geometry, id, links, properties, stac_extensions, stac_version, type

See the available layers in this collection

it_obj <- s_obj |>
  collections("earthenv_landcover") |> items() |>
  get_request() |> items_fetch()
it_obj
## ###Items
## - features (12 item(s)):
##   - class_9
##   - class_8
##   - class_7
##   - class_6
##   - class_5
##   - class_4
##   - class_3
##   - class_2
##   - class_12
##   - class_11
##   - class_10
##   - class_1
## - assets: data
## - item's fields: 
## assets, bbox, collection, geometry, id, links, properties, stac_extensions, stac_version, type

See the properties of the first item (layer) in the collection

it_obj[['features']][[1]]$properties
## $class
## [1] "9"
## 
## $datetime
## [1] "2000-01-01T00:00:00Z"
## 
## $`proj:epsg`
## [1] 4326
## 
## $description
## [1] "Urban/Built-up"
## 
## $full_filename
## [1] "consensus_full_class_9.tif"

Summary of items

df<-data.frame(id=character(),datetime=character(), description=character())
for (f in it_obj[['features']]){
  df<-rbind(df,data.frame(id=f$id,datetime=f$properties$datetime,description=f$properties$description))
}
kable(df)
id datetime description
class_9 2000-01-01T00:00:00Z Urban/Built-up
class_8 2000-01-01T00:00:00Z Regularly Flooded Vegetation
class_7 2000-01-01T00:00:00Z Cultivated and Managed Vegetation
class_6 2000-01-01T00:00:00Z Herbaceous Vegetation
class_5 2000-01-01T00:00:00Z Shrubs
class_4 2000-01-01T00:00:00Z Mixed/Other Trees
class_3 2000-01-01T00:00:00Z Deciduous Broadleaf Trees
class_2 2000-01-01T00:00:00Z Evergreen Broadleaf Trees
class_12 2000-01-01T00:00:00Z Open Water
class_11 2000-01-01T00:00:00Z Barren
class_10 2000-01-01T00:00:00Z Snow/Ice
class_1 2000-01-01T00:00:00Z Evergreen/Deciduous Needleleaf Trees

Access the first item with the Terra package. Here, we will select the layer representing the percentage of “Evergreen/Deciduous Needleleaf Trees”.

library(stars)
lc1<-rast(paste0('/vsicurl/',it_obj[['features']][[12]]$assets$data$href))
plot(lc1)

Select on only a section of the raster

bbox<-st_bbox(c(xmin = -76, xmax = -70, ymax = 54, ymin = 50), crs = st_crs(4326))
lc2 <- lc1 |> crop(bbox)

Visualise it

pal <- colorRampPalette(c("black","darkblue","red","yellow","white"))
plot(lc2,breaks=seq(0,100,10),col=pal(10))

Save it on your computer in Cloud Optimized GeoTiff format.

writeRaster(lc2,'~/lc3.tif',filetype='COG',gdal=c('COMPRESS=DEFLATE'), overwrite=TRUE)

Note that for a variable with categorical values, saving is a bit more complex.

lc1 |> crop(bbox) |> writeRaster('~/lc1.tif',filetype='COG',gdal=c('COMPRESS=DEFLATE','OVERVIEW_RESAMPLING=MODE','OVERVIEWS=IGNORE_EXISTING', overwrite=TRUE))

Use of GDALCUBES

Filtrer en fonction des propriétés des items et créer une collection Filter items based on their properties and create a collection. Here, we will select only the layers representing the percentage of “Evergreen/Deciduous Needleleaf Trees”, “Evergreen Broadleaf Trees”, “Deciduous Needleleaf Trees” and “Deciduous Broadleaf Trees”.

st <- stac_image_collection(it_obj$features, asset_names=c('data'), property_filter = function(f){f[['class']] %in% c('1','2','3','4')},srs='EPSG:4326')
st
## Image collection object, referencing 4 images with 1 bands
## Images:
##      name left top bottom right            datetime       srs
## 1 class_4 -180  90    -56   180 2000-01-01T00:00:00 EPSG:4326
## 2 class_3 -180  90    -56   180 2000-01-01T00:00:00 EPSG:4326
## 3 class_2 -180  90    -56   180 2000-01-01T00:00:00 EPSG:4326
## 4 class_1 -180  90    -56   180 2000-01-01T00:00:00 EPSG:4326
## 
## Bands:
##   name offset scale unit nodata image_count
## 1 data      0     1                       4

Construct a cube to process or visualize the data. Note that this cube can be in a different CRS and resolution than the original items/files. However, the temporal dimension must capture the temporal frame of the item. dt is expressed as a time period. P1D is a period of one day, P1M is a period of one month, P1Y is a period of one year. Resampling methods should be adapted to the type of data. For categorical data, use “mode” or “nearest”. For continuous data, use “bilinear”. Aggregation is only relevant when multiple rasters overlap.

Here, we will sum the four forest categories using aggregation=“sum”. We will also change the reference system of the data to use Quebec Lambert (EPSG:32198) and set the resolution to 1km.

bbox<-st_bbox(c(xmin = -483695, xmax = -84643, ymin = 112704 , ymax = 684311), crs = st_crs(32198))

v <- cube_view(srs = "EPSG:32198", extent = list(t0 = "2000-01-01", t1 = "2000-01-01",
                                                left = bbox$xmin, right =bbox$xmax, top = bbox$ymax, bottom = bbox$ymin), dx=1000, dy=1000, dt="P1D", aggregation = "sum", resampling = "mean")

Combine the collection and the cube_view to create a raster cube. This cube can then be processed or visualized. Note that the resulting cube is not a raster file, but an object in R that can be manipulated with the gdalcubes package. The data is only read and processed when necessary, so it is possible to work with large datasets without having to download them entirely.

lc_cube <- raster_cube(st, v)

Save the resulting file on your computer.

lc_cube |> write_tif('~/',prefix='lc2',COG=TRUE, creation_options=list('COMPRESS'='DEFLATE'))
lc_cube |> plot(zlim=c(0,100),col=pal(10))

Use the dataset “Accessibility from cities”, keeping the same CRS and extent. Here, we will visualize the layer representing the travel time to the nearest city with more than 10,000 inhabitants. Note that this layer is in minutes, so we will use a different color palette and breaks for visualization.

it_obj <- s_obj |>
  collections("accessibility_to_cities") |> items() |>
  get_request() |> items_fetch()
v <- cube_view(srs = "EPSG:32198", extent = list(t0 = "2015-01-01", t1 = "2015-01-01",
                                                left = bbox$xmin, right =bbox$xmax, top = bbox$ymax, bottom = bbox$ymin), dx=1000, dy=1000, dt="P1D", aggregation = "mean", resampling = "bilinear")
for (i in 1:length(it_obj$features)){
  it_obj$features[[i]]$assets$data$roles='data'
}
st <- stac_image_collection(it_obj$features)
lc_cube <- raster_cube(st, v)
lc_cube |> plot(col=heat.colors)

Use the CHELSA dataset on climatologies and create a map of the means for the months of June, July and August from 2010 to 2019. Here, we will visualize the layer representing the average temperature (tas) for these months and years. Note that this layer is in degrees Celsius multiplied by 10.

it_obj <- s_obj |>
  stac_search(collections = "chelsa-monthly", datetime="2010-06-01T00:00:00Z/2019-08-01T00:00:00Z") |> get_request() |> items_fetch()

v <- cube_view(srs = "EPSG:32198", extent = list(t0 = "2010-06-01", t1 = "2019-08-31",
                                                left = bbox$xmin, right =bbox$xmax, top = bbox$ymax, bottom = bbox$ymin),
               dx=1000, dy=1000, dt="P10Y",
               aggregation = "mean",
               resampling = "bilinear")

for (i in 1:length(it_obj$features)){
  names(it_obj$features[[i]]$assets)='data'
  it_obj$features[[i]]$assets$data$roles='data'
}
anames=unlist(lapply(it_obj$features,function(f){f['id']}))
st <- stac_image_collection(it_obj$features, asset_names = 'data',  property_filter = function(f){f[['variable']] == 'tas' & (f[['month']] %in% c(6,7,8)) })
c_cube <- raster_cube(st, v)
c_cube |> plot(col=heat.colors)

c_cube |> write_tif('~/',prefix='chelsa-monthly',creation_options=list('COMPRESS'='DEFLATE'))