Berechnen des topografischen Robustheitsindex in ArcGIS Desktop?

Kann jemand den topografischen Robustheitsindex in ArcGIS Desktop ohne Zugriff auf die Befehlszeile von ArcInfo Workstation berechnen ?

Der topographische Rauheitsindex (TRI) ist eine Messung, die von Riley et al. (1999) entwickelt wurde, um den Höhenunterschied zwischen benachbarten Zellen eines digitalen Höhenrasters auszudrücken. Der Prozess berechnet im Wesentlichen den Höhenunterschied von einer mittleren Zelle und die acht Zellen, die es unmittelbar umgeben. Dann quadriert es jeden der acht Höhenunterschiedswerte, um sie alle positiv zu machen, und mittelt die Quadrate. Der topographische Robustheitsindex wird dann abgeleitet, indem die Quadratwurzel dieses Durchschnitts genommen wird und entspricht der durchschnittlichen Höhenänderung zwischen einem beliebigen Punkt in einem Raster und seiner Umgebung. " - Aus einem Aml-Drehbuch von Jeffrey Evans

Ich würde empfehlen, außerhalb von ArcGIS zu suchen.) Sehr einfach mit der kostenlosen GDAL-Software: http://www.gdal.org/gdaldem.html

gdaldem TRI input_dem output_TRI_map

Oder wenn Sie es in saga gis vorziehen würden: http://www.saga-gis.org/saga_modules_doc/ta_morphometry/ta_morphometry_16.html

Lassen Sie uns eine kleine (nur kleine) Algebra machen.

Sei x der Wert im zentralen Quadrat; Lassen Sie x_i, i = 1, .., 8 die Werte in den benachbarten Quadraten indizieren; und sei r der topographische Robustheitsindex. Dieses Rezept besagt, dass r ^ 2 gleich der Summe von (x_i - x) ^ 2 ist. Zwei Dinge, die wir leicht berechnen können, sind: (i) die Summe der Werte in der Nachbarschaft, gleich s = Summe {x_i} + x; und (ii) die Summe der Quadrate der Werte gleich t = Summe {x_i ^ 2} + x ^ 2. (Dies sind Schwerpunktstatistiken für das ursprüngliche Raster und für sein Quadrat.)

Das Erweitern der Quadrate ergibt

r ^ 2 = Summe {(x_i - x) ^ 2}

= Summe {x_i ^ 2 + x ^ 2 - 2 * x * x_i}

= Summe {x_i ^ 2} + 8 * x ^ 2 - 2 * x * Summe {x_i}

= [Summe {x_i ^ 2} + x ^ 2] + 7 * x ^ 2 - 2 * x * [Summe {x_i} + x - x]

= t + 7 * x ^ 2 - 2 * x * [Summe {x_i} + x] + 2 * x ^ 2

= t + 9 * x ^ 2 - 2 * x * s .

Stellen Sie sich zum Beispiel eine Nachbarschaft vor

1 2 3
4 5 6
7 8 9

Hier ist x = 5, s = 1 + 2 + ... + 9 = 45 und t = 1 + 4 + 9 + ... + 81 = 285. Dann

(1-5) ^ 2 + (2-5) ^ 2 + ... + (9-5) ^ 2 = 16 + 9 + 4 + 1 + 1 + 4 + 9 + 16 = 60 = r ^ 2

und die algebraische Äquivalenz sagt

60 = r ^ 2 = 285 + 9 * 5 ^ 2 -2 * 5 * 45 = 285 + 225 - 450 = 60, was überprüft.

Der Workflow ist daher:

Gegeben ein DEM.

• Berechnen Sie s = Fokalsumme (über 3 x 3 quadratische Nachbarschaften) von [DEM].

• Berechne DEM2 = [DEM] * [DEM].

• Berechnen Sie t = Fokalsumme (über 3 x 3 quadratische Nachbarschaften) von [DEM2].

• Berechnen Sie r2 = [t] + 9 * [DEM2] - 2 * [DEM] * [s].

Rückgabe r = Sqrt ([r2]).

Dies besteht aus insgesamt 9 Netzoperationen , die alle schnell sind. Sie können problemlos im Raster-Rechner (ArcGIS 9.3 und früher), in der Befehlszeile (alle Versionen) und in Model Builder (alle Versionen) ausgeführt werden.

Übrigens handelt es sich nicht um eine "durchschnittliche Höhenänderung" (da Höhenänderungen positiv und negativ sein können): Es handelt sich um eine quadratische mittlere Höhenänderung. Es entspricht nicht dem "topografischen Positionsindex", der unter http://arcscripts.esri.com/details.asp?dbid=14156 beschrieben ist und (laut Dokumentation) x - (s - x) / 8 entspricht. Im obigen Beispiel ist der TPI gleich 5 - (45-5) / 8 = 0, während der TRI, wie wir gesehen haben, Sqrt (60) ist.

Der TRI von Riley et al. (1999) ist die Quadratwurzel der summierten quadratischen Abweichungen. Dies kommt einer nicht skalierten Varianz sehr nahe. Wenn Sie eine Implementierung von Rileys TRI wünschen, folgen Sie bitte der von @whuber beschriebenen Methodik (die von @ user3338736 bereitgestellte Methodik verallgemeinerte die Metrik auf das Maximum im Fenster und repräsentiert nicht die zellenweise Variation).

In unserer ArcGIS-Toolbox "Geomorphometrie und Verlaufsmetriken" gibt es eine TRI-Variante , die die Varianz eines angegebenen Fensters angibt. Ich finde das flexibler und vertretbarer. Es gibt auch einige andere Metriken für die Oberflächenkonfiguration, einschließlich Robustheit und Dissektion.

-Bearbeiten: Die folgenden Informationen sind falsch. Bitte lesen Sie den Beitrag von whuber und erläutern Sie den korrekten Vorgang .....

TRI (Riley 1999) und TPI (Jenness 2002) sind ähnlich, aber unterschiedlich.

So berechnen Sie TRI und TPI mit ArcGIS 10.x ...

Schritt 1: Verwenden Sie das Tool "Focal Statistics", um zwei neue Raster-Datasets aus einem DEM zu erstellen.

Raster 1 "MAX") Nachbarschaft: Rechteck, Höhe: 3, Breite: 3, Einheiten: Zelle, Statistiktyp: Maximum

Raster 2 "MIN") Nachbarschaft: Rechteck, Höhe: 3, Breite: 3, Einheiten: Zelle, Statistiktyp: Minimum

Schritt 2: Verwenden Sie den Raster-Rechner, um die folgenden Funktionen für die gerade erstellten 2 Raster-Datasets auszuführen.

Für TRI: SquareRoot (Abs ((Square ("% MAX%") - Square ("% MIN%")))

Für TPI: ("% Input DEM%" - "% MIN%") / ("% MAX%" - "% MIN%")

Hier ist ein Python-Beispielcode, der aus einem Modell exportiert wurde, das ich für TRI erstellt habe.

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# script.py
# Created on: 2014-03-06 08:56:13.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: script <Input_raster> <TRI_Raster> 
# Description: 
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Script arguments
Input_raster = arcpy.GetParameterAsText(0)

TRI_Raster = arcpy.GetParameterAsText(1)
if TRI_Raster == '#' or not TRI_Raster:
    TRI_Raster = "C:\\Users\\Documents\\ArcGIS\\Default.gdb\\rastercalc1" # provide a default value if unspecified

# Local variables:
MIN = Input_raster
MAX = Input_raster

# Process: 3x3Max
arcpy.gp.FocalStatistics_sa(Input_raster, MAX, "Rectangle 3 3 CELL", "MAXIMUM", "DATA")

# Process: 3x3Min
arcpy.gp.FocalStatistics_sa(Input_raster, MIN, "Rectangle 3 3 CELL", "MINIMUM", "DATA")

# Process: Raster Calculator
arcpy.gp.RasterCalculator_sa("SquareRoot(Abs((Square(\"%MAX%\") - Square(\"%MIN%\"))))", TRI_Raster)

Das klingt sehr nach dem Topografischen Positionsindex, den ich kürzlich für eines meiner Projekte verwendet habe. Es gibt ein ArcScript auf der ESRI-Support-Seite, eine Topografie-Toolbox auf der ESRI Resource Center-Seite und weitere Informationen zum Vorgang auf der Jenness Enterprises- Seite.