Obliczanie indeksu wytrzymałości topograficznej w ArcGIS Desktop?

Czy ktoś wie, jak obliczyć Indeks odporności topograficznej w ArcGIS Desktop bez dostępu do linii poleceń ArcInfo Workstation?

„Wskaźnik wytrzymałości topograficznej (TRI) to pomiar opracowany przez Riley i wsp. (1999) w celu wyrażenia wielkości różnicy wysokości między sąsiednimi komórkami cyfrowej siatki wysokości. Proces ten zasadniczo oblicza różnicę wartości wysokości z komórki środkowej i ośmioma komórkami bezpośrednio go otaczającymi. Następnie obciąża każdą z ośmiu wartości różnicy wysokości, aby wszystkie były dodatnie, i uśrednia kwadraty. Wskaźnik chropowatości topograficznej jest następnie uzyskiwany na podstawie pierwiastka kwadratowego z tej średniej i odpowiada średniej zmianie wysokości między dowolnym punktem na siatce a otaczającym obszarem. ” - z arcriptu amla autorstwa Jeffreya Evansa

Polecam szukać poza ArcGIS) Bardzo łatwe w użyciu darmowego oprogramowania gdal: http://www.gdal.org/gdaldem.html

gdaldem TRI input_dem output_TRI_map

Lub jeśli wolisz w saga gis: http://www.saga-gis.org/saga_modules_doc/ta_morphometry/ta_morphometry_16.html

Zróbmy małą (tylko trochę) algebrę.

Niech x będzie wartością w centralnym kwadracie; niech x_i, i = 1, .., 8 indeksuje wartości w sąsiednich kwadratach; i niech r będzie topograficznym wskaźnikiem odporności. Ten przepis mówi, że r ^ 2 równa się sumie (x_i - x) ^ 2. Dwie rzeczy, które możemy łatwo obliczyć, to (i) suma wartości w sąsiedztwie, równa s = Suma {x_i} + x; oraz (ii) suma kwadratów wartości, równa t = Suma {x_i ^ 2} + x ^ 2. (Są to ogniskowe statystyki oryginalnej siatki i jej kwadratu).

Rozszerzanie kwadratów daje

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

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

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

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

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

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

Na przykład rozważ okolicę

1 2 3
4 5 6
7 8 9

Tutaj x = 5, s = 1 + 2 + ... + 9 = 45, a t = 1 + 4 + 9 + ... + 81 = 285. Następnie

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

a algebraiczna równoważność mówi

60 = r ^ 2 = 285 + 9 * 5 ^ 2 -2 * 5 * 45 = 285 + 225 - 450 = 60, co sprawdza.

Przepływ pracy jest zatem następujący:

Biorąc pod uwagę DEM.

• Oblicz s = suma ogniskowa (ponad 3 x 3 dzielnice kwadratowe) [DEM].

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

• Oblicz t = suma ogniskowa (ponad 3 x 3 dzielnice kwadratowe) z [DEM2].

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

Zwraca r = Sqrt ([r2]).

Składa się z 9 operacji siatki w całości , z których wszystkie są szybkie. Są one łatwo przeprowadzane w kalkulatorze rastrowym (ArcGIS 9.3 i wcześniejszych), w wierszu poleceń (wszystkie wersje) i w Konstruktorze modeli (wszystkie wersje).

BTW, nie jest to „średnia zmiana wysokości” (ponieważ zmiany wysokości mogą być dodatnie i ujemne): jest to podstawowa zmiana średniej wysokości kwadratu. To jest , nie są równe z „topograficznych wskaźnika położenia” opisanej w http://arcscripts.esri.com/details.asp?dbid=14156 , który (zgodnie z dokumentacją) jest równa x - y - (x) / 8. W powyższym przykładzie TPI wynosi 5 - (45-5) / 8 = 0, podczas gdy TRI, jak widzieliśmy, to Sqrt (60).

Riley i wsp., (1999) TRI jest pierwiastkiem kwadratowym z sumowanych kwadratowych odchyleń. Jest to bardzo zbliżone do nieskalowanej wariancji. Jeśli chcesz zaimplementować TRI Riley, postępuj zgodnie z metodologią przedstawioną przez @whuber (metodologia podana przez @ user3338736 uogólniła metrykę na maksimum w oknie i nie reprezentuje komórki według wariantu komórki).

Mam odmianę TRI w naszym Geomorfometrii i Gradient Metrics ArcGIS Toolbox, która jest wariantem określonego okna. Uważam to za bardziej elastyczne i uzasadnione. Istnieją również inne parametry konfiguracji powierzchni, w tym chropowatość i rozwarstwienie.

-Edytuj: poniższe informacje są nieprawidłowe. Proszę zobaczyć post whuber wyjaśniający prawidłowy proces .....

TRI (Riley 1999) i TPI (Jenness 2002) są podobne, ale różne.

Aby obliczyć TRI i TPI za pomocą ArcGIS 10.x ...

Krok 1: Użyj narzędzia Focal Statistics, aby utworzyć 2 nowe zestawy danych rastrowych z DEM.

Raster 1 „MAX”) Okolica: prostokąt, wysokość: 3, szerokość: 3, jednostki: komórka, typ statystyki: maksymalna

Raster 2 „MIN”) Okolica: prostokąt, wysokość: 3, szerokość: 3, jednostki: komórka, typ statystyki: minimum

Krok 2: Użyj kalkulatora rastrowego, aby wykonać następujące funkcje na 2 właśnie utworzonych zestawach danych rastrowych.

Dla TRI: SquareRoot (Abs ((Kwadrat („% MAX%”) - Kwadrat („% MIN%”))))

Dla TPI: („% Input DEM%” - „% MIN%”) / („% MAX%” - „% MIN%”)

Oto przykładowy kod Pythona wyeksportowany z modelu zbudowanego dla TRI ....

# -*- 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)

Brzmi to bardzo podobnie do wskaźnika pozycji topograficznej, procesu, którego ostatnio użyłem w jednym z moich projektów. Na stronie wsparcia ESRI znajduje się ArcScript , przybornik Topografia na stronie Centrum zasobów ESRI oraz kilka dodatkowych informacji na temat procesu na stronie Jenness Enterprises .