Cálculo do índice de robustez topográfico no ArcGIS Desktop?

15

Alguém sabe como calcular o Índice de Robustez Topográfico no ArcGIS Desktop sem acesso à linha de comando ArcInfo Workstation?

"O índice de robustez topográfico (TRI) é uma medida desenvolvida por Riley, et al. (1999) para expressar a quantidade de diferença de elevação entre células adjacentes de uma grade de elevação digital. O processo calcula essencialmente a diferença nos valores de elevação de uma célula central e as oito células que o cercam imediatamente. Depois, esquadrinha cada um dos oito valores de diferença de elevação para torná-los positivos e calcula a média dos quadrados. O índice de robustez topográfico é então derivado da raiz quadrada dessa média e corresponde à mudança de elevação média entre qualquer ponto de uma grade e a área circundante ". - de um aml arcscript de Jeffrey Evans

Matt Wilson
fonte
depende da versão do ArcGIS arcscripts.esri.com/details.asp?dbid=12646 alguma discussão dos fóruns anteriores forums.esri.com/Thread.asp?c=93&f=982&t=145448 desmarcada, mas a pesquisa continha o termo jennessent.com /arcgis/surface_area.htm

Respostas:

18

Eu recomendaria procurar fora do ArcGIS) Muito fácil usando o software gdal gratuito: http://www.gdal.org/gdaldem.html

gdaldem TRI input_dem output_TRI_map

Ou, se preferir, em saga gis: http://www.saga-gis.org/saga_modules_doc/ta_morphometry/ta_morphometry_16.html

johanvdw
fonte
11
Eu sempre aprecio ver soluções que não sejam do ArcGIS para problemas do ArcGIS :-). Isso é uma questão de princípio, não antagonismo ao ArcGIS em particular. Deve-se evitar ficar preso a uma única solução de software: além de ser profissionalmente arriscada, é intelectualmente sufocante.
whuber
Eu sei que pedi uma solução específica para o arcgis, mas estou aceitando esta por causa de sua franqueza. Os utilitários da GDAL são fáceis de adquirir e instalar, universalmente reconhecidos como os melhores da classe, e o comando para gerar esse produto específico é a definição de simplicidade.
mate Wilkie
18

Vamos fazer um pouco (apenas um pouco) de álgebra.

Seja x o valor no quadrado central; seja x_i, i = 1, .., 8 indexe os valores nos quadrados vizinhos; e seja r o índice de robustez topográfico. Esta receita diz que r ^ 2 é igual à soma de (x_i - x) ^ 2. Duas coisas que podemos calcular facilmente são (i) a soma dos valores na vizinhança, igual a s = Sum {x_i} + x; e (ii) a soma dos quadrados dos valores, igual a t = Soma {x_i ^ 2} + x ^ 2. (Essas são estatísticas focais para a grade original e seu quadrado.)

Expandir os quadrados dá

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

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

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

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

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

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

Por exemplo, considere um bairro

1 2 3
4 5 6
7 8 9

Aqui, x = 5, s = 1 + 2 + ... + 9 = 45, e t = 1 + 4 + 9 + ... + 81 = 285. Então

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

e a equivalência algébrica diz

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

O fluxo de trabalho, portanto, é:

Dado um DEM.

  • Calcular s = soma focal (acima de 3 x 3 bairros quadrados) de [DEM].

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

  • Calcular t = soma focal (mais de 3 x 3 vizinhanças quadradas) de [DEM2].

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

Retorno r = Sqrt ([r2]).

Isso consiste em 9 operações de grade no total , todas rápidas. Eles são facilmente executados na calculadora raster (ArcGIS 9.3 e anterior), na linha de comando (todas as versões) e no Model Builder (todas as versões).

BTW, essa não é uma "alteração de elevação média" (porque as alterações de elevação podem ser positivas e negativas): é uma alteração de elevação quadrada média da raiz. É não igual ao "índice de posição topográfico" descrito no http://arcscripts.esri.com/details.asp?dbid=14156 , que (de acordo com a documentação) é igual a x - (S) - x / 8. No exemplo acima, o TPI é igual a 5 - (45-5) / 8 = 0, enquanto o TRI, como vimos, é Sqrt (60).

whuber
fonte
1
Obrigado Bill. Aprecio ver as especificidades de como uma ferramenta ou operação funciona. A partir disso, qualquer pessoa com um investimento adequado em tempo e energia intelectual pode construir um novo aparato para realizar esse trabalho usando as ferramentas que possui em mãos. São informações como essas que farão do GIS.se um serviço útil a longo prazo.
mate Wilkie
1
+1 Ótima explicação. Eu acho que isso significa que uma superfície íngreme, mas lisa, pode ter um TRI maior que uma superfície plana, mas esburacada.
21411 Kirk Kuykendall
1
@ Kirk Isso está correto. Existem maneiras de remover o efeito da inclinação local para obter um índice de robustez "relativa", se desejar. Embora eu não tenha elaborado os detalhes, acredito que subtrair algum múltiplo universal de (c * a) ^ 2 de r2 - onde c é o tamanho da célula e a é a inclinação (como subida / descida, não como um ângulo ou por cento) - deve fazer o truque.
whuber
@whuber como sempre, suas respostas contêm uma quantidade incrível de conhecimento !! Apenas uma pergunta, por favor: isso significa que não é possível calcular o TRI das células que estão localizadas nas extremidades da varredura? Por não estarem cercados por células vizinhas ao redor?
Marco
1
@marco O TRI pode ser estimado mesmo nas células de limite. Conforme indicado na pergunta, ele deve ser expresso como média, e não como soma, dividindo os valores que eu dou aqui por 9. Nas células de fronteira, o valor "9" na fórmula e no denominador precisa ser substituído pelo número de valores não nulos em suas vizinhanças 3X3: 6 para células de borda, 4 para células de canto. Uma grade desses valores pode ser obtida a partir da soma focal da grade indicadora dos valores originais (possui 1 em todas as células não NoData e 0 em outros lugares). Use essa grade no lugar da constante "9" nas fórmulas.
whuber
3

O Riley et al., (1999) TRI é a raiz quadrada dos desvios quadrados somados. Isso está muito próximo da variação não escalonada. Se você deseja uma implementação do TRI de Riley, siga a metodologia descrita por @whuber (a metodologia fornecida por @ user3338736 generalizou a métrica ao máximo na janela e não representa a célula por variação de célula).

Tenho uma variação do TRI em nossa Caixa de ferramentas do ArcGIS, Geomorfometria e métricas de gradiente, que é a variação de uma janela especificada. Acho isso mais flexível e justificável. Existem também outras métricas de configuração de superfície, incluindo rugosidade e dissecção.

Jeffrey Evans
fonte
obrigado Jeffrey. Por alguma razão, essa página está em branco, exceto pelo título no Firefox, que está correto no Chrome; pode ser uma das minhas extensões. É um prazer informar que pelo menos os scripts funcionam inalterados na versão 10.2.2 (que testei de qualquer maneira).
Matt Wilkie
1

-Editar: as informações abaixo estão incorretas. Por favor, veja o post do whuber explicando o processo correto .....

TRI (Riley 1999) e TPI (Jenness 2002) são semelhantes, mas diferentes.

Para calcular TRI e TPI usando o ArcGIS 10.x ...

Etapa 1: use a ferramenta Estatísticas focais para criar 2 novos conjuntos de dados raster a partir de um DEM.

Raster 1 "MAX") Bairro: Retângulo, Altura: 3, Largura: 3, Unidades: Célula, Tipo de estatística: Máximo

Raster 2 "MIN") Bairro: Retângulo, Altura: 3, Largura: 3, Unidades: Célula, Tipo de estatística: Mínimo

Etapa 2: use a Calculadora de varredura para executar as seguintes funções nos 2 conjuntos de dados de varredura que você acabou de criar.

Para TRI: SquareRoot (Abs ((Quadrado ("% MAX%") - Quadrado ("% MIN%")))))

Para TPI: ("% DEM% de entrada%" - "% MIN%") / ("% MAX%" - "% MIN%")

Aqui está um exemplo de código Python exportado de um modelo que criei para o 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
fonte
Este não é o TRI descrito na pergunta. De fato, ele não pode ser considerado como uma medida da "robustez", porque muda quando você apenas muda o dado vertical. Por exemplo, seu TRI de uma vizinhança 3x3 com valores (1,2, ..., 9) seria sqrt (9 ^ 2-1 ^ 2) = 8,9, mas adicionando 100 aos valores (o que apenas altera o dado sem alterar a forma da superfície) resulta em sqrt (109 ^ 2-101 ^ 2) = 41.
whuber
0

Isso se parece muito com o Índice de Posição Topográfica, um processo que usei recentemente para um dos meus projetos. Há um ArcScript na página de suporte da ESRI, uma caixa de ferramentas Topografia na página do ESRI Resource Center e mais algumas informações sobre o processo na página Jenness Enterprises .

Don Meltz
fonte
2
O TPI é uma métrica muito diferente da rugosidade. Por favor, não vamos usar os mesmos de forma intercambiável. Eu acredito que o Índice de Posição Topográfica é tradicionalmente calculado como [dem - focalmean (dem)].
Jeffrey Evans