Vous utilisez SRTM Global DEM pour le calcul de pente?

17

J'ai téléchargé SRTM GDEM (~ 90 km de résolution).

J'utilise ArcGIS 10.

J'ai essayé d'utiliser l'analyste spatial pour calculer la pente.

Cependant, je ne peux pas calculer la pente.

Les valeurs de sortie n'ont que deux plages 0 et 0,1-90.

Je ne sais pas vraiment quel est le problème?

user2543
la source
Cela dépend de l'endroit où vous analysez dans le monde. Il existe différentes projections pour chaque emplacement. Où examinez-vous?
djq
5
La résolution est en fait de ~ 90 m, pas de ~ 90 km.
Akheloes
Juste un commentaire, si vous êtes en maintenance pour Desktop, vous pouvez vous connecter à ArcGIS Online et utiliser leurs services d'élévation (sans besoin d'extension NA). La couche de pente est libre d'utiliser comme couche de référence. En Australie, nous avons les données 1 seconde SRTM (~ 30m res) blogs.esri.com/esri/arcgis/2014/07/11/…
Simon

Réponses:

29

Cela semble être un bon endroit pour décrire un moyen simple, rapide et plus que raisonnablement précis de calculer les pentes pour un DEM global .

Des principes

Rappelons que la pente d'une surface en un point est essentiellement le plus grand rapport "montée" / "course" rencontré à tous les relèvements possibles à partir de ce point. Le problème est que lorsqu'une projection présente une distorsion d'échelle, les valeurs de "run" seront incorrectement calculées. Pire encore, lorsque la distorsion d'échelle varie avec le relèvement - ce qui est le cas pour toutes les projections qui ne sont pas conformes - la façon dont la pente varie avec le relèvement sera incorrectement estimée, ce qui empêchera une identification précise du rapport élévation / course maximal (et faussera le calcul de l'aspect).

Nous pouvons résoudre ce problème en utilisant une projection conforme pour garantir que la distorsion d'échelle ne varie pas avec le relèvement, puis en corrigeant les estimations de pente pour tenir compte de la distorsion d'échelle (qui varie d'un point à un autre sur la carte). L'astuce consiste à utiliser une projection conforme globale qui permet une expression simple pour sa distorsion d'échelle.

La projection de Mercator fait l'affaire: en supposant que l'échelle est correcte à l'équateur, sa distorsion est égale à la sécante de la latitude. Autrement dit, les distances sur la carte semblent être multipliées par la sécante. Cela fait que tout calcul de pente calcule l'élévation: (sec (f) * run) (qui est un rapport), où f est la latitude. Pour corriger cela, nous devons multiplier les pentes calculées par sec (f); ou, de manière équivalente, divisez-les par cos (f). Cela nous donne la recette simple:

Calculez la pente (comme montée: course ou pourcentage) à l'aide d'une projection de Mercator, puis divisez le résultat par le cosinus de la latitude.

Workflow

Pour ce faire, avec une grille indiquée en degrés décimaux (comme un DEM SRTM), procédez comme suit:

  1. Créez une grille de latitude. (Ceci est juste la grille de coordonnées y.)

  2. Calculez son cosinus.

  3. Projetez à la fois le DEM et le cosinus de la latitude à l'aide d'une projection de Mercator dans laquelle l'échelle est vraie à l'équateur.

  4. Si nécessaire, convertissez les unités d'élévation pour qu'elles correspondent aux unités des coordonnées projetées (généralement des mètres).

  5. Calculez la pente du DEM projeté sous forme de pente pure ou de pourcentage ( pas sous forme d'angle).

  6. Divisez cette pente par la grille de cosinus (latitude) projetée.

  7. Si vous le souhaitez, reprojetez la grille de pente à tout autre système de coordonnées pour une analyse ou une cartographie plus poussée.

Les erreurs dans les calculs de pente vont jusqu'à 0,3% (car cette procédure utilise un modèle de terre sphérique plutôt qu'un modèle ellipsoïdal, qui est aplati de 0,3%). Cette erreur est sensiblement plus petite que les autres erreurs qui entrent dans les calculs de pente et peut donc être négligée.


Calculs entièrement globaux

La projection Mercator ne peut gérer aucun des deux pôles. Pour travailler dans des régions polaires, pensez à utiliser une projection stéréographique polaire avec une échelle réelle au pôle. La distorsion d'échelle est égale à 2 / (1 + sin (f)). Utilisez cette expression à la place de sec (f) dans le flux de travail. Plus précisément, au lieu de calculer une grille de cosinus (latitude), calculez une grille dont les valeurs sont (1 + sin (latitude)) / 2 ( modifier : utilisez -latitude pour le pôle Sud, comme indiqué dans les commentaires). Procédez ensuite exactement comme avant.

Pour une solution globale complète, envisagez de diviser la grille terrestre en trois parties - une autour de chaque pôle et une autour de l'équateur -, en effectuant un calcul de pente séparément dans chaque partie en utilisant une projection appropriée et en mosaïquant les résultats. Un endroit raisonnable pour diviser le globe est le long de cercles de latitude à des latitudes de 2 * ArcTan (1/3), ce qui est d'environ 37 degrés, car à ces latitudes les facteurs de correction Mercator et stéréographique sont égaux les uns aux autres (ayant une valeur commune 5/4) et ce serait bien de minimiser la taille des corrections apportées. Pour vérifier les calculs, les grilles doivent être en accord très étroit où elles se chevauchent (de minuscules quantités d'imprécision en virgule flottante et les différences dues au rééchantillonnage des grilles projetées devraient être les seules sources de divergences).

Les références

John P. Snyder, Projections cartographiques - Un manuel de travail . USGS Professional Paper 1395, 1987.

whuber
la source
7
Je me trouve dans la position, comme je le fais si souvent, de remercier encore une fois whuber pour avoir décrit une solution ainsi que pour avoir avancé le raisonnement qui la construit. Mon chapeau vous est retiré monsieur.
Matt Wilkie
Merci @matt. Je ne voulais pas suggérer plus tôt que votre réponse (maintenant supprimée) devrait être annulée: en fait, je l'avais votée parce que vous partagiez un lien vers une référence USGS intéressante qui pourrait être utile à de nombreux lecteurs. (Mon commentaire ne critiquait qu'un passage secondaire dans ce journal, pas le journal lui-même.)
whuber
ahh. Merci pour la clarification. J'ai rétabli la réponse, en faisant confiance aux gens qui ont suffisamment d'informations devant eux maintenant pour faire un choix éclairé :)
matt wilkie
2
Venant d'un milieu français, il m'a fallu un certain temps pour traduire la terminologie nécessaire pour mieux comprendre cette excellente réponse, j'ai donc pensé que supprimer ce lien était une bonne aide pour les débutants comme moi: webhelp.esri.com/arcgisdesktop/9.2/…
Akheloes
C'est une excellente approche et j'ai déjà utilisé votre solution pour générer un raster de pente global. Un indice de l'expérience pratique: les valeurs de latitude au sud de l'équateur étant négatives, vous devez utiliser la valeur de latitude absolue dans l'équation suivante: (1 + sin (latitude)) / 2
Saleika
18

Réponse originale

Je suppose que les unités horizontales de votre raster sont en degrés ou en secondes d'arc. Vous devez reprojeter ce raster à une projection spatiale où vos unités horizontales et verticales sont les mêmes (c.-à-d., Si les unités verticales sont en mètres, alors je suggère d'utiliser UTM, qui a des unités horizontales de mètres).

Pour reprojeter un raster avec ArcCatalog / ArcGIS, regardez dans:

ArcToolbox> Outils de gestion des données> Projections et transformations> Raster> Project Raster

Choisissez une référence spatiale projetée qui couvre votre région d'intérêt, par exemple, essayez une zone UTM. Il existe de nombreuses autres options qui sont mieux documentées dans le manuel . Remarque, vous ne pouvez pas créer un jeu de données de pente pour la Terre entière (si c'est ce que vous essayez de faire).

Meilleure réponse, en utilisant GDAL avec une échelle

Maintenant que les données SRTM sont disponibles dans le monde entier, je peux réellement voir et travailler avec les fichiers. L' gdaldemutilitaire de GDAL peut calculer la pente et l'ombrage à l'aide d'une option d' échelle pour un rapport des unités verticales à l'horizontale. Le manuel recommande 111120 m / ° pour quelque chose comme les carreaux SRTM. Ainsi, par exemple, à partir d'un shell OSGeo4W:

$ gdaldem slope -s 111120 -compute_edges N44E007.hgt N44E007_slope.tif

L' -compute_edgesoption rend les bords plus transparents, si vous souhaitez assembler quelques carreaux ensemble. Ou calculez des tuiles pour une grande région. L'inconvénient de la technique de "l'échelle" est que les distances dans les directions EW et NS ne sont pas égales, sauf à l'équateur, donc pour les carreaux plus proches des pôles, il pourrait y avoir quelques fausses représentations étranges de la pente.

Mike T
la source
Il vaut la peine de souligner votre dernier commentaire: c'est une mauvaise solution pour les points non proches de l'équateur. Ce n'est pas une mince affaire de «fausses déclarations étranges»: les résultats seront grossièrement erronés, surtout dans des endroits plus proches des Polonais que de l'Équateur. La documentation pour les gdaldemétats "Pour les emplacements non proches de l'équateur, il serait préférable de reprojeter votre grille en utilisant gdalwarp avant d'utiliser gdaldem." Malheureusement, cela ne fonctionnera pas pour les ensembles de données couvrant le globe, sauf si vous les divisez en petits morceaux (74 zones UTM, peut-être?), Les projetez, calculez les pentes et mosaïquez les résultats.
whuber
7

Autrement dit, il n'y en a pas. Par définition, un système de coordonnées basé sur les degrés n'est pas projeté. Dans le langage courant, nous disons que WGS84 est une projection "géographique", mais c'est faux, juste pour plus de commodité.

Je pense que je me souviens d'avoir lu un logiciel ou un processus pour travailler avec précision avec des modèles d'altitude dans un espace géographique non projeté, mais je ne peux pas le localiser pour le moment. Dans tous les cas, il aurait été expérimental ou construit vous-même à partir d'un processus de type code.


Ahhh, l'a trouvé: Développement d'un ensemble de données sur les pentes mondiales pour l'estimation de l'occurrence des glissements de terrain résultant des tremblements de terre (USGS). La page 4 décrit bien le problème

... la longueur d'un degré varie en fonction de sa position latitudinale. À l'équateur, un bloc d'un degré par un degré est raisonnablement carré lorsqu'il est converti en unités de mètres (111 321 mètres dans la direction x par 110 567 mètres dans la direction y ... mais plus près des pôles les distances dans le La direction x devient plus petite en fonction du cosinus de latitude, en raison de la convergence des méridiens. La plupart des packages SIG, ArcGIS inclus, fonctionnent uniquement sur des pixels carrés, et donc en utilisant un facteur pour ajuster les dimensions x, y ou z à une unité commune n'est pas possible.

L'article continue en décrivant les calculs spécifiques et les outils logiciels ( , , ) qu'ils ont utilisés pour contourner ce problème fondamental. Le document n'inclut pas le code, mais si on le lui demande gentiment, il pourrait le partager. Quoi qu'il en soit, je demanderais probablement où sont les résultats, étant l'USGS, il est probablement déjà en ligne quelque part. :)

Matt Wilkie
la source
1
La suggestion de cet article selon laquelle une projection équidistante azimutale pourrait être utilisée pour calculer les pentes est erronée et fausse. Il donnera en effet des pentes correctes près de l'origine de la projection, mais elles aussi deviendront progressivement moins précises à mesure que la distance à l'origine augmente.
whuber
merci de l'avoir signalé. Lecteurs, assurez-vous de lire également gis.stackexchange.com/a/40464/108 , pour l'équilibre
Matt Wilkie
2

Les paramètres DEM mondiaux (où la plupart des formules sont basées sur l'hypothèse d'espace euclidien) peuvent être dérivés efficacement en utilisant le système EQUI7 GRID (Bauer-Marschallinger et al. 2014). EQUI7 GRID divise le monde en 7 zones terrestres, toutes projetées dans un système de projection équidistant avec une perte de précision minimale. Voir un exemple de DEM global à 250 m de résolution dans la grille EQUI7. Vous trouverez ici quelques exemples de code qui montrent comment dériver des paramètres DEM globaux à l'aide de SAGA GIS. Une fois que vous avez dérivé les paramètres DEM dans le système EQUI7 GRID, vous pouvez retransformer toutes les cartes en longlatcoordonnées WGS84 , puis créer un mosaïque global à l'aide de GDAL.

Tom Hengl
la source
Pourriez-vous expliquer comment cela répond à la question? Si vous proposez d'utiliser des projections équidistantes pour les calculs de pente, veuillez noter que c'est une mauvaise solution en raison des grandes distorsions métriques relatives qui se produisent lorsque l'on s'éloigne du centre de la projection. Bien que concentrer sept de ces projections sur les masses terrestres contribue à atténuer ce problème, ce n'est toujours pas le meilleur choix.
whuber
L'article de Bauer-Marschallinger et al. (2014) explique pourquoi ces projections ont été choisies pour représenter les masses terrestres mondiales (elles sont supposées avoir une perte de précision minimale). Je conviens que toute projection 2D conduira éventuellement à des déformations, mais à ma connaissance, EQUI7 est un bon compromis entre perte de précision et commodité (algèbre 2D). Cela dit, les hexagones sont à nouveau utilisés pour représenter les surfaces terrestres mondiales (bien que l'analyse DEM avec des hexagones 3D soit toujours lourde).
Tom Hengl
Merci pour la référence. Son résumé suggère qu'il résout un problème très différent, celui de «minimiser le suréchantillonnage des données locales apparaissant lors de la projection d'images satellites génériques sur une grille matricielle régulière». Cela ne signifie pas qu'il fonctionnera bien à d'autres fins, telles que l'estimation des pentes.
whuber
Bien sûr, EQUI7 ne résout pas absolument le problème de l'estimation précise des pentes locales, mais c'est probablement une solution plus élégante que d'utiliser la projection de Mercator suggérée ci-dessus. En fin de compte, si l'on souhaite estimer les pentes avec une précision parfaite, les seules options sont probablement (1) d'utiliser des projections locales (équidistantes) pour des carreaux de plus petite taille (par exemple 100 x 100 km) avec un chevauchement de 10 à 20%, comme mentionné également dans Verdin et al. (2007) ou (2) pour utiliser la grille hexagonale ( paquet dggridR ).
Tom Hengl
Le problème n'est pas la précision - il consiste à produire des pentes et des aspects systématiquement biaisés. Étant donné que les projections équidistantes déforment différentiellement les directions orthogonales aux géodésiques originaires de leurs centres, les aspects seront toujours erronés (bien que raisonnablement précis près des centres où toute distorsion est faible) et les erreurs dans les pentes augmenteront rapidement avec la distance. L'utilisation de nombreuses projections locales fonctionnera bien sûr, mais c'est tout le contraire de l'élégance que vous appréciez.
whuber
-2

La pente est montée / course. Calculez la montée et la course au calcul et vous avez votre réponse. Il est simple de calculer la distance entre les coordonnées géographiques. Cela introduira moins d'erreur de rééchantillonnage par rapport à la conversion en UTM, etc.

THK
la source