Suite

Nuage de points, corrélation de dépendance et r carré entre deux rasters

Nuage de points, corrélation de dépendance et r carré entre deux rasters


J'ai deux rasters de la même région et de la même résolution et je veux tracer dans un nuage de points leurs valeurs et calculer la ligne la mieux ajustée avec r-carré. Une idée de comment puis-je le faire sur Arcmap ou QGIS ou extraire les valeurs des rasters dans un fichier afin que je puisse l'importer dans Excel? Malheureusement, je n'utilise pas R ni python.


Vous pouvez utiliser le Corrélation d'images, Régression d'image, et Tracé de l'espace des entités dans les outils d'analyse géospatiale Whitebox pour y parvenir. Cela fonctionne sur des images entières, avec la mise en garde que lorsque vous avez des tailles d'échantillon très élevées (c'est-à-dire des millions de pixels), même de très petites différences donneront une signification statistique, qui peut ne pas être physiquement significative.

Voici un exemple de sortie de régression d'image pour deux bandes d'images Landsat 8 :


Faire une régression entièrement dans un SIG est un peu un tour de force, mais il est possible (et a peut-être un certain mérite lorsqu'il est incorporé dans un flux de travail plus long qui serait autrement interrompu par l'utilisation d'une plate-forme de calcul statistique.)

Chaque fois que vous disposez d'un vaste ensemble de données appariées (X, Y), telles que celles fournies par deux rasters colocalisés, alors l'ajustement par les moindres carrés ordinaire est mieux trouvé avec une série d'étapes simples (qui évitent les erreurs en virgule flottante pouvant survenir lors de l'application de formules conventionnelles) :

  1. Trouvez les moyennes de X et Y. Appelez-les respectivement mX et mY. (Ordinairement, ceux-ci seraient représentés comme des rasters ayant des valeurs constantes dans leurs cellules.)

  2. Plus récent X et Y en soustrayant leurs moyennes de chacun. (C'est-à-dire, changez X en X-mX et Y en Y-mY).

  3. Trouvez les moyennes des au carré valeurs recentrées : ce sont les écarts des rasters. (Encore une fois, ceux-ci seraient représentés comme des rasters constants.)

  4. Divisez les valeurs recentrées par les racines carrées de leurs variances. (Appelez ces racines carrées sX et sY, respectivement.) Les rasters résultants ont été standardisés. La plupart de leurs valeurs se situeront entre -2 et 2, quelques-unes dépassant ces limites (parfois bien au-delà).

  5. Multipliez les deux rasters standardisés. La moyenne de ces produits est le « bêta » de la régression : c'est le coefficient normal. Son carré est R-carré. Étant donné qu'il s'agit de la réponse souhaitée, si vous ne l'utilisez pas pour la prédiction ou un autre traitement, vous préférerez peut-être l'obtenir dans un tableau (sur une ligne) plutôt que de le représenter sous forme de raster.

Incidemment, la pente de la droite de régression de Y par rapport à X est obtenue en reconvertissant les valeurs standardisées vers les valeurs d'origine : cela multipliera Y par sY et X par sX. Ainsi, bêta doit être multiplié par sY et divisé par sX. La pente de la droite de régression de X par rapport à Y est également obtenue en divisant bêta par sY et en la multipliant par sX. (Ces deux pentes ne seront pas égales à moins que bêta ne soit à la valeur extrême de 1 ou -1.)


Seuls deux types d'opérations raster sont impliqués dans tout cela : les statistiques zonales (pour faire la moyenne des valeurs sur un raster) et l'arithmétique locale (pour soustraire des valeurs de chaque cellule, mettre les valeurs au carré, prendre des racines carrées, multiplier les valeurs dans deux grilles, etc. ). Pour effectuer les statistiques zonales, créez d'abord une grille de zones en superposant les rasters X et Y et en définissant les valeurs résultantes sur une constante (telle que 0 ou 1); chaque cellule sans valeur à la fois dans X et Y sera nulle. Cela fait de l'ensemble du raster une zone unique sur laquelle toutes les moyennes seront effectuées.


Pour dessiner le nuage de points, envisagez un échantillonnage à partir des rasters afin de limiter les points à un nombre gérable.


Habituellement, vous aurez suffisamment d'informations pour un tel calcul en utilisant quelques milliers de points, et Excel peut être limité si vous utilisez l'ensemble de données complet.

Ainsi, vous pouvez générer un ensemble de points aléatoires (échantillonnage aléatoire simple ou grille régulière avec une origine aléatoire), extraire les valeurs de pixels pour chaque point, puis exporter la table attributaire vers Excel pour un calcul ultérieur.

Dans ArcGIS, vous aurez besoin de « résille » ou de « générer des points aléatoires » pour la première étape, puis « extraire plusieurs valeurs à un point » (licence d'analyste spatial requise) et enfin « table à table »


Nuage de points, corrélation de dépendance et r-carré entre deux rasters - Systèmes d'Information Géographique

Comme indiqué dans le chapitre précédent (à partir de la page 59), l'interrogation des grilles de données dans les SIG peut vous donner une première indication des relations entre différents facteurs. Pour avoir une idée plus précise il faut faire une analyse statistique sur les données disponibles. L'analyse de régression est un outil statistique simple pour examiner les corrélations entre deux ou plusieurs types de données et est la technique statistique la plus couramment utilisée en biologie halieutique. D'abord un rappel des mathématiques de base de l'analyse de régression.

15.1. Régression linéaire

La régression linéaire est une technique pour quantifier la relation, qui peut être vue dans un graphique fait entre deux variables. Par exemple, la figure 15.1 présente la relation entre le nombre de pêcheurs et le nombre de filets maillants dans les différents villages autour du lac Kadim. Cela vous montre une relation quand le nombre de pêcheurs augmente, le nombre de filets maillants dans les villages augmente aussi, c'est ce qu'on appelle une relation positive.

FIGURE 15.1
La relation entre le nombre de pêcheurs et le nombre
de filets maillants dans les différents villages autour du lac Kadim

FIGURE 15.2
Relation entre le nombre de pêcheurs et leurs prises annuelles de carpes

Cependant, si vous regardez les prises annuelles de carpes des pêcheurs individuels (Figure 15.2), vous voyez une relation négative où les prises diminuent à mesure que le nombre de pêcheurs augmente.

Il est bon de savoir qu'il existe une relation linéaire positive entre le nombre de pêcheurs et le nombre de filets maillants au lac Kadim, mais notre objectif est de décrire cette relation par un modèle mathématique ou une équation :

Dans cette équation y , le nombre de filets maillants, est la variable sur l'axe vertical du graphique ou la variable dépendante, tandis que x , le nombre de pêcheurs, représente la variable sur l'axe horizontal ou la variable indépendante. La valeur a (qui peut être négative, positive ou nulle) est appelée l'interception, tandis que la valeur b (qui peut être positive ou négative) est appelée ‘pente’ ou ‘coefficient de régression’. La question est de savoir comment calculer les valeurs de a et b . Vous ne vous soucierez pas des détails, mais dans tous les manuels de statistiques, vous verrez que a et b peuvent être calculés avec les équations suivantes :

x et y sont les valeurs des différentes paires x et y, n est le nombre de paires, est la valeur moyenne de y et est la valeur moyenne de x.

Jetez un œil aux données suivantes présentées dans le tableau 15.1.

TABLEAU 15.1
Données pour le calcul de la régression entre le nombre de pêcheurs et le nombre de filets maillants au lac Kadim

Certains calculs de base des moyennes et des sommes sont fournis au bas du tableau. En utilisant ces valeurs et les formules ci-dessus, nous pouvons calculer notre pente de régression et les paramètres d'interception :

et a = 479-9,7412 * 62 = -127.321

À partir de là, nous pouvons décrire la relation y = 127,3 + 9,74x ou en mots :

Nombre de filets maillants = -127,3 + 9,74 * Nombre de pêcheurs.

Notez que la pente de la droite (9.74) est un nombre positif, indiquant qu'il s'agit d'une relation positive. Cela concorde avec notre interprétation visuelle de la figure 15.1.

Dans le tableau 15.2, les données pour le nombre de pêcheurs et leurs captures annuelles de carpes (PUE) sont fournies. Calculer la relation de régression

TABLEAU 15.2
Relation entre le nombre de pêcheurs et la CPUE des carpes du lac Kadim

Résultat de l'analyse de régression : y = 1 465 - 8,6627x.

De nos jours, les analyses de régression sont devenues plus faciles car elles sont incluses dans tous les tableurs tels que Lotus 1-2-3 et Microsoft Excel. Dans Microsoft Excel, l'analyse de régression est effectuée dans des graphiques créés avec les ensembles de données.

Faisons notre exemple de Lake Kadim dans Microsoft Excel :

1. Démarrez Microsoft Excel, ouvrez la feuille de calcul ‘Lake Kadim regression analysis.xls’, à partir du dossier 㢳_Lake_Kad_regr’. Vous voyez l'ensemble de données avec le nombre de pêcheurs, leur CPUE et deux graphiques.

2. Activez un graphique en cliquant dessus.

3. Accédez à Graphique/Ajouter une ligne de tendance via la barre de menu. La fenêtre Ajouter une ligne de tendance s'ouvrira (Figure 15.3) et vous sélectionnerez linéaire en cochant sa case. Cliquez ensuite sur l'onglet Options dans la fenêtre Ajouter une ligne de tendance et cochez Afficher l'équation sur le graphique et Afficher la valeur R au carré sur le graphique (Figure 15.4). Cliquez sur OK.

FIGURE 15.3
La fenêtre Ajouter une ligne de tendance

FIGURE 15.4
Vérification de l'affichage de l'équation sur le graphique et de l'affichage du R au carré sur le graphique

Le graphique du graphique s'affiche à nouveau. Maintenant avec une ligne droite (la régression calculée) et la relation y = -8,6627x + 1465,8, que vous avez calculée précédemment affichée dans le graphique (Figure 15.5).

FIGURE 15.5
La régression Microsoft Excel entre le nombre de pêcheurs et leur CPUE dans le lac Kadim

Dans le graphique, vous voyez également une autre valeur : R 2 = 0,9164. R au carré, ou le coefficient de détermination, est le carré du coefficient de corrélation R. C'est une mesure de l'association linéaire entre deux ensembles de données et reflète la quantité de variation de la variable dépendante qui peut être expliquée par la variation de la variable indépendante. variable. La valeur des valeurs R-carré est comprise entre 0 (reflétant absolument aucune relation linéaire entre les variables) et 1 (indiquant une corrélation parfaite). Par exemple, si R au carré = 0,25, nous pourrions dire que la variance de la variable indépendante explique 25 % de la variation de la variable dépendante. Plus la valeur 1 est proche, plus la corrélation entre les deux variables sera élevée. R-carré est calculé comme (Figure 15.6) :

FIGURE 15.6
L'équation pour calculer le R au carré

Cependant, une valeur élevée de R ne signifie pas que la droite de régression est toujours statistiquement valide. La façon officielle de voir cela est d'effectuer un test t sur le coefficient de régression b et de tester la valeur t calculée, ou d'effectuer une ANOVA (ou analyse de la variance) et de tester la valeur calculée de la statistique F . Cela peut être fait dans n'importe quel logiciel statistique.

15.2. Régression avec un script Avenue dans ArcView

Contrairement à Microsoft Excel, l'analyse de régression ne fait pas partie intégrante d'ArcView. Heureusement, ArcView fournit un langage de codage interne qui nous permet d'écrire ce type de fonctions. De nombreuses fonctions personnalisées, y compris l'analyse de régression, ont déjà été écrites par les utilisateurs d'ArcView et sont disponibles pour un usage public (par exemple, voir ArcScripts sur http://gis.esri.com/arcscripts/index.cfm).

Un exemple de script de régression est inclus avec les données sur le CD qui permet une analyse de régression linéaire. La méthode disponible est plutôt basique, ne peut être appliquée que sur des fichiers de formes et a une fonction de nuage de points limitée. Cette application vous sera présentée en utilisant l'exemple du lac Kadim et en réalisant une analyse de régression entre le nombre de pêcheurs et le nombre de filets maillants dans les villages autour du lac Kadim.

1. Démarrez ArcView, ouvrez un nouveau projet et une nouvelle vue. Ajoutez à la vue les thèmes (à partir du dossier 㢳_Lake_Kad_regr’): ‘Pais pesca country.shp’, ‘Lake kadim bound.shp’, ‘Lake Kadim data.shp’ et ‘Fishing village lac kadim.shp’.

2. Vérifiez la projection et le répertoire de travail.

3. Vous devez d'abord ajouter le script. Fermez la vue et ouvrez un nouveau script dans la fenêtre Projet (Figure 15.7).

FIGURE 15.7
Ouvrir un nouveau script

4. Vous arrivez dans la fenêtre Script, où vous devez ouvrir le script avenue ‘bvreg.ave’ Il s'agit d'un fichier texte et pour votre commodité il est placé dans le même dossier que le fichier Thème, 㢳_Lake_Kad_regr’ . Allez à S cript/Load Text File. (Figure 15.8).

5. La fenêtre de chargement du script apparaît. Accédez au sous-répertoire approprié et sélectionnez le script ‘bvreg.ave’ et cliquez sur OK (Figure 15.9).

FIGURE 15.8
Chargement d'un fichier texte de script

FIGURE 15.9
Sélection du script de régression

Pour exécuter le script, vous devez connaître quelques astuces. Vous devez d'abord compiler le script. ‘Compilation’ signifie simplement qu'ArcView recherche les erreurs dans le texte du code, puis convertit le code dans un format qu'ArcView peut exécuter directement. Après avoir cliqué sur le ‘bouton Compiler’ (Figure 15.10), vous verrez le bouton à côté avec la personne qui court devient actif. C'est le bouton pour exécuter le script. Si vous cliquez sur ce bouton, vous obtiendrez un message d'erreur : ‘A(n) le projet ne reconnaît pas la demande GetActiveThemes’ (Figure 15.11).

FIGURE 15.10
Compiler un script

FIGURE 15.11
Message d'erreur de régression

La raison principale pour laquelle nous recevons ce message d'erreur est qu'ArcView ne sait pas où chercher les données pour effectuer l'analyse. Le moyen le plus simple de le dire est d'ouvrir la vue, de trouver le thème qui contient les données pour l'analyse de régression et de rendre ce thème actif. La façon de faire cela semble compliquée, mais ne vous inquiétez pas, cela fonctionne et une fois que vous l'avez fait plusieurs fois vous connaissez l'astuce. Fondamentalement, vous devez mosaïquer la fenêtre pour que la fenêtre Vue et la fenêtre Script soient visibles. Ensuite, vous pouvez passer de l'un à l'autre.

6. Allez dans la fenêtre de script vers W indow/T ile via la barre de menu (Figure 15.12). Vous obtiendrez deux fenêtres sur votre écran. Ouvrez la vue (Figure 15.13).

FIGURE 15.12
Mosaïque dans la vue de script

FIGURE 15.13
Script en mosaïque et fenêtre d'affichage

7. Après avoir ouvert votre vue, placez à nouveau les fenêtres en mosaïque via W indow/T ile dans la barre de menu. Vous devriez maintenant voir trois fenêtres ouvertes dans votre projet, et votre fenêtre d'affichage devrait être votre fenêtre active. Vous pouvez le voir parce que la barre de fenêtre de cette fenêtre est bleue. Votre vue doit avoir quatre thèmes répertoriés. Cliquez sur les mots ‘Fishing village lake kadim.shp’ (pas sur la case à cocher) pour activer le thème. ArcView saura maintenant quel thème contient les données pour l'analyse de régression (Figure 15.14).

8. Revenez à la fenêtre de script en cliquant quelque part à l'intérieur, puis cliquez sur le bouton Exécuter le script (Figure 15.15).

FIGURE 15.14
Activation du thème de données pour la régression

FIGURE 15.15
Exécuter le script

9. La fenêtre de régression bivariée apparaîtra. Vous devez d'abord indiquer la variable indépendante (X): ‘Pêcheurs’ (Figure 15.16), après avoir sélectionné ‘Pêcheurs’, cliquez sur OK. Dans la fenêtre suivante, vous indiquez la variable dépendante (Y), ‘Gill_nets’ (Figure 15.17), et cliquez sur OK .

FIGURE 15.16
Sélection de la variable indépendante

FIGURE 15.17
Sélection de la variable dépendante

10. La fenêtre des résultats de la régression bivariée apparaît et fournit les résultats de la régression. Dans ce cas, les résultats sont similaires au calcul que nous avons effectué auparavant, cliquez sur OK. Dans la fenêtre suivante, il vous sera demandé si vous souhaitez créer un nuage de points. Cliquez sur Oui et le nuage de points avec la régression calculée apparaît (Figure 15.19).

FIGURE 15.18
Résultats de la régression bivariée

FIGURE 15.19
Le nuage de points de la régression calculée

15.2.1. Analyse de régression des données du lac Kadim à l'aide d'un script avenue

Faire l'analyse de régression entre le nombre de pêcheurs et la CPUE avec les données du thème ‘Villages de pêcheurs du lac Kadim’ et les comparer avec les résultats de l'analyse réalisée sous Excel.

Faites un certain nombre d'analyses de régression avec les données du lac thématique Kadim et remplissez le tableau 15.3.

TABLEAU 15.3
Résultats de l'analyse de régression des données brutes du lac Kadim


Vous pouvez facilement trouver une réponse à cette question en utilisant simplement Google ou un manuel de statistiques de base. Par exemple, j'ai basé cette réponse sur un chapitre d'un livre intitulé Fundamental Statistics for the Behavioral Sciences de Howell (2013).

R2, le coefficient de corrélation au carré, explique la force de la relation entre les deux variables dans votre nuage de points. Supposons que vous ayez deux variables, X (prédicteur) et Y (résultat), il y a beaucoup de variabilité dans Y. Une partie de cette variabilité sera liée à votre variable prédictive, X, mais certaines seront également du bruit, également appelé Erreur.

Si X (disons la taille) est un bon prédicteur de Y (disons le poids), alors une grande partie de la variabilité de la taille sera associée à la variabilité du poids. Cela signifie essentiellement qu'une partie de la raison pour laquelle les gens varient en poids sera parce que les gens diffèrent dans leur taille. Donc, une partie des raisons pour lesquelles les gens diffèrent par leur poids est qu'elles diffèrent par leur taille.

R2 = variation de Y (dans notre exemple poids) expliquée par X (dans notre exemple taille) / Variation de Y (poids)

Compte tenu de l'équation ci-dessus, R2 est égal au pourcentage de la variabilité du poids (Y), que la taille (X) est capable de prédire ou d'expliquer. Dans votre cas, la valeur R2 signifie que votre prédicteur explique moins de 1 % de la variabilité de votre variable de résultat.


10.5 MAUP

10.5.1 Obtenir les données

En utilisant une mise à jour récente, nous pourrions également changer cela un peu avec le package pins() qui télécharge et met en cache l'URL et vérifie automatiquement les modifications de fichier, en ne téléchargeant à nouveau que si nécessaire !

pins() vous permet également de partager facilement des données, consultez le billet de blog RStudio pins pour plus d'informations.

  1. Prenez les données téléchargées et filtrez-les en fonction du nom de fichier qui contient : Borough OU Ward_ ET .shp en utilisant grepl()
  1. Maintenant, lisez les deux fichiers en utilisant map() qui applique une fonction (ici st_read() du package sf) à une liste. map() vient du package purrr qui étend la capacité de programmation fonctionnelle de R.

Pour mapper ou accéder à chaque fichier de formes individuel, il suffit de…

  1. Et pour OSM, nous le téléchargerons depuis geofabrik, vous pouvez également utiliser l'interface de programmation d'applications (API) OSM, mais le nombre de points téléchargeables par appel est limité, vous devrez donc faire quelque chose d'un peu plus compliqué pour obtenir l'ensemble de Londres… cependant, j'ai fourni un exemple de l'appel api.

10.5.2 Projeter les données

10.5.3 Fonctions

  1. Souvenez-vous maintenant de notre fonction que nous avons créée pour joindre nos données Airbnb (ou hôtels) à la couche des arrondissements de Londres… faisons à nouveau la fonction…

10.5.4 Boucles

  1. D'accord, mais nous voulons obtenir le nombre de points Airbnb par quartier et arrondissement de Londres… comment pouvons-nous faire cela ?… enfin manuellement bien sûr… comme ça…

Mais nous pouvons également automatiser cela en utilisant une boucle (une boucle while ou for). J'ai utilisé une boucle while ici, car lorsque j'ai fait mon MSc, vous n'étiez pas en mesure de mettre une boucle for dans une boucle for. Je crois que cela a maintenant changé, mais à cause de cela, un jour, j'ai dû tout changer, j'utilise toujours par défaut une boucle while.

Dites-nous ce qu'est déjà une boucle ?

Une boucle vous permet de parcourir quelque chose en ajoutant 1 (ou n'importe quelle valeur) à chaque fois… par exemple, regardons une boucle de base. Vous devez tout exécuter dans la boucle à la fois de while au >. Si vous ne créez qu'un Rscript normal, vous pouvez définir des points d'arrêt - le code s'arrêtera chaque fois qu'il atteint le point d'arrêt dans la boucle. Vous ne pouvez pas le faire pour le moment avec des morceaux de code RMarkdown, je développe normalement la boucle en dehors de celle-ci en boucle, puis je la rassemble.

Cette boucle renvoie les valeurs de 1 à 5, car nous avons commencé avec une valeur de 1, puis ajouté 1 pour faire 2. Il est resté en dessous de 6, donc le code s'est exécuté à nouveau en imprimant le 2 puis en a ajouté 1 pour faire 3 et ainsi de suite. Comme nous l'avons spécifié moins de 6, cela s'est arrêté là.

Nous pouvons également enregistrer ces résultats dans différentes variables, mais nous devons commencer par créer une liste (ou une base de données / tout ce dont vous avez besoin) pour les enregistrer dans

Ici, nous utilisons la variable basicloop pour indexer notre liste vide.. donc chaque fois que nous ajoutons 1, cela change la valeur de l'index….voyez ce que je veux dire…

Bon, alors comment allons-nous appliquer cela à nos données. Nous avons deux fichiers .shp (arrondissements et quartiers) dans une liste auxquels nous voulons appliquer notre fonction.

Tout d'abord, configurons la longueur à laquelle s'arrêter, faisons une liste vide et un point de départ. Comme nos données sont dans une liste, nous voulons juste la longueur de celle-ci (sous forme de nombre).

Eh bien, c'est une question difficile… généralement, les boucles étaient considérées comme inefficaces dans R, mais je n'ai pas encore trouvé de méthode qui vous permettra d'incrémenter une variable tout en gardant une autre cohérente… je veux dire par là comme ce que nous avons fait dans notre boucle. Nous avons modifié les données spatiales jointes aux données Airbnb — les données Airbnb sont restées les mêmes. Il existe quelques autres fonctions comme mapply() , sapply() et tapply() qui fonctionnent dans le même format que lapply() mais à ma connaissance, elles incrémenteront toujours toutes les variables en même temps.

En termes de clarté du code, je vais me référer à la section 21.5 sur l'itération dans Wickham $ Grolemund (2017)…

« Certaines personnes vous diront d'éviter les boucles for car elles sont lentes. Ils ont tort ! (Eh bien au moins, ils sont plutôt obsolètes, car les boucles n'ont pas été lentes depuis de nombreuses années). Le principal avantage de l'utilisation des « autres fonctions » n'est pas la vitesse, mais la clarté : elles rendent votre code plus facile à écrire et à lire. »

Cela dit, dans ce cours, vous n'êtes pas marqué sur la qualité ou l'efficacité de votre code. Bien sûr, je veux que vous écriviez du bon code, mais si vous écrivez vraiment du code pour votre mission ou tout projet futur, mon conseil est le même. Obtenez quelque chose qui fonctionne et répond au schéma de notation (lisez le schéma de notation !) et puis si vous avez le temps, améliorez-le plus tard. Il y a aussi une section plus loin sur l'écriture de code avancé dans le devoir.

10.5.5 Mappage avancé (encore)

  1. D'accord, nous pouvons donc en quelque sorte voir la différence entre les niveaux spatiaux (arrondissements et quartiers) mais regardons de plus près à l'intérieur de Westminster… voici le «préambule» de la carte du dépliant…
  1. Maintenant, cartographions-le en utilisant ce que nous venons de spécifier… j'ai ajouté quelques fonctionnalités de plus que dans la carte, ce qui est pratique

Regardez autour de la carte… L'arrondissement de Westminster utilise une échelle tenant compte de toutes les autres valeurs d'arrondissement de Londres, tandis que l'échelle de quartier est spécifique à Westminster. Utilisez le code suivant pour explorer les valeurs…


Les coordonnées géographiques

Les coordonnées géographiques des échantillons constituent une source d'informations inestimable, allant de l'affichage de leur distribution spatiale à la récupération des variables environnementales. Chaque fois que vous travaillez sur le terrain, l'utilisation d'un GPS est le meilleur moyen d'enregistrer les coordonnées des échantillons. En tant que tel, nous recommandons fortement d'enregistrer l'emplacement de chaque échantillon, au lieu de l'emplacement du centroïde d'une population par exemple. Premièrement, il permet une récupération plus précise des valeurs environnementales. Deuxièmement, attribuer le même emplacement à plusieurs échantillons invoque une pseudo-réplication, un biais statistique qui doit être traité dans une analyse plus approfondie. Troisièmement, les coordonnées des individus proches permettent une mesure correcte de la dispersion, en utilisant par exemple une relation génétique par paire avec la distance. Concernant les appareils GPS, les GPS standards, et dans une moindre mesure les smartphones, sont suffisamment précis dans la plupart des cas. Cependant, des appareils plus précis, tels que le DGPS (GPS différentiel), sont recommandés pour les études à l'échelle locale dans lesquelles les échantillons sont situés à moins de quelques mètres l'un de l'autre : la précision de l'emplacement doit rester dans la résolution spatiale du grain.

Lorsque les coordonnées GPS ne sont pas enregistrées, il est toujours possible d'approcher les emplacements des échantillons à l'aide d'images satellites ou en encodant l'adresse de l'emplacement (géo-référencement ou géocodage), mais avec une précision moindre. Dans le premier cas, la création d'une nouvelle couche vectorielle superposée sur une image satellite ou une carte en ligne (voir Section suivante) permet de récupérer les coordonnées des échantillons à partir d'un emplacement approximativement connu (par exemple, un carrefour, un arbre, une rivière Docs.QGIS, 2014). Pour ce dernier cas, des plugins ont été développés pour lire un fichier délimité par du texte contenant des adresses (par exemple, votre propre adresse personnelle) que vous souhaitez localiser (par exemple le plugin MMQGIS dans QGIS, Mangomap, 2012 MMQGIS Plugin, 2012). Il faut noter que chaque ligne doit contenir l'adresse, la ville, l'état et le pays.

Une autre considération essentielle est le choix du système de référence de coordonnées approprié. En effet, les appareils GPS affichent les coordonnées d'un point en valeurs de latitude et de longitude, généralement dans le système géodésique mondial (WGS84). Il s'agit d'un système de référence global dans lequel la Terre est représentée par un ellipsoïde, et chaque position à la surface est définie par deux angles au centre de la Terre : la latitude et la longitude. Cependant, les systèmes projetés pour lesquels un emplacement géographique est converti de l'ellipsoïde (distances exprimées en degrés) à un emplacement correspondant sur une surface bidimensionnelle (x et y exprimés en mètres) sont préférés pour les analyses. Il est important de noter que, bien qu'il existe des systèmes mondiaux couvrant l'ensemble de la planète, chaque pays ou région possède son propre système de coordonnées qui est localement plus précis que le système mondial. Là où il n'existe pas de système national de projection, il est toujours possible d'utiliser le système de coordonnées Universal Transverse Mercator (UTM), un système de coordonnées projetées couvrant l'ensemble du globe et le divisant en soixante zones longitudinales de largeur 6 (Dmap, 1993). Même si les logiciels SIG traitent généralement des systèmes de projection différents, la reprojection manuelle de toutes les couches dans le même système de projection local est recommandée pour éviter les incompatibilités potentielles (voir section suivante). Cependant, différents SIG peuvent ne pas utiliser exactement le même nom pour un système de coordonnées. Par conséquent, pour faciliter l'identification des systèmes de coordonnées à travers la diversité des logiciels SIG, la base de données EPSG (European Petroleum Survey Group) (EPSG, 1985) est une base de données largement utilisée référençant tous les systèmes de coordonnées projetés, mis en œuvre dans chaque SIG et leur fournissant un ID unique (Maling, 1992), par exemple, EPSG : 4326 correspond au système de référence WGS84.


Comparaison des jeux de données

La qualité des MNT sous-jacents, en termes de précision verticale et spatiale, influence directement la qualité des couches géomorphométriques nouvellement développées. Dans cette section, nous évaluons la sensibilité des couches géomorphométriques par rapport à la précision DEM et nous la décrivons en trois sous-sections organisées comme suit : i) Évaluation de la projection géographique, où nous montrons les artefacts possibles qui pourraient découler du calcul de la variable géomorphométrique sous Système géodésique WGS84 et projection Equi7 ii) MERIT-DEM versus couches géomorphométriques dérivées de 3DEP-1, où nous décrivons la divergence des couches géomorphométriques les plus courantes obtenues à partir de MERIT-DEM et 3DEP-1 DEM iii) MERIT-DEM vs élévation LiDAR comparaison, qui décrit l'influence de la suppression du biais de hauteur des arbres dans MERIT-DEM en utilisant le DTM et le DSM obtenus à partir de LiDAR. Globalement, ces analyses mettent en évidence la qualité des couches géomorphométriques dérivées de MERIT et permettent d'identifier les erreurs potentielles dans les MNT.

Évaluation de la projection géographique

Des projections cartographiques sont nécessaires pour cartographier la surface de la Terre sur un plan 2D, ce qui est particulièrement pertinent pour les MNT. Quelle que soit la projection utilisée, un certain type de distorsion se produira dans la carte résultante, mais la sélection d'une projection appropriée devrait, en principe, minimiser l'étendue et le type de distorsions 36 . En général, les distorsions de la carte diminuent à mesure que la zone géographique est réduite, par exemple, lors du passage d'une échelle mondiale à une échelle continentale ou régionale. De plus, les distorsions augmentent au fur et à mesure que l'on se déplace le long d'une surface, loin du centre de projection. Cette distorsion est une propriété incontournable des projections cartographiques, et il est important d'évaluer son effet sur tout type d'analyse spatiale, notamment sur celles réalisées à grande échelle. Pour évaluer l'effet des distorsions de la carte sur les variables géomorphométriques, nous avons analysé les variations de pente sous deux emplacements géographiques ayant des distorsions de surface différentes.

La pente est définie comme le taux de changement d'altitude le long de la direction de l'écoulement de l'eau, et calculée en utilisant une fenêtre mobile de 3 × 3 cellules. Le taux de changement peut être exprimé en pourcentage du changement d'altitude sur 100 mètres. Afin d'avoir le même poids dans le X et oui directions sur le taux de changement, la taille de la cellule doit avoir la même dimension dans le X et oui directions. Ce n'est pas le cas lorsque le système de coordonnées géographiques est utilisé et en particulier lorsque le gradient longitudinal (en oui dimension) reste constant par rapport au gradient latitudinal (en X dimension), qui diminue en s'éloignant de l'équateur.

La pente du terrain étant l'une des variables géomorphométriques les plus utilisées, tous les logiciels de SIG et de télédétection ont adopté des algorithmes pour la calculer. Cependant, aucun de ces algorithmes n'utilise une procédure de correction pour tenir compte de la distorsion de la grille dans le X et oui dimension. En d'autres termes, ils traitent les données de latitude et de longitude comme une matrice sur une grille carrée. Au lieu de cela, il est établi que les méridiens convergent vers les pôles et que la distance des cercles parallèles ne varie que légèrement. Ainsi, une grille carrée à l'équateur devient une grille rectangulaire aux latitudes plus élevées. Par conséquent, pour quantifier l'influence des changements dans la dimension latitudinale de X, nous avons comparé les valeurs de pente sous WGS84 avec celles sous Equi7 pour deux zones d'étude de 500 × 500 cellules de grille, avec MERIT-DEM comme couche de base statique. Cette procédure peut être illustrée à l'aide d'un MNT simulé sous deux emplacements distincts, l'un dans la zone subtropicale et l'autre dans la zone subarctique. Néanmoins, pour démontrer un effet potentiel réel, nous sélectionnons une seule zone de MERIT-DEM.

Les graphiques présentés dans la Fig. 2 montrent la différence dans les calculs de pente comme résultat direct de l'utilisation de la latitude et de la longitude sur WGS84 avec celles d'une grille carrée sur Equi7. Pour comparer le même MERIT-DEM sous WGS84 et Equi7 à deux endroits distincts, nous transposons (c'est-à-dire, coordonnées de décalage Equi7 - note : aucune reprojection) la zone subtropicale MERIT-DEM (centre de l'image : longitude -83.26, latitudes 9.05 au Costa Rica) (Fig. 2g) vers une zone subarctique (centre de l'image : longitude -38,19, latitudes 72,80 au Groenland) (Fig. 2a), sous Equi7. Cela produit un simple déplacement le long de l'axe nord-sud sans modification de la valeur du pixel d'élévation (pas d'interpolation). Après avoir reprojeté le MERIT-DEM sur WGS84 (Fig. 2b, h), nous avons calculé la pente (Fig. 2f, d) et l'avons ensuite reprojetée à Equi7 (flèches bleues sur la Fig. 2) pour comparer les résultats à l'aide du scatter. parcelles. La figure 2i correspond à la corrélation de pente de la zone subarctique, tandis que la figure 2j correspond à la zone subtropicale). Dans chaque nuage de points (figures 2i,j), la ligne rouge représente une relation 1:1, tandis que la ligne noire est un modèle de régression ajusté entre les variables. Les variations dans les calculs de pente entre les deux systèmes sont minimes dans la zone subtropicale (Fig. 2j), car la zone d'étude est adjacente à l'équateur. Cependant, dans la zone subarctique, les variations sont significativement différentes, car la pente calculée sous WGS84 est sous-estimée par rapport à la pente Equi7. En effet, les pentes exposées à l'est ou à l'ouest verront leur pente considérablement sous-estimée en raison de l'étirement des X dimensions dans le sens est-ouest (noter tous les points notamment sous la ligne rouge de la Fig. 2i). En revanche, les versants nord et sud peuvent être modérés correctement estimés (noter tous les points proches de la ligne rouge de la Fig. 2i).

Évaluation du biais de projection. Graphical representation of the difference in terrain slope calculation due to the effect of using the World Geodetic System 1984 (WGS84) (raster panels right-hand side) compared to the Equi7 projection (raster panels left-hand side). A study area located in the subtropical zone (image centre: longitude −83.26, latitudes 9.05) was used to subset the MERIT-DEM for an area of 500 × 500 grid cells (g). This area has been transposed to a subarctic zone (image centre: longitude −38.19, latitudes 72.80) under the Equi7 projection (une). After having been reprojected to WGS84 (b,h), the variable slope was calculated in the four conditions (c–f), and then reprojected back to Equi7 for comparisons (see blue line-arrows). The scatter plots on the right-hand side shows the WGS84-MERIT slope (d,f) vs. the MERIT slope under the Equi7 projection (c,e), respectively for the subarctic zone (i) and for the subtropical zone (j). The red lines represent the 1:1 relationship and black lines represent a linear model between the two axes. The slope calculated under WGS84 in the subarctic zone is clearly underestimated compared to the one calculated under the Equi7 projection.

Similar to slope, all geomorphometric variables are influenced by underlying grid distortions. In particular, the slope is influenced by both length and angular distortions, as are all of the other geomorphometric variables listed under the first and second derivatives group. In contrast, the ruggedness geomorphometric variables are influenced more by areal distortions because of elevation differences at the pixel level. These results emphasise the importance of computing the geomorphometric variables under the Equi7 projection.

In conclusion, it is not that the WGS84 geodetic datum is wrong and distorted but its treatment of latitudinal and longitudinal grids as squares is erroneous, as in the Plate Carrée projection. Consequently, the calculation of any geomorphometric variables under WGS84 should be avoided.

3DEP-1 versus MERIT-DEM comparison

For geomorphometrical and hydrographical applications, the elevation difference between two DEMs is important since any application will be contingent on the values of the derived geomorphometric variables, for example, impacting on the delineation of streams and catchments. In the following sections, we analyse the difference in the elevation values between the 3DEP-1 and MERIT DEMs, as well as their derived geomorphometric variables. The 3DEP-1 is a LiDAR-based DEM and given its high accuracy, can be used as a reference elevation that has negligible errors. Initially, we compare the elevation difference between 3DEP-1 and the MERIT-DEM, and subsequently we analyse how the differences in DEMs influence the derived geomorphometric variables.

Comparing DEMs using the Elevation Deviation Index (EDI)

The elevation difference, or deviation, at pixel level between two DEMs can be expressed as

where X et y are the elevation values in each single pixel i. The (>_) value is equal to 0 if the two DEMs have the same elevation. The overall raster of (>_) values represents the Δ surface.

To identify areas where the deviation is stronger, the deviation at each pixel needs to be considered with the surrounding elevation pixel values (Xi+1). In our case, we label 3DEP-1 as y and MERIT-DEM as X. Therefore, considering a circular window of 23 × 23 pixels that slides across y, it is possible to obtain the standard deviation, which estimates the local elevation roughness. Mathematically, the standard deviation of yi in a moving circular window is expressed as:

If we integrate the (>_) and its surrounding standard deviation, we obtain the Elevation Deviation Index (EDI), which is defined as the ratio

EDI represents the relative deviation over the surrounding elevation variability in the moving window. The component k, in above equation, is used to prevent the situation that areas completely flat, with σi = 0, will produce infinite values of the EDI. In our case, we set k = 0.1, which is a very small value compared to the calculated σi, even in quasi-flat areas. k does not influence the σi and consequently the overall performance of the EDI.

For example, a local elevation difference of 1 m will create a higher index in flat areas compared to mountain regions and the index can be positive or negative with respect to the (>_) . The EDI can be used to select zones where the elevation difference between the 3DEP-1 and the MERIT-DEM is substantial considering the roughness of the surrounding areas. Hence, flat areas will be more sensitive to the EDI compared to steep, mountainous areas. We expect that areas with extreme EDI will be more prone to deviating stream networks compared to zones with EDI close to 0 (see Fig. 3).

Elevation Deviation Index. Elevation Deviation Index EDI (c) obtained as the ratio of elevation difference (une: 3DEP-1 - MERIT) and elevation standard deviation calculated using a moving window of 5 × 5 pixels of the 3DEP-1 (b). The coordinates reported in une are in Equi7 and expressed in metres. The study area refers to a zone of 18.4 × 20 km, which has a high level of forest cover, and is located in Alberta, Canada, close to Jasper National Park - image centre 118.25°W 53.29°N. The same area is used in the Fig. 4 to assess the geomorphometric variables.

When comparing DEMs with unknown or significant errors, the standard deviation of the mean of the two DEMs can be calculated. Besides, this standard deviation is a measure of roughness, and the window size (Eq. 2) reflects neighbouring influences. A large window size will produce a larger standard deviation and thus lower EDI, on average. The moving window size can be adjusted with respect to the resolution of the DEMs or on the basis of the surrounding roughness. The EDI can be applied on a global scale by comparing different DEMs, and highlighting areas where the DEMs have discrepancies.

Figure 3 shows EDI and its components for an area of 18.4 × 20 km. The extreme Δi values (black and blue areas in une) do not necessarily produce extreme EDI. With respect to the EDI, in the largest part of the study area, the DEMs are in agreement (yellow - green colour) and located in zones with a high level of roughness. On the contrary, in flat areas (blue colour in b) the EDI can reach extreme values (blue and dark red colour in c).

Comparing the continuous geomorphometric variables

To compare the geomorphometric variables derived from the 3DEP-1 and MERIT-DEM under the same scale unit, we normalise the difference expressed as a Δ surface. Hence, we deal with the difference (for example pcurv-3DEP-1 - pcurv-MERIT) by scaling all positive values to fall between 0 and +1, and negative values to fall between -1 and 0. As a result, the difference value at 0 remains at 0 when scaled (e.g. 0, 9 scaled to 0, 1 −3, 0 scaled to −1, 0).

Consequently, the normalised Δ surface derivative from each geomorphometric variable can be compared having the same unit and can be used to assess the sensitivity of the variables to the differences in DEM elevation. In fact, in instances where the Δ surface has a value close to 0, this suggests that a geomorphometric variable is not strongly influenced by the difference between the two DEMs. In contrast, in instances where the Δ surface has several pixels with negative or positive values, this means that they are influenced by the DEM’s difference.

Figure 4 shows an overview of the normalised Δ surface for each geomorphometric variables. Two elevation plots (Fig. 4a,b) show the 3DEP-1 and MERIT DEMs and relative scatter plot (Fig. 4c). The elevation difference (Fig. 4d) shows values ranging from −113 m to +216 m. The largest values of difference are located close to the peak areas, and the smallest values are concentrated in the valley areas. Figure 4e shows the normalised version. In contrast, the other plots show the spatial variability of the geomorphometric difference, expressed with normalised values. Values close to −1 and +1 mean high sensitivity to elevation difference, and conversely, values close to 0 mean less sensitivity. In general, the overall correlation between 3DEP-1 and MERIT-DEM is very high, with the blue line representing a fitted regression model, which is very close to the 1:1 red line. These results are in line with other studies that evaluate the accuracy of the MERIT-DEM 17,19 .

Normalised difference maps. Normalised difference maps represented as Δ surface, for each geomorphometric variable derived from 3DEP-1 minus MERIT-DEM. To compare the geomorphometric variables under the same scale unit, the difference has been scaled from -1 to 1 (minimum and maximum stretching) keeping the 0 value (no difference) at the 0 position. The bottom plot reports the normalised difference as mean (blue line) and standard deviation values (orange vertical lines) of the maps. The mean and standard deviation plot helps to identify which geomorphometric variables are more sensitive to variation in the DEMs (e.g. high sensitivity for variables derived from slope and aspect). The coordinates u in Equi7 are expressed in metres and refer to a study area of 18.4 × 20 km, which has a high level of forest cover, and is located in Alberta, Canada, close to Jasper National Park - image centre 118.25°W 53.29°N.

The differences in elevation are mainly due to the radar beam’s shadow, which is usually evident in steep terrain. In fact, in Fig. 4d, the area with high elevation difference is located in the south-west corner with Δ values larger than 200 m. The greatest relative deviations are in high-relief areas for most of the geomorphometric variables except for the compound topographic index, convergence, and sine and cosine of aspect (see Fig. 4j,x,l,m). For these exceptions, the relative deviations are greatest in low-relief areas, especially in the valleys. Indeed, flat areas are very sensitive to the DEMs accuracy and slight variations in the elevation can switch the aspect to the opposite direction.

The behaviour between the compound topographic index (cti) and stream power index (spi) differs due to the logarithmic scale used in the cti. Consequently, the deviation of the cti is visible when there is a small variation of the flow accumulation. On the other hand, the spi does not employ a logarithmic scale, and the deviation is only evident when there is a drastic change in the flow accumulation areas that are adjacent to stream locations.

Visually, the aspect-sine, the aspect-cosine, slope, eastness, northness and the convergence are geomorphometric features that are very sensitive to differences between the two DEMs (see Fig. 4f,g,h,x).

To support the visual assessment of the maps in Fig. 4 with numerical values, we analysed the normalised Δ surface by plotting the mean values and standard deviation for the positive and negative values (see Fig. 4). In fact, similar patterns can also be seen for any of the aforementioned variables with high standard deviation (see vertical lines). The aspect is very sensitive to DEM differences in both steep terrain and flat areas. Even the derived sine and cosine Δ surface show black and blue areas (+1 and −1, respectively) in the central valley (see Fig. 4i,m). Note that these areas are not apparent in the other variables. Just as 1st partial derivatives have been used to detect artefacts in DEMs 19 , the Δ surface of the aspect-sine and aspect-cosine can be used to highlight areas where the two DEMs show differences in elevation.

Comparing the categorical geomorphometric variables

To assess the sensitivity of pattern delineation of the geomorphological forms derived from MERIT-DEM and 3DEP, we compare the geomorphological classification agreement for an area of 300 × 300 km (3000 × 3000 picels at 100 m spatial resolution) in South Dakota, USA. Figure 5 shows two raster plots (a,b) with a similar pattern at large scale but when observed in finer detail, there are differences in the classification at the pixel level (see Fig. 5a,b magnified circle). A common way to analyse the differences between two classifications is the calculation of a so-called confusion (or error) matrix 62 . The confusion matrix displays the probabilities with which pixels belonging to a certain class in one product appear in the same or a different class in the compared product. A confusion matrix can therefore be used to illustrate not only the degree to which the two classifications agree but also reveal how likely a class is misclassified.

Geomorphological forms maps and confusion matrices. The geomorphological forms have been computed for a study area of 3000 × 3000, 100 m pixels in South Dakota (USA) derived from MERIT (a–c) and 3DEP-1 (b–d), respectively. The confusion matrix values are expressed in percentages of the MERIT-DEM classes, with the sum of vertical values equal to 100 (d) and of the 3DEP-1 classes, with the sum of the horizontal values equal to 100 (c). The sum of the values in the blue boxes is equal to 100, and so on for each row (c) and column (d).

In order to allow a numerical comparison of the geomorphological classifications, we calculate two confusion matrices among the 10 classes in each product (see Fig. 5). One is expressed as percentage of MERIT-DEM classes such that the sum within each row is equal to 100 (see Fig. 5c). Considering 3DEP-1 as a reference product, this would give the “user accuracy” of the MERIT-DEM geomorphometric classes. For instance, almost 70% of the flat pixels in MERIT-DEM also appear to be flat in 3DEP-1, while 10% are overlapping either with foothill or shoulder which are indeed likely spatial neighbours to flat. This is possibly an indication for some co-registration or interpolation issues affecting the two products.

The other matrix is calculated to display the percentages with respect to the 3DEP-1 classification, i.e. each column will sum up to 100 (see Fig. 5d). It shows the likelihood with which a 3DEP-1 class appears in the same or other classes of MERIT, which is called “producer accuracy”, e.g. 86% of flat pixels in 3DEP-1 are also correctly classified as flat in MERIT. Additionally, the most likely confusion here is again with foothill ou alors shoulder classes (both around 6%), which support the above assumption of a co-registration issue. Another interesting finding is the widespread confusion between the summit et ridge classer.

The congruence of summit pixels between the two products is in fact less likely than their respective confusion with ridges in the other. In addition, there are at least three times as many summit pixels detected in 3DEP-1 compared to MERIT. Similar anomalies occur for the morphologically inverse classes depression et valley. The reason for these results could either be due to an increased richness of detail (or actual resolution) offered by 3DEP-1 or a higher level of noise (though the latter being less likely given its high level of detail). Nevertheless, this preliminary analysis shows that the underlying DEM data yield significantly different geomorphometric characteristics and that the confusion matrix allows these differences to be numerically expressed.

Comparing MERIT-DEM vs LiDAR elevation

Last return points in LiDAR data, which penetrate dense vegetation, are used to extract the DTM, whereas the first returns that hit the canopy of vegetation are used to derive the DSM. In our case, the LiDAR DTM and DSM were used to assess the quality of the tree height removal procedure carried out for the MERIT-DEM. Figure 6 reports the DTM vs. MERIT-DEM (red points) and the DSM vs. MERIT-DEM (blue points) of four study areas in USA with a high forest cover. In the four scatter plots, it is possible to distinguish the height difference between the DTM and DSM. The MERIT-DEM dataset has been corrected for the tree height bias 17 , and consequently the MERIT-DEM elevation values are expected to be closer to the LiDAR DTM than those of the LiDAR DSM. The LiDAR DTM vs. MERIT, and the LiDAR DSM vs. MERIT-DEM differences are analytically quantified by the linear model depicted in the scatter plot of Fig. 6. Where there is no vegetation or low vegetation (e.g. agricultural plains - Fig. 6c,d bare ground mountain tops - Fig. 6a,b), the DTM and DSM have almost identical values, denoted by an overlap of red and blue points. The presence of similar elevation values in the DTM and DSM causes a convergence of the linear models (blue and red lines the linear model functions are reported at the bottom of each scatter plot). The convergence in the lower parts of the plot (lower elevations) corresponds with the landscape (flat areas) and the absence of forest cover, which contributes to very similar values in the DTM and DSM. Whereas, at higher elevations, where there is increased forest cover, there is a greater difference between the DTM and DSM. This phenomena is more evident in Fig. 6c,d. On the contrary, Fig. 6a,b show a more parallel trend of blue and red lines, which is due to forest cover that is more equally distributed along the relief. It is important to note that the the scatter plots are plotted with different elevation ranges (x- and y-axis), and therefore the deviation of the red and blue line appears more evident in plots with lower elevation range (Fig. 6c). The regression coefficients for 3DEP-1 DTM vs. MERIT-DEM are slightly below 1 with an intercept value ranging from 8 to 34 m. On the other hand, the 3DEP-1 DSM vs. MERIT-DEM correlation has a regression coefficient larger than 1 and an intercept ranging from -154 to 46. The mean difference between DTM and DSM for Fig. 6a–d study areas ranges from 12 to 21 m. Overall, these values demonstrate the strong correlation between the MERIT-DEM and the LiDAR DTM.

MERIT-DEM vs LiDAR-DEM. Comparison of MERIT-DEM with the LiDAR DSM and LiDAR DTM for four study areas, represented by their scatter-plots and their relative linear models.

Identification of artefacts

It is important to note that the use and application of MERIT is based on the understanding that while it currently represents the best quality DEM available, it still contains errors and artefacts, which have not been corrected within the context of this research, as this was considered beyond the scope of the specific research objective. Consequently these errors cascade into the new Geomorpho90m dataset. The effect of these errors are mainly due to stripes that were recurrent in the AW3D product. These are visible in flat areas, where the artefact error is larger than the delta between pixels (e.g. slope). For instance, stripe artefacts can be found at the following locations in western Russia 63 as well as in central Russia 64 .


3 Answers 3

The answer is no, there is no such regular relationship between $R^2$ and the overall regression p-value, because $R^2$ depends as much on the variance of the independent variables as it does on the variance of the residuals (to which it is inversely proportional), and you are free to change the variance of the independent variables by arbitrary amounts.

As an example, consider any set of multivariate data $((x_, x_, ldots, x_, y_i))$ with $i$ indexing the cases and suppose that the set of values of the first independent variable, $<>>$, has a unique maximum $x^*$ separated from the second-highest value by a positive amount $epsilon$. Apply a non-linear transformation of the first variable that sends all values less than $x^* - epsilon/2$ to the range $[0,1]$ and sends $x^*$ itself to some large value $M gg 1$. For any such $M$ this can be done by a suitable (scaled) Box-Cox transformation $x o a((x-x_0)^lambda - 1)/(lambda-1))$, for instance, so we're not talking about anything strange or "pathological." Then, as $M$ grows arbitrarily large, $R^2$ approaches $1$ as closely as you please, regardless of how bad the fit is, because the variance of the residuals will be bounded while the variance of the first independent variable is asymptotically proportional to $M^2$.

You should instead be using goodness of fit tests (among other techniques) to select an appropriate model in your exploration: you ought to be concerned about the linearity of the fit and of the homoscedasticity of the residuals. And don't take any p-values from the resulting regression on trust: they will end up being almost meaningless after you have gone through this exercise, because their interpretation assumes the choice of expressing the independent variables did not depend on the values of the dependent variable at all, which is very much not the case here.


3 Answers 3

$R^2$ being zero for flat lines is not a bias: is just what it is intended to be. I'll give two reasons.

First reason is that $R^2$ measures if there is a linear relationship between two variables. A relationship would mean that the values we can expect from one variable depend on the value of the other variable. If the variables are related by a flat line, the value of one variable is always the same, it doesn't depend on the value of the other variable and there is no linear relationship between both variables. $R^2$ being zero or very small just shows it.

The second reason is that a flat line is a line just because of scale. If points are not exactly aligned, changing the scale of the vertical axe shows the difference between the points and the line and the scatterplot becomes a cloud and not a line.

Here is a graphical example:

In the bottom and right plots, points are clearly aligned, while upper left plot shows no relationship between the variables. However, all of them are just the same plot with different graphical scales.


One way to deal with this is with alpha blending, which makes each point slightly transparent. So regions appear darker that have more point plotted on them.

This is easy to do in ggplot2 :

Another convenient way to deal with this is (and probably more appropriate for the number of points you have) is hexagonal binning:

And there is also regular old rectangular binning (image omitted), which is more like your traditional heatmap:


6 Answers 6

The estimated value of the slope does not, by itself, tell you the strength of the relationship. The strength of the relationship depends on the size of the error variance, and the range of the predictor. Also, a significant $p$-value doesn't tell you necessarily that there is a strong relationship the $p$-value is simply testing whether the slope is exactly 0. For a sufficiently large sample size, even small departures from that hypothesis (e.g. ones not of practical importance) will yield a significant $p$-value.

Of the three quantities you presented, $R^2$, the coefficient of determination, gives the greatest indication of the strength of the relationship. In your case, $R^ <2>= .089$, means that $8.9\%$ of the variation in your response variable can be explained a linear relationship with the predictor. What constitutes a "large" $R^2$ is discipline dependent. For example, in social sciences $R^2 = .2$ might be "large" but in controlled environments like a factory setting, $R^2 > .9$ may be required to say there is a "strong" relationship. In most situations $.089$ is a very small $R^2$, so your conclusion that there is a weak linear relationship is probably reasonable.


Voir la vidéo: Korrelaatiokerroin ja ennuste