Calcul de l'indice de robustesse topographique dans ArcGIS Desktop?

Quelqu'un sait-il comment calculer l' index de robustesse topographique dans ArcGIS Desktop sans accéder à la ligne de commande ArcInfo Workstation?

"L'indice de rugosité topographique (TRI) est une mesure développée par Riley, et al. (1999) pour exprimer la quantité de différence d'élévation entre les cellules adjacentes d'une grille d'élévation numérique. Le processus calcule essentiellement la différence des valeurs d'élévation d'une cellule centrale et les huit cellules qui l'entourent immédiatement. Ensuite, il évalue chacune des huit valeurs de différence d'élévation pour les rendre toutes positives et fait la moyenne des carrés. L'indice de robustesse topographique est ensuite dérivé en prenant la racine carrée de cette moyenne, et correspond au changement d'élévation moyen entre n'importe quel point sur une grille et ses environs. " - extrait d'un aml arcscript de Jeffrey Evans

Je recommanderais de regarder en dehors d'ArcGIS) Très facile à utiliser avec le logiciel gratuit gdal: http://www.gdal.org/gdaldem.html

gdaldem TRI input_dem output_TRI_map

Ou si vous le préférez dans la saga gis: http://www.saga-gis.org/saga_modules_doc/ta_morphometry/ta_morphometry_16.html

Faisons un peu (juste un peu) d'algèbre.

Soit x la valeur du carré central; soit x_i, i = 1, .., 8 indexent les valeurs dans les carrés voisins; et soit r l'indice de rugosité topographique. Cette recette dit que r ^ 2 est égal à la somme de (x_i - x) ^ 2. Deux choses que nous pouvons calculer facilement sont (i) la somme des valeurs dans le voisinage, égale à s = ​​Sum {x_i} + x; et (ii) la somme des carrés des valeurs, égale à t = Sum {x_i ^ 2} + x ^ 2. (Ce sont des statistiques focales pour la grille d'origine et pour son carré.)

L'élargissement des carrés donne

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

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

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

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

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

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

Par exemple, considérons un quartier

1 2 3
4 5 6
7 8 9

Ici, x = 5, s = 1 + 2 + ... + 9 = 45 et t = 1 + 4 + 9 + ... + 81 = 285. Ensuite

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

et l'équivalence algébrique dit

60 = r ^ 2 = 285 + 9 * 5 ^ 2 -2 * 5 * 45 = 285 + 225 - 450 = 60, ce qui vérifie.

Le workflow est donc:

Étant donné un DEM.

• Calculer s = somme focale (sur 3 x 3 voisinages carrés) de [DEM].

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

• Calculer t = somme focale (sur 3 x 3 voisinages carrés) de [DEM2].

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

Renvoie r = Sqrt ([r2]).

Cela consiste en 9 opérations de grille au total, toutes rapides. Ils sont facilement exécutés dans la calculatrice raster (ArcGIS 9.3 et versions antérieures), la ligne de commande (toutes les versions) et Model Builder (toutes les versions).

BTW, ce n'est pas un "changement d'élévation moyen" (car les changements d'élévation peuvent être positifs et négatifs): c'est un changement d'élévation quadratique moyen. Il n'est pas égal à l '"indice de position topographique" décrit sur http://arcscripts.esri.com/details.asp?dbid=14156 , qui (selon la documentation) est égal à x - (s - x) / 8. Dans l'exemple ci-dessus, le TPI est égal à 5 ​​- (45-5) / 8 = 0 alors que le TRI, comme nous l'avons vu, est Sqrt (60).

Le TRI de Riley et al., (1999) est la racine carrée des écarts au carré sommés. Ceci est très proche de la variance non mise à l'échelle. Si vous souhaitez une implémentation de Riley TRI, veuillez suivre la méthodologie décrite par @whuber (la méthodologie fournie par @ user3338736 a généralisé la métrique au maximum dans la fenêtre et ne représente pas la variation cellule par cellule).

J'ai une variation de TRI dans notre boîte à outils ArcGIS Geomorphometry & Gradient Metrics qui est la variance d'une fenêtre spécifiée. Je trouve cela plus flexible et plus justifiable. Il existe également d'autres mesures de configuration de surface, notamment la rugosité et la dissection.

-Edit: les informations ci-dessous sont incorrectes. Veuillez consulter l'article par whuber expliquant le processus correct .....

TRI (Riley 1999) et TPI (Jenness 2002) sont similaires, mais différents.

Pour calculer TRI et TPI à l'aide d'ArcGIS 10.x ...

Étape 1: utilisez l'outil Statistiques focales pour créer 2 nouveaux jeux de données raster à partir d'un DEM.

Raster 1 "MAX") Quartier: Rectangle, Hauteur: 3, Largeur: 3, Unités: Cellule, Type de statistique: Maximum

Raster 2 "MIN") Quartier: Rectangle, Hauteur: 3, Largeur: 3, Unités: Cellule, Type de statistique: Minimum

Étape 2: utilisez la calculatrice raster pour exécuter les fonctions suivantes sur les 2 jeux de données raster que vous venez de créer.

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

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

Voici un exemple de code Python exporté à partir d'un modèle que j'ai construit pour 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)

Cela ressemble beaucoup à l'indice de position topographique, un processus que j'ai utilisé récemment pour l'un de mes projets. Il y a un ArcScript sur la page de support ESRI, une boîte à outils Topographie sur la page du centre de ressources ESRI et quelques informations supplémentaires sur le processus sur la page Entreprises Jenness .