¿Cómo calcular el índice de resistencia topográfica en ArcGIS Desktop?

15

¿Alguien sabe cómo calcular Índice de resistencia topográfica en ArcGIS Desktop sin acceso a la línea de comandos de la estación de trabajo ArcInfo?

"El índice de rugosidad topográfica (TRI) es una medida desarrollada por Riley, y otros (1999) para expresar la cantidad de diferencia de elevación entre celdas adyacentes de una cuadrícula de elevación digital. El proceso esencialmente calcula la diferencia en los valores de elevación de una celda central y las ocho celdas que la rodean inmediatamente. Luego, ajusta cada uno de los ocho valores de diferencia de elevación para hacerlos todos positivos y promedia los cuadrados. El índice de rugosidad topográfica se obtiene al tomar la raíz cuadrada de este promedio. y corresponde al cambio de elevación promedio entre cualquier punto en una cuadrícula y su área circundante ". - desde un aml arcscript de Jeffrey Evans

    
pregunta matt wilkie 11.02.2011 - 01:43

5 respuestas

18

Recomendaría mirar fuera de ArcGIS) Muy fácil de usar el software gratuito gdal: enlace

gdaldem TRI input_dem output_TRI_map

O si lo prefieres en la saga gis: enlace

    
respondido por el johanvdw 11.02.2011 - 08:46
18

Hagamos un poco (solo un poco) de álgebra.

Sea x el valor en el cuadrado central; hagamos que x_i, i = 1, .., 8 indexe los valores en los cuadrados vecinos; y sea r el índice de rugosidad topográfica. Esta receta dice que r ^ 2 es igual a la suma de (x_i - x) ^ 2. Dos cosas que podemos calcular fácilmente son (i) la suma de los valores en la vecindad, igual a s = Suma {x_i} + x; y (ii) la suma de los cuadrados de los valores, igual a t = Suma {x_i ^ 2} + x ^ 2. (Estas son estadísticas focales para la cuadrícula original y para su cuadrado.)

La expansión de los cuadrados da

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 .

Por ejemplo, considera un vecindario

1 2 3
4 5 6
7 8 9

Aquí, x = 5, s = 1 + 2 + ... + 9 = 45, y t = 1 + 4 + 9 + ... + 81 = 285. Entonces

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

y la equivalencia algebraica dice

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

El flujo de trabajo por lo tanto es:

Dado un DEM.

  • Calcular s = Suma focal (más de 3 x 3 vecindades cuadradas) de [DEM].

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

  • Calcular t = suma focal (sobre 3 x 3 vecindades cuadradas) de [DEM2].

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

Volver r = Sqrt ([r2]).

Esto consiste en 9 operaciones de red in toto , todas rápidas. Se realizan fácilmente en la calculadora ráster (ArcGIS 9.3 y anteriores), la línea de comandos (todas las versiones) y Model Builder (todas las versiones).

Por cierto, esto no es un "cambio de elevación promedio" (porque los cambios de elevación pueden ser positivos y negativos): es un cambio de elevación del cuadrado medio de la raíz. Es no igual al "índice de posición topográfica" descrito en enlace , que (según la documentación) es igual a x - (s - x) / 8. En el ejemplo anterior, el TPI es igual a 5 - (45-5) / 8 = 0 mientras que el TRI, como vimos, es Sqrt (60).

    
respondido por el whuber 11.02.2011 - 05:07
3

Riley et al., (1999) TRI es la raíz cuadrada de las desviaciones cuadradas sumadas. Esto está muy cerca de la varianza sin escala. Si desea una implementación del TRI de Riley, siga la metodología descrita por @whuber (la metodología proporcionada por @ user3338736 generalizó la métrica al máximo en la ventana y no representa la variación de celda por celda).

Tengo una variación de TRI en nuestro Geomorfometría & Gradient Metrics ArcGIS Toolbox que es la varianza de una ventana específica. Me parece más flexible y justificable. También hay algunas otras métricas de configuración de superficie que incluyen rugosidad y disección.

    
respondido por el Jeffrey Evans 10.06.2014 - 18:16
1

-Editar: la siguiente información es incorrecta. Por favor vea la publicación de whuber explicando el proceso correcto .....

TRI (Riley 1999) y TPI (Jenness 2002) son similares, pero diferentes.

Para calcular TRI y TPI utilizando ArcGIS 10.x ...

  

Paso 1: use la herramienta Focal Statistics para crear 2 nuevos datasets ráster desde un DEM.

     

Raster 1 "MAX") Barrio: Rectángulo, Altura: 3, Ancho: 3, Unidades:   Celda, Tipo de estadísticas: Máximo

     

Raster 2 "MIN") Barrio: Rectángulo, Altura: 3, Ancho: 3, Unidades:   Celda, Tipo de estadísticas: Mínimo

     

Paso 2: Usa la Calculadora ráster para realizar las siguientes funciones en los 2 conjuntos de datos ráster que acabas de crear.

     

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

     

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

Aquí se muestra un código de Python de muestra exportado desde un modelo que construí para 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)
    
respondido por el user3338736 06.03.2014 - 18:09
0

Esto se parece mucho al índice de posición topográfica, un proceso que utilicé recientemente para uno de mis proyectos. Hay un ArcScript en la página de soporte de ESRI, un Caja de herramientas de topografía en la página del Centro de recursos de ESRI, y más información sobre el proceso en Jenness Enterprises .

    
respondido por el Don Meltz 11.02.2011 - 04:07

Lea otras preguntas en las etiquetas