Cálculo del índice de robustez topográfico en ArcGIS Desktop?

15

¿Alguien sabe cómo calcular el índice de rugosidad topográfica en ArcGIS Desktop sin acceso a la línea de comandos ArcInfo Workstation?

"El índice de rugosidad topográfica (TRI) es una medida desarrollada por Riley et al. (1999) para expresar la cantidad de diferencia de elevación entre las 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 lo rodean inmediatamente. Luego cuadra 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 deriva tomando la raíz cuadrada de este promedio, y corresponde al cambio promedio de elevación entre cualquier punto de una cuadrícula y su área circundante ". - de un archivo aml de Jeffrey Evans

wilkie mate
fuente
depende de la versión de ArcGIS arcscripts.esri.com/details.asp?dbid=12646 alguna discusión de foros anteriores foros.esri.com/Thread.asp?c=93&f=982&t=145448 sin verificar pero la búsqueda contenía el término jennessent.com /arcgis/surface_area.htm

Respuestas:

18

Recomendaría mirar fuera de ArcGIS) Muy fácil de usar el software gratuito de gdal: http://www.gdal.org/gdaldem.html

gdaldem TRI input_dem output_TRI_map

O si lo prefiere en saga gis: http://www.saga-gis.org/saga_modules_doc/ta_morphometry/ta_morphometry_16.html

johanvdw
fuente
11
+1 Siempre agradezco ver soluciones que no son de ArcGIS a los problemas de ArcGIS :-). Esto es una cuestión de principios, no un antagonismo hacia ArcGIS en particular. Se debe evitar quedar encerrado en una única solución de software: no solo es riesgoso profesionalmente, sino que es intelectualmente sofocante.
whuber
Sé que pedí una solución específica de arcgis, pero acepto esta por su franqueza. Las utilidades de GDAL son fáciles de adquirir e instalar, universalmente reconocidas como las mejores de su clase, y el comando para generar este producto en particular es la definición de simplicidad.
Matt Wilkie
18

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

Sea x el valor en el cuadrado central; dejemos que x_i, i = 1, .., 8 indexen 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 = Sum {x_i} + x; y (ii) la suma de cuadrados de los valores, igual a t = Sum {x_i ^ 2} + x ^ 2. (Estas son estadísticas focales para la cuadrícula original y para su cuadrado).

Expandir 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, considere 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, lo que comprueba.

El flujo de trabajo por lo tanto es:

Dado un DEM.

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

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

  • Calcule t = suma focal (más de 3 x 3 vecindades cuadradas) de [DEM2].

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

Devuelve r = Sqrt ([r2]).

Esto consiste en 9 operaciones de cuadrícula en toto , todas las cuales son rápidas. Se llevan a cabo fácilmente en la calculadora ráster (ArcGIS 9.3 y versiones anteriores), la línea de comando (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 cuadrático medio raíz. Se no igual al "índice de posición topográfica" descrito en http://arcscripts.esri.com/details.asp?dbid=14156 , que (de acuerdo con 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).

whuber
fuente
1
Gracias Bill Aprecio ver los detalles de cómo funciona una herramienta u operación. A partir de esto, cualquier persona con una inversión adecuada de tiempo y energía intelectual puede construir un nuevo aparato para llevar a cabo este trabajo utilizando las herramientas que tienen a mano. Es información como esta la que hará que GIS.se sea un servicio útil a largo plazo.
Matt Wilkie
1
+1 Gran explicación. Supongo que esto significa que una superficie empinada pero lisa podría tener un TRI más alto que una superficie plana pero con baches.
Kirk Kuykendall
1
@ Kirk Eso es correcto. Hay formas de eliminar el efecto de la pendiente local para obtener un índice de robustez "relativa" si lo desea. Aunque no he resuelto los detalles, creo que restando algunos múltiplos universales de (c * a) ^ 2 de r2, donde c es el tamaño de la celda y a es la pendiente (como subida / carrera, no como un ángulo o por ciento) - debería hacer el truco.
whuber
@whuber como siempre, ¡tus respuestas contienen una increíble cantidad de conocimiento! Solo una pregunta, por favor: ¿esto significa que no es posible calcular el TRI de las celdas ubicadas en los bordes del ráster? ¿Debido a que no están rodeados por células vecinas por todas partes?
marco
1
@marco El TRI se puede estimar incluso en las celdas límite. Como se indica en la pregunta, debe expresarse como un promedio, en lugar de una suma, dividiendo los valores que doy aquí por 9. En las celdas límite, el valor "9" en la fórmula y en el denominador debe ser reemplazado por el cantidad de valores no nulos dentro de sus vecindarios 3X3: 6 para celdas de borde, 4 para celdas de esquina. Se puede obtener una cuadrícula de dichos valores a partir de la suma focal de la cuadrícula del indicador de los valores originales (tiene 1 en todas las celdas que no son NoData y 0 en otras partes). Use esa cuadrícula en lugar de la constante "9" en las fórmulas.
whuber
3

El TRI de Riley et al. (1999) 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 celda por celda).

Tengo una variación de TRI en nuestra caja de herramientas de Geomorfometría y Gradient Metrics ArcGIS que es la varianza de una ventana específica. Esto 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.

Jeffrey Evans
fuente
gracias Jeffrey Por alguna razón, esa página está en blanco, excepto el título en Firefox, pensé que está bien en Chrome; podría ser una de mis extensiones Me complace informar al menos que los scripts funcionan sin cambios en 10.2.2 (los que probé de todos modos).
Matt Wilkie
1

-Editar: la información a continuació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 de Estadísticas Focales para hacer 2 nuevos datasets ráster a partir de un DEM.

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

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

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

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

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

Aquí hay un código Python de muestra exportado de 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)
usuario3338736
fuente
Este no es el TRI descrito en la pregunta. De hecho, no se puede considerar que mida la "robustez" en absoluto, porque cambia cuando simplemente cambia el dato vertical. Por ejemplo, su TRI de una vecindad de 3x3 con valores (1,2, ..., 9) sería sqrt (9 ^ 2-1 ^ 2) = 8.9, pero sumando 100 a los valores (lo que simplemente cambia el dato sin cambiando la forma de la superficie) da sqrt (109 ^ 2-101 ^ 2) = 41.
whuber
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, una caja de herramientas de Topografía en la página del Centro de Recursos de ESRI y algo más de información sobre el proceso en la página de Jenness Enterprises .

Don Meltz
fuente
2
TPI es una métrica muy diferente a la rugosidad. Por favor, no sigamos el camino de usarlos indistintamente. Creo que el índice de posición topográfica se calcula tradicionalmente como [dem - focalmean (dem)].
Jeffrey Evans