Menghitung Indeks Kekasaran Topografi di ArcGIS Desktop?

Adakah yang tahu cara menghitung Indeks Kekasaran Topografi di ArcGIS Desktop tanpa akses ke baris perintah ArcInfo Workstation?

"Indeks kekasaran topografi (TRI) adalah pengukuran yang dikembangkan oleh Riley, et al. (1999) untuk menyatakan jumlah perbedaan ketinggian antara sel-sel yang berdekatan dari jaringan elevasi digital. Proses ini pada dasarnya menghitung perbedaan dalam nilai ketinggian dari sel pusat. dan delapan sel yang mengelilinginya, kemudian mengkuadratkan masing-masing dari delapan nilai perbedaan ketinggian untuk membuat semuanya positif dan rata-rata kuadrat. Indeks kekasaran topografi kemudian diturunkan dengan mengambil akar kuadrat dari rata-rata ini, dan sesuai dengan perubahan elevasi rata-rata antara titik di grid dan sekitarnya. " - dari skrip aml oleh Jeffrey Evans

Saya akan merekomendasikan untuk melihat keluar ArcGIS) Sangat mudah menggunakan perangkat lunak gdal gratis: http://www.gdal.org/gdaldem.html

gdaldem TRI input_dem output_TRI_map

Atau jika Anda lebih suka di saga gis: http://www.saga-gis.org/saga_modules_doc/ta_morphometry/ta_morphometry_16.html

Mari kita lakukan aljabar kecil (hanya sedikit).

Biarkan x menjadi nilai di alun-alun pusat; misalkan x_i, i = 1, .., 8 indeks nilai-nilai di kotak tetangga; dan biarkan r menjadi indeks kekasaran topografi. Resep ini mengatakan r ^ 2 sama dengan jumlah (x_i - x) ^ 2. Dua hal yang dapat kita hitung dengan mudah adalah (i) jumlah nilai di lingkungan, sama dengan s = Jumlah {x_i} + x; dan (ii) jumlah kuadrat dari nilai-nilai, sama dengan t = Jumlah {x_i ^ 2} + x ^ 2. (Ini adalah statistik fokus untuk kisi asli dan kuadratnya.)

Memperluas kotak memberi

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

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

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

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

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

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

Misalnya, pertimbangkan lingkungan

1 2 3
4 5 6
7 8 9

Di sini, x = 5, s = 1 + 2 + ... + 9 = 45, dan t = 1 + 4 + 9 + ... + 81 = 285. Kemudian

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

dan kesetaraan aljabar mengatakan

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

The alur kerja karena itu:

Diberi DEM.

• Hitung s = Jumlah fokus (lebih dari 3 x 3 lingkungan persegi) dari [DEM].

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

• Hitung t = Jumlah fokus (lebih dari 3 x 3 lingkungan persegi) dari [DEM2].

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

Kembali r = Sqrt ([r2]).

Ini terdiri dari 9 operasi grid di toto , yang semuanya cepat. Mereka mudah dilakukan dalam kalkulator raster (ArcGIS 9.3 dan sebelumnya), baris perintah (semua versi), dan Model Builder (semua versi).

BTW, ini bukan "perubahan ketinggian rata-rata" (karena perubahan ketinggian bisa positif dan negatif): ini adalah perubahan rata-rata ketinggian akar kuadrat. Itu tidak sama dengan "indeks posisi topografi" yang dijelaskan di http://arcscripts.esri.com/details.asp?dbid=14156 , yang (menurut dokumentasi) sama dengan x - (s - x) / 8. Dalam contoh di atas, TPI sama dengan 5 - (45-5) / 8 = 0 sedangkan TRI, seperti yang kita lihat, adalah Sqrt (60).

The Riley et al., (1999) TRI adalah akar kuadrat dari deviasi kuadrat yang dijumlahkan. Ini sangat dekat dengan varian yang tidak berskala. Jika Anda menginginkan implementasi Riley's TRI maka silakan ikuti metodologi yang diuraikan oleh @whuber (metodologi yang disediakan oleh @ user3338736 menggeneralisasikan metrik secara maksimal di jendela dan tidak mewakili sel dengan variasi sel).

Saya memiliki variasi TRI di Geomorfometri & Gradien Metrik ArcGIS Toolbox kami yang merupakan varian dari jendela yang ditentukan. Saya menemukan ini lebih fleksibel dan dapat dibenarkan. Ada juga beberapa metrik konfigurasi permukaan lainnya termasuk rugosity dan diseksi.

-Edit: informasi di bawah ini salah. Silakan lihat posting oleh whuber menjelaskan proses yang benar .....

TRI (Riley 1999) dan TPI (Jenness 2002) serupa, tetapi berbeda.

Untuk menghitung TRI dan TPI menggunakan ArcGIS 10.x ...

Langkah 1: Gunakan alat Statistik Focal untuk membuat 2 set data raster baru dari DEM.

Raster 1 "MAX") Lingkungan: Persegi Panjang, Tinggi: 3, Lebar: 3, Unit: Sel, Tipe statistik: Maksimum

Raster 2 "MIN") Lingkungan: Persegi Panjang, Tinggi: 3, Lebar: 3, Unit: Sel, Tipe statistik: Minimum

Langkah 2: Gunakan Kalkulator Raster untuk melakukan fungsi-fungsi berikut pada 2 dataset raster yang baru saja Anda buat.

Untuk TRI: SquareRoot (Abs ((Kotak ("% MAX%") - Kotak ("% MIN%")))))

Untuk TPI: ("% Input DEM%" - "% MIN%") / ("% MAX%" - "% MIN%")

Berikut ini contoh kode Python yang diekspor dari model yang saya buat untuk 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)

Ini terdengar sangat mirip dengan Indeks Posisi Topografis, suatu proses yang saya gunakan baru-baru ini untuk salah satu proyek saya. Ada ArcScript di halaman dukungan ESRI, kotak alat Topografi di halaman Pusat Sumber Daya ESRI, dan beberapa info lebih lanjut tentang proses di halaman Jenness Enterprises .