Ich habe SRTM GDEM (~ 90 km Auflösung) heruntergeladen.

Ich verwende ArcGIS 10.

Ich habe versucht, mithilfe eines räumlichen Analysten die Steigung zu berechnen.

Allerdings kann ich die Steigung nicht berechnen.

Der Ausgabewert hat nur zwei Bereiche 0 und 0.1-90.

Ich bin mir nicht sicher, wo das Problem liegt.

Dies scheint ein guter Ort zu sein, um eine einfache, schnelle und mehr als einigermaßen genaue Methode zur Berechnung von Steigungen für ein global ausgedehntes DEM zu beschreiben .

Prinzipien

Es sei daran erinnert, dass die Neigung einer Oberfläche an einem Punkt im Wesentlichen das größte Verhältnis von "Anstieg" zu "Lauf" ist, das bei allen möglichen Peilungen von diesem Punkt aus auftritt. Das Problem ist, dass bei einer Projektion mit Skalenverzerrung die Werte von "run" falsch berechnet werden. Schlimmer noch, wenn die Maßstabsverzerrung mit der Peilung variiert - was bei allen nicht konformen Projektionen der Fall ist -, wird die Abweichung der Steigung mit der Peilung falsch geschätzt, wodurch eine genaue Ermittlung des maximalen Anstiegs - Lauf - Verhältnisses (und eine Verdrehung des Maßstabs) verhindert wird Berechnung des Aspekts).

Wir können dies lösen, indem wir eine konforme Projektion verwenden, um sicherzustellen, dass die Skalierungsverzerrung nicht mit der Peilung variiert, und dann die Neigungsschätzungen korrigieren , um die Skalierungsverzerrung zu berücksichtigen (die von Punkt zu Punkt auf der Karte variiert). Der Trick besteht darin, eine globale konforme Projektion zu verwenden, die einen einfachen Ausdruck für die Skalenverzerrung ermöglicht.

Die Mercator-Projektion passt zur Rechnung: Wenn der Maßstab am Äquator korrekt ist, entspricht seine Verzerrung der Sekante des Breitengrads. Das heißt, Entfernungen auf der Karte scheinen mit dem Sekanten multipliziert zu werden. Dies veranlasst jede Steigungsberechnung, einen Anstieg zu berechnen: (sec (f) * run) (was ein Verhältnis ist), wobei f der Breitengrad ist. Um dies zu korrigieren, müssen wir die berechneten Steigungen mit sec (f) multiplizieren. oder gleichwertig durch cos (f) dividieren. Dies gibt uns das einfache Rezept:

Berechnen Sie die Steigung (als Anstieg: Laufen oder Prozent) mit einer Mercator-Projektion und dividieren Sie das Ergebnis durch den Kosinus des Breitengrads.

Arbeitsablauf

Führen Sie dazu mit einem in Dezimalgraden angegebenen Raster (z. B. einem SRTM-DEM) die folgenden Schritte aus:

1. Erstellen Sie ein Breitengradraster. (Dies ist nur das y-Koordinatengitter.)

2. Berechnen Sie den Cosinus.

3. Projekt sowohl die DIE und die Kosinus der Breite einer Mercator - Projektion unter Verwendung von in dem Maßstab am Äquator wahr ist.

4. Falls erforderlich, rechnen Sie die Höheneinheiten so um, dass sie mit den Einheiten der projizierten Koordinaten (normalerweise Meter) übereinstimmen.

5. Berechnen Sie die Steigung des projizierten DEM entweder als reine Steigung oder als Prozent ( nicht als Winkel).

6. Teilen Sie diese Steigung durch das projizierte Kosinusgitter (Breitengrad).

7. Projizieren Sie das Steigungsgitter zur weiteren Analyse oder Zuordnung gegebenenfalls auf ein anderes Koordinatensystem.

Die Fehler in den Steigungsberechnungen betragen bis zu 0,3% (da bei diesem Verfahren ein sphärisches Erdmodell anstelle eines ellipsoiden Modells verwendet wird, das um 0,3% abgeflacht ist). Dieser Fehler ist wesentlich kleiner als andere Fehler, die in die Steigungsberechnung einfließen und daher vernachlässigt werden können.

Vollständig globale Berechnungen

Die Mercator-Projektion kann nicht mit beiden Polen umgehen. Verwenden Sie für Arbeiten in Polarregionen eine polare stereografische Projektion mit echter Skalierung am Pol. Die Skalenverzerrung beträgt 2 / (1 + sin (f)). Verwenden Sie diesen Ausdruck anstelle von sec (f) im Workflow. Anstatt ein Kosinusgitter (Breitengrad) zu berechnen, berechnen Sie insbesondere ein Gitter mit den Werten (1 + sin (Breitengrad)) / 2 ( Bearbeiten : Verwenden Sie -latitude für den Südpol, wie in den Kommentaren erläutert). Gehen Sie dann genauso vor wie zuvor.

Um eine vollständige globale Lösung zu erhalten, sollten Sie das Erdgitter in drei Teile aufteilen - einen um jeden Pol und einen um den Äquator -, in jedem Teil eine separate Steigungsberechnung mit einer geeigneten Projektion durchführen und die Ergebnisse mosaikieren. Ein vernünftiger Ort, um den Globus zu spalten, liegt in Breitengraden von 2 * ArcTan (1/3), was ungefähr 37 Grad entspricht, da in diesen Breitengraden der Mercator- und der Stereographic-Korrekturfaktor gleich sind (einen gemeinsamen Wert haben) von 5/4) und es wäre schön, die Größe der vorgenommenen Korrekturen zu minimieren. Zur Überprüfung der Berechnungen sollten sich die Gitter in sehr enger Übereinstimmung befinden, wo sie überlappen (winzige Mengen an Gleitkommaungenauigkeit und Differenzen aufgrund der Neuabtastung der projizierten Gitter sollten die einzigen Ursachen für Abweichungen sein).

Verweise

John P. Snyder, Kartenprojektionen - Ein Arbeitshandbuch . USGS Professional Paper 1395, 1987.

Ursprüngliche Antwort

Ich vermute, die horizontalen Einheiten für Ihr Raster sind entweder in Grad oder Bogensekunden. Sie müssen dieses Raster in eine räumliche Projektion umprojizieren, bei der die horizontalen und vertikalen Einheiten identisch sind (dh wenn die vertikalen Einheiten in Metern angegeben sind, empfehle ich die Verwendung von UTM mit horizontalen Einheiten in Metern).

So projizieren Sie ein Raster mit ArcCatalog / ArcGIS neu:

ArcToolbox> Datenverwaltungstools> Projektionen und Transformationen> Raster> Projektraster

Wählen Sie einen projizierten Raumbezug, der Ihre Region von Interesse abdeckt, z. B. eine UTM-Zone. Es gibt viele andere Optionen, die am besten im Handbuch dokumentiert sind . Beachten Sie, dass Sie kein Steigungs-Dataset für die gesamte Erde erstellen können (wenn Sie dies versuchen).

Bessere Antwort unter Verwendung von GDAL mit einer Skala

Jetzt, da SRTM-Daten global verfügbar sind, kann ich die Dateien tatsächlich sehen und damit arbeiten. Das gdaldemDienstprogramm von GDAL kann Neigung und Schatten mithilfe einer Skalierungsoption für ein Verhältnis von vertikalen zu horizontalen Einheiten berechnen . Das Handbuch empfiehlt 111120 m / ° für etwas wie SRTM-Fliesen. Zum Beispiel aus einer OSGeo4W-Shell:

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

Mit dieser -compute_edgesOption werden die Kanten nahtloser, wenn Sie einige Kacheln zusammenfügen möchten. Oder berechnen Sie Kacheln für eine große Region. Der Nachteil der "Skalentechnik" besteht darin, dass die Abstände in EW- und NS-Richtung mit Ausnahme des Äquators nicht gleich sind. Daher kann es bei Kacheln, die näher an den Polen liegen, zu merkwürdigen Falschdarstellungen der Neigung kommen.

Einfach ausgedrückt, es gibt keinen. Per Definition wird ein auf Grad basierendes Koordinatensystem nicht projiziert. Im allgemeinen Sprachgebrauch sagen wir, dass WGS84 eine "geografische" Projektion ist, aber das ist nicht wahr, nur aus Bequemlichkeitsgründen.

Ich erinnere mich, dass ich etwas über eine Software oder einen Prozess zum genauen Arbeiten mit Höhenmodellen im nicht projizierten geografischen Raum gelesen habe, aber ich kann ihn derzeit nicht finden. In jedem Fall wäre es ein Experiment gewesen oder würde es selbst aus einem Code-Prozess erstellt.

Ahhh, gefunden: Entwicklung eines globalen Slope-Datensatzes zur Abschätzung des Auftretens von Erdrutschen infolge von Erdbeben (USGS). Seite 4 beschreibt das Problem gut

... die Länge eines Grades variiert je nach Breitengrad. Am Äquator ist ein Block von einem Grad zu einem Grad einigermaßen quadratisch, wenn er in Einheiten von Metern umgerechnet wird (111.321 Meter in x-Richtung und 110.567 Meter in y-Richtung ... aber näher an den Polen liegen die Abstände in Die x-Richtung wird aufgrund der Konvergenz der Meridiane als Funktion des Breitengrads kleiner. Die meisten GIS-Pakete, einschließlich ArcGIS, werden nur mit quadratischen Pixeln ausgeführt und verwenden daher einen Faktor zum Anpassen der x-, y- oder z-Dimensionen an eine gemeinsame einheit ist nicht möglich.

In diesem Artikel werden die spezifischen Berechnungen und Softwaretools ( gdal , python , numpy ) beschrieben, mit denen dieses grundlegende Problem umgangen wurde . Das Papier enthält nicht den Code, aber wenn sie freundlich gefragt werden, können sie teilen. Auf jeden Fall würde ich wahrscheinlich nur fragen, wo die Ergebnisse sind. Als USGS ist es wahrscheinlich schon irgendwo online. :)

Globale DEM-Parameter (wobei die meisten Formeln auf der Annahme eines euklidischen Raums basieren) können mithilfe des EQUI7-GRID-Systems effizient abgeleitet werden (Bauer-Marschallinger et al. 2014). EQUI7 GRID teilt die Welt in 7 Landgebiete auf, die alle mit einem äquidistanten Projektionssystem mit minimalem Präzisionsverlust projiziert werden. Siehe ein Beispiel für globales DEM bei einer Auflösung von 250 m im EQUI7-GITTER. Hier finden Sie einen Beispielcode , der zeigt, wie Sie globale DEM-Parameter mit SAGA GIS ableiten. Sobald Sie die Ableitung der DEM-Parameter im EQUI7-GRID-System abgeschlossen haben, können Sie alle Karten in WGS84- longlatKoordinaten zurücktransformieren und dann mit GDAL ein globales Mosaik erstellen.

Piste ist Steigen / Laufen. Berechnen Sie rise und run und Sie haben Ihre Antwort. Es ist einfach, die Entfernung zwischen geografischen Koordinaten zu berechnen. Dies führt zu weniger Resampling-Fehlern im Vergleich zur Umstellung auf UTM usw.