Méthodes de cartographie

Dean Evans and David Iles

Résumé

Afin de générer des cartes normalisées à l’échelle de la province illustrant l’abondance relative des oiseaux nicheurs en Saskatchewan, nous avons réalisé une modélisation statistique à partir des données recueillies dans le cadre de l’Atlas des oiseaux nicheurs de la Saskatchewan. Étant donné que la couverture et l’effort de recensement variaient d’une parcelle à l’autre et que de nombreuses parcelles situées au nord n’avaient pas été recensées, nous avons utilisé la modélisation statistique pour corriger ces variations et estimer l’abondance dans les zones non échantillonnées.

Notre approche comprenait quatre étapes principales :

1. Filtrage des données spécifiques à chaque espèce appropriées pour l’analyse de l’abondance ;
2. Ajustement de modèles statistiques pour estimer les tendances spatiales d’abondance et les associations d’habitats, en tenant compte des différences dans l’effort et le moment des enquêtes ;
3. Générer des prévisions basées sur des modèles concernant l’abondance relative et la probabilité d’observation à travers la province, pour un effort de recensement standardisé ;
4. Examiner et valider les résultats du modèle, tant sur le plan quantitatif que qualitatif, afin d’en garantir la fiabilité.

Toutes les analyses ont été réalisées avec la version 4.4.0 de R. Le code et la documentation sont disponibles à l’adresse suivante : https://github.com/deanrobertevans/SaskAtlas/.

Sources de données et filtrage

Nous avons utilisé trois sources principales de données d’observation des espèces :

• Points d’écoute : points d’écoute standardisés de 3, 5 ou 10 minutes au cours desquels les observateurs ont enregistré tous les oiseaux qu’ils ont vus ou entendus ;
• Enregistrements acoustiques : d’une durée maximale de 10 minutes, transcrits par des experts afin d’obtenir un décompte des espèces similaire à celui obtenu par points d’écoute (mais sans détection visuelle) ;
• Dénombrements stationnaires et transects linéaires : listes complétées dans chaque parcelle de 10 × 10 km indiquant les espèces détectées et celles qui ne l’ont pas été au cours de l’inventaire. Seules les listes complètes (c’est-à-dire celles dans lesquelles les observateurs ont indiqué toutes les espèces observées) ont été incluses. Les listes de transects linéaires indiquent que le participant à l’atlas s’est déplacé dans la parcelle tout en enregistrant les oiseaux et comprennent des mesures associées de la distance parcourue et du temps écoulé. Les dénombrements stationnaires sont similaires aux points d’écoute, mais les durées ne sont pas standardisées..

Toutes les données ont été extraites des plateformes NatureCounts et WildTrax, puis recoupées afin d’éliminer les duplications. Les recensements ont été inclus dans l’analyse de l’abondance s’ils :

• A eu lieu pendant la saison de reproduction de l’espèce considérée et pendant la période où celle-ci était censée être la plus active vocalement ;
• Étaient référencées spatialement et temporellement (c’est-à-dire qu’elles comprenaient les coordonnées GPS et les informations de date/heure) ;
• Représentaient des protocoles d’enquête valides et complets (par exemple, des points d’écoute de 3 à 10 minutes ou des listes complètes) ;
• Respectaient les exigences minimales en matière d’assurance qualité (par exemple, si elles n’étaient pas signalées comme aberrantes).

Même si un recensement n’a pas été utilisé pour modéliser l’abondance relative, il a tout de même contribué à d’autres produits de l’atlas, comme les résumés sur la présence et la nidification des espèces.

Covariables

Nous avons modélisé l’abondance des oiseaux à l’aide d’un ensemble de covariables environnementales, résumées dans un rayon d’un kilomètre autour de chaque site de recensement. Ces covariables comprenaient :

• Couverture terrestre : Pourcentage de couverture « urbaine », « agricole », « humide », « forestière à feuilles aiguillées » et « herbeuse/arbustive » tiré de l’ensemble de données sur la couverture terrestre du Canada en 2020 (Ressources naturelles Canada, 2020) ;
• Caractéristiques forestières : pourcentage de couverture forestière de conifères et de feuillus, ainsi que fermeture de la canopée, à partir de l’ensemble de données SCANFI (Guindon et al., 2024) ;
• Topographie : altitude moyenne et pente dérivées des cartes numériques d’altitude de la Saskatchewan obtenues à l’aide du paquet R « elevatr » (Hollister 2023) ;
• Climat : température annuelle moyenne et précipitations (normales 1991-2020 ; Fick & Hijmans 2017 ; Wang et al. 2016) ;
• Hydrologie : distance moyenne par rapport à l’eau et récurrence des zones humides (Pekel et al. 2016).

Notre analyse comprenait au total 14 variables biophysiques. Afin de réduire la multicolinéarité et de simplifier l’interprétation, nous avons effectué une analyse en composantes principales (ACP) de ces variantes environnementales. Nous avons retenu les sept composantes principales (CP) les plus importantes comme prédicteurs dans les modèles, qui expliquaient ensemble 84 % de la variation de l’habitat entre les sites. Chaque CP décrivait différents gradients d’habitat dans la zone d’étude. L’axe CP 1, par exemple, décrivait la transition entre une forêt dense (canopée dense, pourcentage élevé de feuillus et de conifères, température plus basse) et des paysages ouverts (canopée clairsemée, couverture agricole et prairiale élevée, température annuelle plus élevée). Le deuxième axe CP décrivait le contraste entre les zones à forte couverture agricole et celles à forte couverture herbacée/arbustive. Le troisième axe CP décrivait principalement la « récurrence de l’humidité » du paysage, avec un gradient allant des zones basses à forte humidité et récurrence des zones humides, vers les zones à haute altitude et faible récurrence des zones humides. Les autres axes CP avaient des interprétations plus complexes qui décrivaient d’autres gradients d’habitats à travers les sites d’étude (par exemple, les gradients de précipitations, le développement urbain, etc.). Tous les prédicteurs continus ont été modélisés à l’aide de termes quadratiques, ce qui a permis d’établir des relations non linéaires et intermédiaires optimales lorsque les données le permettaient.

Nous avons également inclus des covariables liées à l’observateur afin de prendre en compte les variations spatio-temporelles de l’effort et du calendrier de l’enquête. Ces corrections relatives à l’effort et au moment de l’enquête comprenaient :

1. Une variable « temps écoulé depuis le lever du soleil » (TSS) est incluse à chaque relevé afin de tenir compte des variations d’activité vocale des oiseaux au cours d’une journée. Cette variable est associée à un effet spline curvilinéaire flexible sur l’abondance relative.
2. Pour les points d’écoute et les enregistrements acoustiques, nous avons inclus des compensations de détectabilité calculées à l’aide de l’approche QPAD (Solymos et al., 2013), qui tiennent compte de la durée des relevés (par exemple, 3, 5 ou 10 minutes) et du rayon de détection effectif des oiseaux. Nous avons utilisé les mêmes calculs de compensation pour les points d’écoute humains et les enregistrements acoustiques, car plusieurs études ont montré que les points d’écoute sont similaires pour ces deux méthodes (les enquêtes menées à l’aide d’unités d’enregistrement autonomes (ARU) sont largement équivalentes aux points d’écoute humains pour la plupart des espèces).
3. Pour les listes de transects linéaires, nous avons inclus des effets de splines curvilignes flexibles qui tiennent compte à la fois de la durée du relevé et de la distance parcourue.
4. Pour les dénombrements stationnaires, nous avons inclus un effet spline curviligne flexible de la durée de l’enquête.

Modélisation statistique de l’abondance relative

L’abondance relative de chaque espèce d’oiseau pendant la saison de nidification a été estimée à l’aide d’un cadre de modélisation spatiale intégrant simultanément plusieurs types de relevés (points d’écoute, relevés ARU, listes de transects linéaires et dénombrements stationnaires). Dans cette approche de modélisation conjointe, toutes les sources de données ont contribué à l’estimation des relations entre les espèces et leur habitat, ainsi qu’à la détermination des limites de leur aire de répartition et de leur structure spatiale. Cette approche a été mise en œuvre à l’aide du paquet R « inlabru » (Bachl et al., 2019), qui prend en charge les modèles bayésiens intégrés spatiaux à travers plusieurs vraisemblances partageant des composants de modèle communs.

Structure du modèle conjoint

Le modèle a pris en compte :

• Covariables spatiales : ensemble de variables environnementales supposées influencer la répartition des oiseaux et mesurées à l’aide des axes des composantes principales (voir ci-dessus).
• Autocorrélation spatiale : modélisée à l’aide de champs aléatoires gaussiens (CAG) en utilisant l’approche des équations différentielles partielles stochastiques (EDPS), ce qui permet au modèle de capturer les dépendances spatiales résiduelles dans les données qui n’étaient pas directement liées aux covariables de l’habitat.
• Corrections liées à l’effort : pour tenir compte des différences dans le dénombrement des oiseaux résultant d’études de durée différente (pour les dénombrements stationnaires) et de distance parcourue différente (pour les listes de transects linéaires).
• Le temps écoulé depuis le lever du soleil : modélisé à l’aide d’effets de type spline pour tenir compte des variations de détectabilité des oiseaux au cours d’une journée.
• Corrections de détectabilité : estimées à l’aide de l’échantillonnage par retrait temporel et et de l’échantillonnage à distance, selon l’approche QPAD décrite par Solymos et al. (2013) et Edwards et al. (2023).
• A priori informatifs : les limites estimées de l’aire de répartition des espèces provenant d’eBird (Fink et al. 2023) ont été utilisées comme a priori informatifs, fournissant des contraintes spatiales spécifiques à chaque espèce.

Les trois types de données (points d’écoute/dénombrements acoustiques, dénombrements stationnaires et transects linéaires) ont été reliés par une structure prédictive commune, dans laquelle le prédicteur linéaire de chaque ensemble de données comprenait :

• Une interception spécifique à l’ensemble des données.
• Effet de bordure de l’aire de répartition : modélisé comme une fonction linéaire de la distance par rapport à la bordure estimée de l’aire de reproduction de l’espèce.
• Covariables environnementales : basées sur les axes des composantes principales décrits ci-dessus, incluses sous forme de termes linéaires et quadratiques.
• Un champ aléatoire spatial partagé, estimé à l’aide d’une approche EDPS.

La variable de réponse et la distribution des erreurs utilisées pour la modélisation variaient en fonction de l’ensemble de données :

1. Points d’écoute et enregistrements acoustiques (PC) : modélisés sous forme de dénombrements à l’aide d’une erreur binomiale négative et d’une fonction de liaison logarithmique.
2. Listes de transects linéaires (LT) : modélisées sous forme d’enregistrements binaires présence/absence à l’aide d’une erreur binomiale et d’une fonction de liaison log-log complémentaire.
3. Les dénombrements stationnaires (SC) : également modélisés comme des données de présence/absence à l’aide d’une distribution binomiale et d’une fonction de liaison log-log complémentaire.

Bien que les listes de contrôle complètes indiquaient le nombre d’oiseaux, les analyses préliminaires ont révélé que l’inclusion des informations sur le nombre total provenant de ces listes introduisait une grande incertitude dans les prévisions et réduisait la précision de la validation croisée par rapport aux modèles qui n’utilisaient que les informations sur la présence/absence. Nous avons donc réduit les nombres à la présence/absence en considérant tout nombre supérieur à 0 comme une présence.

Équations de vraisemblance

La vraisemblance des dénombrements dérivés des points d’écoute et des enregistrements acoustiques a donc été modélisée comme suit :

Les termes de cette équation sont :

• : Interception globale pour les relevés de points d’écoute et enregistrements acoustiques.
• : Correction de compensation pour la détectabilité spécifique à l’espèce et à l’effort à partir des modèles QPAD.
• : Effet spline du temps écoulé depuis le lever du soleil pour relevé .
• : Coefficient permettant de saisir les changements d’abondance à proximité des limites de l’aire de répartition des espèces.
• : Distance par rapport à la limite de l’aire de répartition de l’espèce pour relevé .
• : Effet aléatoire spatial (champ gaussien latent) pour l’emplacement du relevé .
• : Effets fixes additifs pour les covariables environnementales standardisées et leurs termes quadratiques, pour chacun des axes des composantes principales inclus comme prédicteurs.

Les observations binaires de présence/absence issues des transects linéaires ont été modélisées à l’aide d’une fonction de vraisemblance de Bernoulli, avec un lien log-log complémentaire correspondant à un processus de détection de Poisson latent, comme suit :

Les termes supplémentaires dans cette équation sont :

• : Interception globale pour les listes de transects linéaires.
• : Logarithme de la durée de l’enquête comme compensation pour tenir compte de l’effort fourni.
• : Logarithme de la distance parcourue lors du relevé, utilisé comme compensation pour tenir compte de l’effort fourni.

De même, les dénombrements stationnaires ont été modélisés à l’aide de :

où :

• : Interception globale pour les dénombrements stationnaires
• : Logarithme de la durée de l’enquête (en minutes), tenant compte de l’effort fourni.

Validation et cartographie de l’abondance relative et de la probabilité d’observation

Pour produire des prévisions spatiales de l’abondance relative, nous avons généré des estimations à l’aide de modèles pour l’ensemble de la Saskatchewan, à l’aide d’une grille d’une résolution d’un kilomètre. Pour chaque espèce, nous avons prédit le nombre d’individus détectés au cours de 15 points d’écoute standardisés de cinq minutes par cellule de la grille, en nous basant sur les covariables télédétectées dans cette cellule. Ces prévisions, dérivées du modèle ajusté à l’aide d’un échantillonnage a posteriori (1 000 tirages), étaient initialement sur une échelle logarithmique. Nous avons ensuite exponentié ces prévisions et les avons résumées à l’aide de la médiane et de l’intervalle de crédibilité à 90 % (c’est-à-dire les 5^e et 95^e centiles), fournissant ainsi une mesure de la tendance centrale et de l’incertitude de l’abondance prévue. Nous avons également calculé la largeur de l’intervalle de crédibilité à 90 % afin d’indiquer le degré d’incertitude des prévisions.

Afin de rendre les prévisions plus compréhensibles pour les utilisateurs finaux, nous avons également calculé la probabilité d’observer l’espèce au moins une fois sur 15 points d’écoute, en nous basant sur l’abondance prévue et sur les paramètres ajustés du modèle binomial négatif. Cela a été calculé comme un moins la probabilité d’observer zéro individu sur 15 points d’écoute, ce qui a permis aux utilisateurs d’identifier les zones où ils ont le plus de chances d’observer l’espèce.

Pour évaluer la qualité et la fiabilité de la carte de chaque espèce, nous avons comparé les abondances relatives prévues aux décomptes observés ainsi qu’aux tendances de détection/non-détection dans les données brutes. Cette évaluation a inclus un examen qualitatif par des ornithologues experts connaissant bien les tendances de répartition spécifiques à chaque espèce dans la province, ainsi que des vérifications visuelles visant à détecter tout surajustement ou extrapolation peu plausible.

Pour certaines espèces, plusieurs itérations de raffinement du modèle et de filtrage des données ont été nécessaires afin de garantir la cohérence des cartes de prévision avec les tendances connues en matière d’occurrence. Dans certains cas, cela a nécessité d’ajuster la structure du modèle et de supprimer les données aberrantes. Seules les cartes de prévision ayant satisfait à la fois à l’examen basé sur les données et à celui basé sur l’avis d’experts ont été jugées aptes à être publiées dans l’atlas.

Références

Bachl, F. E., Lindgren, F., Borchers, D. L., & Illian, J. B. (2019). inlabru: an R package for Bayesian spatial modelling from ecological survey data. Methods in Ecology and Evolution, 10(6), 760-766.

Edwards, B. P., Smith, A. C., Docherty, T. D., Gahbauer, M. A., Gillespie, C. R., Grinde, A. R., … & Zlonis, E. J. (2023). Point count offsets for estimating population sizes of north American landbirds. Ibis, 165(2), 482-503.

Fink, D., T. Auer, A. Johnston, M. Strimas-Mackey, S. Ligocki, O. Robinson, W. Hochachka, L. Jaromczyk, C. Crowley, K. Dunham, A. Stillman, C. Davis, M. Stokowski, P. Sharma, V. Pantoja, D. Burgin, P. Crowe, M. Bell, S. Ray, I. Davies, V. Ruiz-Gutierrez, C. Wood, A. Rodewald. 2023. eBird Status and Trends, Data Version: 2023; Released: 2024. Cornell Lab of Ornithology, Ithaca, New York. https://doi.org/10.2173/WZTW8903

Fick, S.E. and R.J. Hijmans, 2017. WorldClim 2: new 1km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37 (12): 4302-4315.

Guindon, L., Manka, F., Correia, D. L., Villemaire, P., Smiley, B., Bernier, P., … & Boulanger, Y. (2024). A new approach for Spatializing the Canadian National Forest Inventory (SCANFI) using Landsat dense time series. Canadian Journal of Forest Research, 54(7), 793-815.

Hollister, J.W. (2023). elevatr: Access Elevation Data from Various APIs. R package version 0.99.0. https://CRAN.R-project.org/package=elevatr/

Pekel, J. F., Cottam, A., Gorelick, N., & Belward, A. S. (2016). High-resolution mapping of global surface water and its long-term changes. Nature, 540(7633), 418-422. Downloaded from: https://global-surface-water.appspot.com/download

Ressources naturelles Canada.  Couverture des terres du Canada 2020.  Téléchargé depuis : https://open.canada.ca/data/en/dataset/ee1https://ouvert.canada.ca/data/fr/dataset/ee1580ab-a23d-4f86-a09b-79763677eb47580ab-a23d-4f86-a09b-79763677eb47

Sólymos, P., Matsuoka, S. M., Bayne, E. M., Lele, S. R., Fontaine, P., Cumming, S. G., … & Song, S. J. (2013). Calibrating indices of avian density from non‐standardized survey data: making the most of a messy situation. Methods in Ecology and Evolution, 4(11), 1047-1058.

Wang T., Hamann A., Spittlehouse D., & Carroll C. (2016) Locally Downscaled and Spatially Customizable Climate Data for Historical and Future Periods for North America. PLoS ONE 11, e0156720. doi:10.1371/journal.pone.0156720. Downloaded from: https://adaptwest.databasin.org/pages/adaptwest-climatena/