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

15

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

matowe wilkie
źródło
zależy od wersji ArcGIS arcscripts.esri.com/details.asp?dbid=12646 trochę dyskusji z poprzednich forów forums.esri.com/Thread.asp?c=93&f=982&t=145448 niezaznaczone, ale wyszukiwanie zawierało termin jennessent.com /arcgis/surface_area.htm

Odpowiedzi:

18

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

johanvdw
źródło
11
+1 Zawsze cenię sobie rozwiązania ArcGIS inne niż ArcGIS :-). Jest to kwestia zasady, a nie antagonizmu w stosunku do ArcGIS. Należy unikać wiązania się z pojedynczym oprogramowaniem: nie tylko jest ono ryzykowne zawodowo, ale także duszące intelektualnie.
whuber
Wiem, że poprosiłem o rozwiązanie specyficzne dla Arcgis, ale akceptuję to ze względu na jego bezpośredniość. Narzędzia GDAL są łatwe do zdobycia i instalacji, powszechnie uznawane za najlepsze w swojej klasie, a polecenie wygenerowania tego konkretnego produktu to definicja prostoty.
matt wilkie
18

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).

Whuber
źródło
1
Dziękuję Bill. Doceniam to, jak działa narzędzie lub operacja. Na tej podstawie każdy, kto ma odpowiednią inwestycję czasu i energii intelektualnej, może zbudować nowy aparat do wykonania tej pracy, korzystając z posiadanych narzędzi. Dzięki takim informacjom GIS.se będzie użyteczną usługą na dłuższą metę.
matt wilkie
1
+1 Świetne wyjaśnienie. Myślę, że to oznacza, że ​​stroma, ale gładka powierzchnia może mieć wyższy TRI niż płaska, ale wyboista powierzchnia.
Kirk Kuykendall
1
@Kirk To prawda. Istnieją sposoby na usunięcie efektu lokalnego nachylenia w celu uzyskania indeksu „względnej” odporności, jeśli chcesz. Chociaż nie opracowałem szczegółów, uważam, że odejmując jakąś uniwersalną wielokrotność (c * a) ^ 2 od r2 - gdzie c jest rozmiarem komórki, a a jest nachyleniem (jako wzrost / bieg, a nie jako kąt lub procent) - powinien załatwić sprawę.
Whuber
@ Whuber jak zawsze Twoje odpowiedzi zawierają niesamowitą ilość wiedzy !! Proszę tylko o jedno pytanie: czy to oznacza, że ​​nie można obliczyć TRI komórek znajdujących się na samych krawędziach rastra? Czy nie są otoczeni sąsiednimi komórkami dookoła?
marco
1
@marco TRI można oszacować nawet w komórkach granicznych. Jak wskazano w pytaniu, należy ją wyrazić jako średnią, a nie jako sumę, dzieląc podane tutaj wartości przez 9. W komórkach granicznych wartość „9” we wzorze i w mianowniku należy zastąpić wartością liczba wartości innych niż null w ich sąsiedztwie 3X3: 6 dla komórek krawędzi, 4 dla komórek narożnych. Siatkę takich wartości można uzyskać z sumy ogniskowej siatki wskaźników wartości pierwotnych (ma 1 na wszystkich komórkach innych niż NoData i 0 na innych miejscach). Użyj tej siatki zamiast stałej „9” we wzorach.
Whuber
3

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.

Jeffrey Evans
źródło
dzięki Jeffrey. Z jakiegoś powodu ta strona jest pusta, z wyjątkiem tytułu w Firefoksie, ponieważ w Chrome jest w porządku; może być jednym z moich rozszerzeń. Z przyjemnością informuję przynajmniej o tym, że skrypty działają bez zmian w 10.2.2 (i tak testowałem).
matt wilkie
1

-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)
użytkownik3338736
źródło
To nie jest TRI opisany w pytaniu. W rzeczywistości nie można go w ogóle traktować jako pomiaru „wytrzymałości”, ponieważ zmienia się, gdy przesuwa się jedynie pionową linię odniesienia. Na przykład, twoje TRI sąsiedztwa 3x3 o wartościach (1,2, ..., 9) będzie sqrt (9 ^ 2-1 ^ 2) = 8,9, ale dodasz 100 do wartości (co po prostu zmienia układ odniesienia bez zmiana kształtu powierzchni w ogóle) daje sqrt (109 ^ 2-101 ^ 2) = 41.
whuber
0

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 .

Don Meltz
źródło
2
TPI to zupełnie inna miara niż szorstkość. Proszę, nie idźmy dalej, używając ich zamiennie. Uważam, że Indeks pozycji topograficznej jest tradycyjnie obliczany jako [dem - focalmean (dem)].
Jeffrey Evans,