Вычислить индекс топографической прочности в ArcGIS Desktop?

15

Кто-нибудь знает, как рассчитать индекс топографической стойкости в ArcGIS Desktop без доступа к командной строке ArcInfo Workstation?

«Индекс топографической прочности (TRI) - это измерение, разработанное Riley и др. (1999) для выражения величины разности высот между соседними ячейками цифровой сетки высот. Этот процесс по существу вычисляет разность значений высот из центральной ячейки и восемь ячеек, непосредственно окружающих его. Затем он возводит в квадрат каждое из восьми значений разности высот, чтобы сделать их все положительными, и усредняет квадраты. Индекс топографической прочности затем получается путем извлечения квадратного корня из этого среднего и соответствует среднему изменению высот между любой точкой на сетке и ее окружением. "- из рукописи Амла Джеффри Эванса

Мэтт Уилки
источник
зависит от версии ArcGIS arcscripts.esri.com/details.asp?dbid=12646 некоторого обсуждения предыдущих форумов forums.esri.com/Thread.asp?c=93&f=982&t=145448 бесконтрольно , но поиск содержал термин jennessent.com /arcgis/surface_area.htm

Ответы:

18

Я бы порекомендовал заглянуть за пределы ArcGIS) Очень просто с помощью бесплатного программного обеспечения GDAL: http://www.gdal.org/gdaldem.html

gdaldem TRI input_dem output_TRI_map

Или, если вы предпочитаете его в сагах: http://www.saga-gis.org/saga_modules_doc/ta_morphometry/ta_morphometry_16.html

johanvdw
источник
11
+1 Мне всегда приятно видеть решения ArcGIS, не относящиеся к ArcGIS :-). Это принципиальный вопрос, а не антагонизм к ArcGIS в частности. Следует избегать привязки к одному программному решению: оно не только профессионально рискованно, но и душно интеллектуально.
whuber
Я знаю, что просил конкретное решение arcgis, но я принимаю это из-за его прямоты. Утилиты GDAL легко приобрести и установить, они широко признаны лучшими в своем классе, и команда для создания этого конкретного продукта - определение простоты.
Matt Wilkie
18

Давайте сделаем немного (просто немного) алгебры.

Пусть х будет значением в центральной площади; пусть x_i, i = 1, .., 8 индексируют значения в соседних квадратах; и пусть r будет индекс топографической прочности. Этот рецепт говорит, что r ^ 2 равно сумме (x_i - x) ^ 2. Мы можем легко вычислить две вещи: (i) сумму значений в окрестности, равную s = Sum {x_i} + x; и (ii) сумма квадратов значений, равная t = Sum {x_i ^ 2} + x ^ 2. (Это фокусная статистика для исходной сетки и для ее квадрата.)

Расширение квадратов дает

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

= Сумма {x_i ^ 2 + x ^ 2 - 2 * x * x_i}

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

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

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

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

Например, рассмотрим окрестности

1 2 3
4 5 6
7 8 9

Здесь x = 5, s = 1 + 2 + ... + 9 = 45 и t = 1 + 4 + 9 + ... + 81 = 285. Тогда

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

и алгебраическая эквивалентность говорит

60 = r ^ 2 = 285 + 9 * 5 ^ 2 -2 * 5 * 45 = 285 + 225 - 450 = 60, что проверяет.

рабочий , следовательно , является:

Учитывая DEM.

  • Вычислите s = Фокальная сумма (более 3 x 3 квадратных окрестностей) из [DEM].

  • Вычислить DEM2 = [DEM] * [DEM].

  • Вычислить t = Фокальная сумма (более 3 x 3 квадратных окрестностей) из [DEM2].

  • Вычислить r2 = [t] + 9 * [DEM2] - 2 * [DEM] * [s].

Возвращение r = Sqrt ([r2]).

Это состоит из 9 операций сетки в целом , все из которых являются быстрыми. Они легко выполняются в калькуляторе растра (ArcGIS 9.3 и более ранних версиях), командной строке (все версии) и построителе моделей (все версии).

Кстати, это не «среднее изменение высоты» (потому что изменение высоты может быть положительным и отрицательным): это среднеквадратичное изменение высоты. Он не равен «индексу топографической позиции», описанному по адресу http://arcscripts.esri.com/details.asp?dbid=14156 , который (согласно документации) равен x - (s - x) / 8. В приведенном выше примере TPI равен 5 - (45-5) / 8 = 0, тогда как TRI, как мы видели, равен Sqrt (60).

Whuber
источник
1
Спасибо, Билл. Я ценю понимание специфики работы инструмента или операции. Исходя из этого, любой, у кого есть подходящее вложение времени и интеллектуальной энергии, может создать новый аппарат для выполнения этой работы, используя имеющиеся у него инструменты. Именно такая информация сделает GIS.se полезным сервисом в долгосрочной перспективе.
Мэтт Уилки
1
+1 Отличное объяснение. Я предполагаю, что это означает, что крутая, но гладкая поверхность может иметь более высокое значение TRI, чем плоская, но неровная поверхность.
Кирк Куйкендалл
1
@ Кирк Это правильно. Существуют способы устранить влияние локального уклона, чтобы получить индекс «относительной» прочности, если хотите. Хотя я не проработаны детали, я считаю , что вычитая некоторое универсальное кратное (с * а) ^ 2 из r2 - где с является cellsize и является наклон (как подъем / бега, а не как угол или процентов) - должен сделать свое дело.
whuber
@whuber как всегда ваши ответы содержат удивительное количество знаний! Только один вопрос, пожалуйста: означает ли это, что невозможно рассчитать TRI ячеек, которые расположены на самых краях растра? Из-за того, что они не окружены соседними камерами?
Марко
1
@marco TRI можно оценить даже в граничных ячейках. Как указано в вопросе, оно должно быть выражено как среднее, а не сумма, путем деления значений, которые я даю здесь, на 9. В граничных ячейках значение «9» в формуле и в знаменателе необходимо заменить на число ненулевых значений в их окрестностях 3X3: 6 для краевых ячеек, 4 для угловых ячеек. Сетка таких значений может быть получена из фокальной суммы индикаторной сетки исходных значений (она имеет 1 во всех ячейках, отличных от NoData, и 0 в других местах). Используйте эту сетку вместо константы «9» в формулах.
whuber
3

Riley et al., (1999) TRI является квадратным корнем из суммированных квадратов отклонений. Это очень близко к немасштабированной дисперсии. Если вы хотите реализовать TRI Riley, пожалуйста, следуйте методологии, изложенной @whuber (методология, предоставленная @ user3338736, обобщает метрику до максимума в окне и не представляет ячейку по вариации ячейки).

У меня есть вариант TRI в нашем ArcGIS Toolbox для геоморфометрии и градиента, который представляет собой дисперсию указанного окна. Я считаю это более гибким и оправданным. Есть также некоторые другие метрики конфигурации поверхности, включая морщинистость и рассечение.

Джеффри Эванс
источник
спасибо Джеффри По какой-то причине эта страница пуста, за исключением заголовка в Firefox, но в Chrome все в порядке; может быть одним из моих расширений. Я рад сообщить, по крайней мере, что скрипты работают без изменений в 10.2.2 (те, которые я тестировал в любом случае).
Мэтт Вилки
1

Редактировать: приведенная ниже информация неверна. Пожалуйста, смотрите пост whuber, объясняющий правильный процесс .....

TRI (Riley 1999) и TPI (Jenness 2002) похожи, но различны.

Для расчета TRI и TPI с помощью ArcGIS 10.x ...

Шаг 1: Используйте инструмент Focal Statistics для создания 2 новых наборов растровых данных из матрицы высот.

Растр 1 "MAX") Соседство: Прямоугольник, Высота: 3, Ширина: 3, Единицы: Ячейка, Тип статистики: Максимум

Растр 2 "MIN") Соседство: Прямоугольник, Высота: 3, Ширина: 3, Единицы: Ячейка, Тип статистики: Минимум

Шаг 2: Используйте калькулятор растров для выполнения следующих функций на 2 наборах растровых данных, которые вы только что создали.

Для TRI: SquareRoot (Abs ((Квадрат ("% MAX%") - Квадрат ("% MIN%"))))

Для TPI: («% Input DEM%» - «% MIN%») / («% MAX%» - «% MIN%»)

Вот пример кода Python, экспортированного из модели, которую я построил для 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)
user3338736
источник
Это не TRI, описанный в вопросе. На самом деле, его вообще нельзя считать измерением «прочности», потому что оно меняется, когда вы просто сдвигаете вертикальную точку отсчета. Например, ваш TRI соседства 3x3 со значениями (1,2, ..., 9) будет иметь значение sqrt (9 ^ 2-1 ^ 2) = 8,9, но при добавлении 100 к значениям (что просто изменяет данные без изменение формы поверхности вообще) дает sqrt (109 ^ 2-101 ^ 2) = 41.
whuber
0

Это очень похоже на Индекс топографической позиции, процесс, который я недавно использовал для одного из моих проектов. На странице поддержки ESRI есть ArcScript , набор инструментов Топография на странице Центра ресурсов ESRI и дополнительная информация о процессе на странице Jenness Enterprises .

Дон мельц
источник
2
TPI - это совсем другой показатель, чем шероховатость. Пожалуйста, давайте не будем идти по пути взаимозаменяемости их использования. Я считаю, что Индекс топографической позиции традиционно рассчитывается как [dem - focalmean (dem)].
Джеффри Эванс