Hacking the topographic ruggedness index

11 The topographic ruggedness index (TRI) is widely adopted for the analysis of digital elevation 12 models, providing information on local surface spatial variability. In this work, the TRI is 13 interpreted according to a geostatistical perspective, highlighting its main characteristics and 14 drawbacks. TRI can be interpreted as an omnidirectional short-range spatial variability index, 15 computed according to a pixel centered perspective. The simplicity and interpretability of 16 the index, free from user-dependent selections, promoted its implementation in several 17 software environments and its application in a wide set of case studies. However, the index 18 has several drawbacks for its application in earth sciences, such as a strong dependency on 19 local slope (it is basically an average adjacent neighbor slope algorithm) and the selection of 20 different lag distances in the computation of spatial variability along the main directions and 21 the diagonal ones. We propose a new metric radial roughness (RRI) in order to solve the 22 main drawbacks of TRI but maintaining its main philosophy (i.e., pixel centered perspective 23 and simplicity of the algorithm). The new index corrects for the differences in lag distances 24 and resolves the dependency on trend using increments of order 2. The code of the index, 25 implemented in R statistical language, and test data are provided with the paper 26 (https://doi.org/10.5281/zenodo.7132160) to promote its implementation in other software 27 environments.


Introduction
The analysis of surface roughness from digital elevation models (DEMs), analogously to the analysis of image texture from imagery (Haralick et al., 1973), represents a key element in many geoenvironmental, geoengineering and remote sensing applications (e.g., McKean and Roering, 2004;Glenn et al., 2006;Wilson et al., 2007;Woodcock et al., 1988, Różycka et al., 2017).But the possible applications can be broader, including for example ecological applications (e.g., Nelleman and Thomsen, 1994) and the study of distribution of human residential development (Vukomanovic and Orr, 2014).It is not our intention to provide an epistemological discussion on the concept of roughness in earth sciences and on the related terminology.In the context of earth sciences (Shepard et al. 2001;Smith, 2014) and geomorphometry (Pike, 2000;Wilson and Gallant, 2000) there is not a consensus on a definition of roughness, and multiple related terms are frequently adopted, such as "fabric", "ruggedness", "rugosity", "texture", "waviness" and others.Here, surface roughness (further referred for brevity as roughness) is intended as synonym of surface texture and hence it is a complex and multiscale characteristic of a surface' spatial variability.Accordingly, this interpretation deviates from the definition provided in surface metrology (e.g., Leach, 2013), where roughness is referred to short-wavelength surface texture features.When the analysis of roughness is performed on a digital elevation model (DEM), represented as a regular grid (i.e., a raster image), there is a complete algorithmic convergence between image texture and surface texture (Haralick et al., 1973;Woodcok et al. 1988;Atkinson and Lewis, 2000;Ojala et al., 2002;Lucier and Stein 2005;Garrigues et al., 2006;Balaguer et al., 2010;Trevisani et al. 2012).Moreover, the same concepts can be adopted to real 3D surfaces as well as to elevation points clouds (Pollyea and Fairley, 2011), given the proper adaptation of algorithms.It should be also highlighted that with usual DEMs, which are a 2.5D representation of solid earth surface (Burrough and McDonnel, 1998) the roughness computed from them is an apparent roughness, given that this typology of DEM represents the projection of a 3D surface on an horizontal plane.
The topographic ruggedness index (TRI, Riley et al., 1999) is a widely adopted measure of short-range roughness.It was initially developed in the context of ecological analysis because many biosphererelated processes are influenced by surface roughness (Nellemann and Thomsen, 1994;Jenness, 2004;Sappington, 2007;Wilson et al., 2007;Hagerty et al., 2011).Even though more statistically robust and informative roughness indexes can be considered (e.g., Herzfeld and Higginson, 1996;Trevisani and Rocca, 2015;Trevisani et al. 2023), the approach underlying TRI is still valuable and it is worth trying to improve its drawbacks while keeping its appealing features.TRI is particularly interesting in earth sciences and ecology applications because of its pixel centered perspective in evaluating roughness.Another characteristic that made popular this index is the extreme simplicity of the algorithm (Wilson et al., 2007) that fostered its implementation, efficiently and easily, in multiple software packages, e.g.GRASS GIS (GRASS, 2023), MICRODEM (USNA, 2023), SAGA GIS, Terra package for R (RSPATIAL, 2023), Whitebox (Whitebox Geospatial, 2023), ArcMap (ESRI, 2023), and others.However, TRI has some severe limitations, often leading to a misuse of the index.A first limitation is the strong dependency of the index to local slope, unlike other indexes (Woodcock, 1977;Guth 1999 and2001;Wilson and Gallant, 2000;Grohmann et al., 2011;Trevisani et al. 2023).Other issues are related to the usage of different pixel distances along the main directions (NS and EW) with respect to the diagonal ones and the use of only 8 samples in the estimation of the index, making it highly sensitive to noise.
Accordingly, the main aim of this technical note is to provide a new implementation of the TRI, that maintains the main philosophy behind it (i.e., simplicity and pixel centered perspective) while improving its main drawbacks.The presented improvements are based on ongoing research on geostatistical-based approaches for roughness analysis, that led to the recent implementation of an efficient and informative geostatistical algorithm (Trevisani and Rocca, 2015;Trevisani et al., 2023;Trevisani, 2023), capable of highlighting specific aspects of roughness.Moreover, the reanalysis of TRI from the geostatistical perspective permits highlighting its key features and finding a clear pathway for its improvement.
The implementation provided here refers to the original TRI implementation relying on a 3x3 kernel; alternative extensions devised for multiscale analysis are not considered.Some approaches, for example, generalize TRI considering larger kernels with weights or adopt smoothing approaches prior to the computation (Wilson et al. 2007).These variations are not considered here because they increase the complexity of the approach, making less convenient its application with respect to other more flexible and controllable approaches based on geostatistics (e.g., Garrigues et al., 2006;Balaguer et al., 2010;Trevisani and Rocca 2015;Trevisani et al, 2023).Moreover, the improvements proposed here can still be adopted for multiscale analysis without modifications, using multilevel upscaling of the original DEM (e.g., Wilson et al., 2007;Lindsay et al. 2019;Newman et al., 2022).It should be also noted that the dependence on slope can be filtered out by calculating TRI on a residual DEM (Guisan et al., 1999;Wilson and Gallant, 2000;Ilich et al., 2021), also known as topographic position index (TPI).TPI also suffers because it does not consider the different calculation distances along the diagonals.However, as discussed in Trevisani et al. (2023), there are different approaches and different calculation parameters that can be selected for deriving the residual DEM, generating some subjectivity in its derivation.The solution adopted here is unequivocal and the modified index, as the original one, doesn't require user-based choices, apart from the resolution of the input DEM.

Methodology
2.1 TRI according to geostatistics TRI has clear connections with geostatistical spatial variability estimators; in fact, differences between pixel values separated by a given lag distance (further referred as directional differences, DDs) are the building blocks of estimators of spatial variability (Isaaks and Srivastava, 1989;Goovaerts, 1997;Chilès and Delfiner 2012), such as the variogram (eq.1), the madogram (eq.2) and the robust version based on the median of absolute differences (MAD, Trevisani and Rocca, 2015). (1) where ()  = (  ) − (  + ).In Eqs. ( 1)-( 2),  is the separation vector (lag) between two locations (u), () is the value of the variable of interest in the location u (e.g., elevation, residual elevation, band intensity, etc.), and () is the number of point pairs with a separation vector  found in the search window considered.Hence, ()  = (  ) − (  + ) represents a given DD.
The main difference between TRI and geostatistical approaches relies on the way in which the DDs are sampled from the local spatial domain considered (e.g., search window or kernel, the two terms are used interchangeably).In the original formulation, TRI is calculated with a 3x3 kernel, considering the eight differences between the central pixel and the external ones (figure 1 left).Accordingly, it represents a radial measure of spatial variability.In a usual geostatistical approach (figure 1, right), with a 3x3 window, the set of included DDs in the estimation would be different and much larger.Moreover, with the geostatistical approaches adapted to geomorphometry (Trevisani and Rocca 2015;Trevisani et al., 2023) the set of samples is still larger and includes a correction for diagonal distances.
In the examples provided in figure 1, we don't consider the issue of the different lag distances between main directions and diagonal ones; this topic will be discussed later.Moreover, it is worth noting that the directional differences of TRI if divided by the lag distances are essentially (discrete) directional derivatives radiating from the center, and hence related to the divergence.In the original formulation of Riley et al. (1999), TRI is defined as the square root of the sum of the squared DDs; if normalized by the number of DDs (i.e., divided by 8), it is the square root of the mean of the squared DDs (eq.3); accordingly, it is analogous to the square root of the variogram (eq.1), apart from the division by 2. Ultimately, it provides an omnidirectional measure of spatial variability for a lag distance of 1.207 px, where px is the non-SI unit pixel.
An improved version, as implemented in Terra package and in most software, computes the index by means of the absolute value of DDs (e.q. 4, Wilson et al., 2007), given that squared DDs are very sensitive to data contamination.Accordingly, it is like an omnidirectional madogram (eq.2).
The last version is convenient from the interpretative point of view because it provides the mean absolute difference in elevation observed by the central pixel respect to its neighbors.However, the average lag distance is 1.207 px, which is a little unpractical from the interpretative point of view; moreover, the different lag distances can generate various biases with anisotropic features depending on their orientation, wavelengths, and spatial variability characteristics.
The evident similarities with geostatistics explain very clearly the dependency of TRI on slope; in fact, in geostatistics, the variables from which to calculate the spatial variability/continuity indexes should be at least in conditions of intrinsic stationarity (Isaaks and Srivastava, 1989).This means that in the considered spatial domain (i.e., a 3x3 window for TRI) there should not be a trend, representing longrange variations.If a trend is present, the estimates of spatial variability would be biased.This is evident with TRI; if this index is applied on a perfectly smooth but sloping face, it will return a nonzero result.In ecological applications the influence of slope can be the intended result; for example, the main factor controlling animal movement could be the overall spatial variability, considering both short-range (the residual) and long-range (the trend) components.However, for most earth sciences related applications, such as for deriving proxies of flow impedance and for geomorphological interpretation, the short-range roughness is the object of the study.For this kind of aim, TRI should not be applied directly to a DEM, otherwise it will represent a proxy of slope more than a measure of short-range roughness (figure 2).For the 2 m resolution DEM represented in figure 2 the correlation between TRI and slope is 0.989.As reported above, one adopted solution the computation of TRI on a detrended surface (Wilson et al. 2007;Ilich et al., 2021); however, the detrending can be performed in different ways and this induces higher complexity and user choices (Trevisani et al., 2023).We present a solution that resolves these issues maintaining the simplicity and the pixel centered perspective of the original algorithm.The current implementation in R is devised for DEMs and imagery represented in a projected coordinate system.For working in geographical coordinates, the algorithm has been adapted in MICRODEM according to the latitude/longitude distance corrections (Guth and Kane, 2021).

The modified algorithms
Accordingly, we present two main modifications of TRI original formulation (eq.4) that, while maintaining the simplicity of the original approach, remove the effect of slope and correct for different lag distances.These modifications lead to the definition of a short-range radial roughness index (RRI).The algorithms are implemented in R using the "Terra" library (RSPATIAL; 2023).However, the algorithms can be implemented easily in any GIS and image analysis environment; the code and test data provided should serve as the basis for implementing it in any environment.For example, these have been also implemented in the open source software MICRODEM, written in the Delphi (object Pascal) language (Guth, 2023).
In terms of software implementation, referring to the layout of pixels indexing for 3x3 kernel in figure 3, the pseudocode of conventional TRI can be expressed as: (5) that is the mean of absolute difference between the central pixel and its eight surroundings.
Figure 3. Pixels indexing layout for 3x3 kernels adopted for representing the pseudocode of TRI.

TRIK2: removing the influence of local slope
The impact of local slope/trend can be filtered out with the same approach based on increments of order-k used in geostatistics (Chilès and Delfiner, 2012) and recently implemented in the algorithm for roughness analysis of Trevisani et al. (2023).In fact, the detrending can be avoided by exploiting the capability of increments (in this case differences) of order K to filter out a polynomial trend of order k-1 (Chilès and Delfiner, 2012), the core idea behind generalized covariances in geostatistics.
For example (figure 4), considering the differences of differences (increments of order 2), it is possible to filter out a trend of order 1, i.e., a planar trend in 2D.The modeling of the trend with a simple planar trend is a reasonable assumption when considering short lag distances.For longer distances quadratic or cubic surfaces should be considered (Burrough and McDonnel, 1998) and hence higher order increments.

Figure 4. Example (1D) illustrating how the differences of order 2 can be computed with one convolution for lags of 1 px (in the formula T indicates transposition).
An analytical example for computing the differences of order 2 centered on the pixel  2 is shown.
Accordingly, the new TRIK2 estimator (in which "K2" is reminiscent of increments of order 2), referring to the pixel indexing layout of figure 5, is derived with the following pseudocode: A graphical representation of the approach is reported in figure 6.It should be noted that TRIK2 can be computed also with a 3x3 kernel (figure 6, orange and red arrows in the central pixel); however, only 4 DDs of order 2 are used (directions N-S, NE-SW, E-W and SE-NW), leading to a potentially noisy metric.With the 5x5 kernel, 12 DDs of order 2 are used, leading to larger number of samples to compute the index and hence more statistical stability respect to the original formulation.As done for TRI, it is worth noting that, given the layout of figure 6, dividing the DDs of order 1 and 2 by the lag distance one obtains the (discrete) directional derivatives of order 2 radiating from the central pixel, and hence there is a relation with the discrete Laplacian operator.Therefore, like the Laplacian operator, TRIK2 is able to detect sharp edges and, in general, high frequency (i.e.short range) features, which is the case of short-range roughness.Clearly, TRIK2, like the Laplacian, is affected by highfrequency noise.On the one hand, it is clear that the measurement of a quality such as roughness requires that noise of comparable frequency, i.e. 1 px -1 , be negligible.On the other hand, the fact that TRIK2 comes from the mean of 12 values, each of which corresponds to a 1D discrete Laplacian operator, as shown in Eq. ( 6), reduces the effect of noise.
Figure 6.The calculation is performed in one convolution, as in figure 4 for the 1D case; however, for illustrative purposes it can be interpreted as two step-convolution: step 1 (left) computation DDs of order 1; step 2, computation of DDs of order 2 from the differences (green squares) computed in step 1. Blue arrows represent directional differences of order 1 with lags of 1 px and the light blue with a lag of √2 pixels.Red arrows represent directional differences of order 2 with lags of 1 px and the orange with a lag of √2 px.
This solution resolves the algorithmic dependency on local slope as shown in figure 7, where it is applied on the same high-resolution DEM of figure 2. For this DEM, the correlation between TRIk2 and slope is 0.753.

Radial roughness index
The potential bias related to the different lag distances between main directions and the diagonal ones can be reduced with the same approach adopted by Trevisani et al. 2015, based on bilinear interpolation.In particular (figure 8) bilinear interpolation is adopted to derive an elevation value at a distance of 1 px along the diagonals.Accordingly, using the weights of bilinear interpolation is possible to modify the code for TRIK2 so as to calculate the differences at the same 1 px distance in all directions (figure 8).The modified algorithm is named radial roughness index (RRI) and its  From the visual perspective RRI (figure 9 top) apparently produces the same patterns to TRIk2 (figure 7); however, in correspondence of anisotropic features there could be marked differences (figure 9 bottom) and most of the time RRI will compute lower values.In fact, for the DEM considered, there are only 0.2% of cases in which the differences between TRIk2 and RRI are lower than -0.02 m.In fact, RRI will provide lower values than TRIk2 in presence of anisotropic features with wavelengths longer than 2√2 px in the direction of maximum spatial variability, depending on their orientation; however, with anisotropic features whose wavelength is shorter than 2√2 px, depending on their orientation, RRI could return higher variability than TRIK2.

Discussion and conclusion
The solution provided is simple and easy to implement, and this should be preferred to the wellknown TRI when the focus is on short-range radial roughness without the effect of slope.RRI minimizes the bias in spatial variability estimation and improves the interpretation of the index, that now is referred to a 1 px lag.The applications of RRI are analogous to those of TRI, which is a measure of topographic heterogeneity.However, RRI highlights the effects of short-wavelength terrain undulations, whereas TRI depends on both local slope and short-wavelength undulations.For the DEM shown in the figures here, RRI decreases the correlation between TRI and slope from 0.989 to a correlation of 0.753 between RRI and slope; the new name emphasizes that it is a new and independent parameter and not just another measure of slope.
From the set of twelve DDs of order 2 (corrected or not for diagonals) one can be tempted to consider other estimators with respect to the mean, such as the median and the minimum.The median could reduce the impact of hotspots and the minimum can be interesting from the ecological perspective, given that very often the governing driver is the lowest impedance found in the surrounding directions.However, the low number of samples limits this kind of application.For example, the adoption of the median as estimator would be analogous to MAD based estimators (Trevisani and Rocca 2015;Trevisani et al., 2023) including MADk2, the version using increments of order 2.However, for MAD and MADk2 (figure 10 bottom right) a much larger set of DDs is generally used for deriving the indexes of roughness (omnidirectional roughness and anisotropy).With MAD, for statistical reasons, it is suggested to use at least a 5 x 5 window, permitting to select 25 DDs in each of the 4 directions (main ones and diagonals); ultimately, for the isotropic short-range roughness 100 DDs are used.MAD approach is more robust to non-stationarity and abrupt transitions (e.g., figure 10).Moreover, this approach permits computing other roughness indexes such as anisotropy and the possibility to select different lag distances.
RRI can be considered a valid option as an omnidirectional short-range roughness metric, when the pixel centered perspective is significant for the study at hand; otherwise, the geostatistical estimators provide a more flexible alternative, with a sound theoretical framework.The implemented code and test files (DEM and outputs presented in this note) are available (Trevisani, 2023), promoting the implementation in readers' preferred software.
Finally, as far as computational efficiency is concerned, calculating RRI for the DEM considered here (716x943 px) takes 3.46 s with a 3.2 MHz 8 th Gen.Intel Core i7 CPU with 16 GB RAM.It is therefore a program that can run on any PC currently available.Moreover, the current implementation can be further improved for computing performance.(Trevisani et al., 2023).The latter has been computed with a lag of 1 px, with a 3x3 (bottom left) and 5x5 (bottom right) square search window.MAD is conceived to be more robust to non-stationarity and abrupt morphological transitions, permitting a sharper representation of roughness.

Figure 1 .
Figure 1.Example of selection of DDs for TRI (left) and for an hypothetical geostatistical estimator of omnidirectional variability (right).Blue arrows represent directional differences with lags of 1 pixel and the light blue with a lag of √2 pixels.

Figure 2 .
Figure 2. Example of contrasting local surface roughness for an Alpine glacial environment (top left: DEM; top right: hillshade).Calculation of conventional TRI (bottom right) highlighting the strong dependence on slope (bottom left).The 2 m resolution DEM has been derived from an airborne Lidar survey (Trentino, Italy).Color scales for slope and TRI are histogram equalized.
Figure 5. Pixel indexing layout for a 5x5 kernel adopted for the pseudocode of TRIK2

Figure 8 .
Figure 8. Example of computation of DDs of order 2 correcting for the different lag distance between the DD in NS direction (1 px) and NE-SW direction (√2 px).

Figure 9 .
Figure 9. Short-range radial roughness computed on the high-resolution DEM with the RRI algorithm (top).The differences with TRIK2 (bottom) can be relevant in presence of anisotropic features.

Figure 10 .
Figure10.For a detail of the test area of figure2, comparison of short-range roughness computed with RRI and with MADk2(Trevisani et al., 2023).The latter has been computed with a lag of 1 px, with a 3x3 (bottom left) and 5x5 (bottom right) square search window.MAD is conceived to be more robust to non-stationarity and abrupt morphological transitions, permitting a sharper representation of roughness.