POSITIVE AND NEGATIVE ROUGHNESS ACCORDING TO LOCAL DIFFERENCES BETWEEN DEM SURFACE AND 3D REFERENCE PLANES

The irregularities of the earth’s surface are quantified by means of roughness measurements using Digital Elevation Models (DEM’s). This article presents a roughness measurement method that is based on the calculation of the difference of altitude existing between a plane passing through the centre of a moving window and the altitude of the DEM surface inside this window. This method differs from the measure of the standard deviation and best fit plane, in the sense that it considers all difference values, positives or negatives. The measurement is done in a 3 × 3 or a 5 × 5 moving window and contemplates inside this window the plane which passes through the centre of the window and the highest pixel located in the border or perimeter of this window. According to the 3D configuration of the DEM surface inside the moving window, the sum of all the differences is positive or negative, allowing to discriminate the local morphology independently of the global roughness. The roughness variable which distinguishes negative and positive values allows to classify accurately landscape units such as watersheds, riverbeds, volcanic assemblages as well as landforms associated with tectonic structures.


INTRODUCTION
The notion of roughness is commonly used in quantitative geomorphology studies. This variable characterizes rock types, performs relief classifications, identifies landslide mechanisms, among others (Barnett et al., 2004, Berti et al., 2013. The roughness extracted from a Digital Elevation Model (DEM) is defined in terms of the variability of the elevation data. The main extraction methods are the root mean square applied to elevation and slope grids (Grohman, 2015), eigenvalue ratios (Cloude, 1999;Cloude et al., 2000), structuring elements (Cao et al., 2015), fractal dimension (Taud and Parrot, 2005), discrete Fourier transform (Bingham and Siegert, 2009), continuous wavelet transform and wavelet lifting scheme (Hani et al., 2011(Hani et al., , 2012, as well as standard deviation in a fit plane (Hobson, 1972;Evans, 1972;Herzfeld et al., 2000;Haneberg et al., 2005,). The results of these measurements depend on the resolution of the DEM and the size of the moving window (Grohmann et al., 2011).
Surface roughness describes elevation variations over a particular scale (Grohmann andHargitai, 2014, Grohman et al., 2015). To define the relationships existing between distinct curvature and various roughness characteristics, Hani et al. (2011Hani et al. ( , 2012 proposed to compute individually all the cells describing a Digital Elevation Model (DEM) and then realize a multiscale analysis. Their research showed that large scale surface roughness can be characterized using the distribution of slope frequencies, but a size increasing of the moving window generally used to measure roughness parameter is not particularly beneficial and in fact small pixels and small moving windows better reveals the surface condition. Regarding altitude acquisition techniques, the parameterization of radar data often * Corresponding author uses the root mean square (RMS) of slopes (Aharonson et al. 1998), but when it comes to laser altimetry sets, the absolute value of median slope is sufficient and provides a better characterization of surface slopes (Kreslavsky and Head, 1999). In fact, the roughness is generally considered as a geomorphological landform component, but landforms are often represented in heterogeneous thematic maps (De Reu et al., 2013;Melelli et al., 2017). For this reason, a numerical assessment provides a better quantitative characterization of a discretized space represented by continuous and regular cells, as it is the case for rasterized DEM.
However, in this research, it is proposed to generate a threedimensional plane from the real altitude values in each window size considering not only the altitude deviation, but also the direction of the reference planes. The results of the algorithm are compared considering 1) two sizes of moving window (3 × 3 and 5 × 5), 2) two different resolutions of DEM and 3) two roughness measurement methods providing results comparable to the best fit plane method: the three-dimensional fractal measurement (Taud and Parrot, 2005) and the calculation of the threedimensional surface (Parrot, 2007).

Study area
The study area is in the southern limit of the Trans-Mexican Volcanic Belt (TMVB), to the east of the state of Michoacán, in the central region of Mexico. This region corresponds to a transition zone between Oligocene (Sierra Madre Occidental) and Miocene to Quaternary volcanism (TMVB) (Blatter and Carmichael, 1998) that overlies Jurassic schists and Cretaceous limestones, marls, and shales of the Balsas River depression (Petersen, 1989;McLeod, 1989) (Fig. 1).
In the central area, the Tuxpan river runs in a tectonic control canyon in a north-south direction, and when it finds the transition from volcanic to sedimentary rocks, it quickly deepens its channel (Fig. 2). The volcanic rocks to the northeast correspond to quaternary basaltic volcanism forming shield volcanoes. The northwest region is a mountainous relief mainly formed by dacite from the Neogene (Siebe et al., 2007).

Figure 1. Situation map and block of the study area.
Underlying this unit, to the southwest, Paleogene rhyolitic rocks with oligomictic conglomerates form a mountainous range in the direction of the Balsas depression. Finally, to the southeast, are the limestone rocks and oligomictic conglomerates of the Upper Cretaceous, as well as sandstones-shales also from the Cretaceous (SGM, 1998).

Data used
Data used correspond to Digital Elevation Models ( Fig. 3) with resolution of 5 and 20 m generated by means of the interpolation of vector data (INEGI, 2015) using contour lines dilation method (Taud et al., 1999).

Overview
The method presented here is based on the difference of altitude existing between a plane passing through the centre of a moving window and the altitude of the DEM surface inside this window. This method differs from the measure of the standard deviation and best fit plane, in the sense that it considers all difference values, positives or negatives. The measurement of these differences is realized in a 3 × 3 or a 5 × 5 moving window and contemplates inside this window the plane which passes through the centre of the window and the highest pixel located in the border or perimeter of the window (Fig. 4). The calculation considers the difference between the altitude values of each pixel of this plane and the altitude values of each corresponding pixel of the DEM in the moving window and the sum of these differences allows characterizing the level of roughness of the local studied surface.
Depending on the result, it is possible to differentiate three classes: the first groups the positive values of the sum, the second ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-3-2021 XXIV ISPRS Congress (2021 edition) those which are negative, and the third group located in the centre of the histogram correspond to flat areas whatever the value of the slope.

Treatments
Inside the moving window, the algorithm first searches, comparing the altitude of the window centre, the highest elevation located on the edge of the window and calculates the difference between the two altitudes. The second step locates the position of this maximum to define the code used to generate the reference plane in the three-dimensional space.
As presented in the figure 5, 8 codes are attributed according to the position of the highest border point: code 1 for the NW, code 2 for the North, code 4 for the NE, code 8 for the East, code 16 for the SE, code 32 for the South, code 64 for the SW and code 128 for the West. This codification used here by analogy, is generally employed to define the configuration of the pixel assemblage surrounding the central pixel of a moving window (Tarboton, 1997;Jenson and Domingue, 1988). The same codes and the same cardinal positions are used in the case of a 5 × 5 moving window, that means that 8 intermediary positions in the border are neglected, but this restriction does not affect the validity of the result.
A quite simple generation of the plan of reference is as follows: the difference of altitude D between the altitude AC of the central point and the highest altitude value AH is considered as well as the code. In the case of a 3 × 3 moving window, three difference values are added or subtracted to the central pixel value to generate the plane of reference. In the case of a 5 × 5 window, 5 difference values are required to produce the plane of reference. Figure 6 provides all the patterns used to define the plane in a 3 × 3 window. In this figure, D corresponds to the value of difference added to AC, D/2 is the half of the difference added to AC, d and d/2 are negative values that correspond to the quantity subtracted to AC, 0 means unchanged AC value.   The sum of the differences within the moving window corresponds to a measure of the surface roughness that can be expressed in the following way: where R(i,j) is the local roughness of the pixel which coordinates are i and j (i for the lines and j for the columns); As is the altitude of the DEM and Ap the altitude of the reference plan in the same point, and m represents the window range. The range m allows defining the size ws of the window side according to the formula: 2 1.
When sometimes, the central pixel corresponds to the greatest hypsometric value, the reference plane is horizontal, and its altitude corresponds to the altitude of the central pixel. In this case the sum of all the negative differences will be evidently negative.
The frequency of the resulting roughness values is reported in a pyramid curve (Fig. 8). The roughness values between a lower value Rlow (in the example of treatment (resolution 5; moving window 3×3) shown in figure 10, Rlow is equal to -32.52) and 0 are normalized between 0 and 128 and values comprised between 0 and an upper value Rup (25.56 in the same example) are normalized between 128 and a maximum normalized value Nmax according to the following formula when Rup is lower than Rlow: in the present example, Nmax is 228.
An inverse calculation is done when Rup is greater than Rlow.
In the case of the figure 9 that corresponds to a treatment applied to a DEM with a resolution of 20 meters and a moving window 5 × 5, Rlow is equal to -371.98 and Rup is 343.23.
Afterwards, a lateral and symmetrical stretching is applied to visualized better the results ( Fig. 9 and Fig. 10).
The minimum value of the curve corresponds to crest lines and the maximum values to riverbeds. The highest point of the curve represents zones corresponding to a plane whatever the value of the slope (Fig. 8).

RESULTS AND DISCUSSION
We present two results here: the first one ( Fig. 9) corresponds to a treatment applied to a DEM with a horizontal resolution of 20 meters and a calculation using a 5 × 5 moving window; the second ( Fig. 10) is the result of a processing applied to a DEM with a resolution of 5 meters using a 3 × 3 moving window.
The process presented in this paper corresponds mainly to a textural treatment that emphasizes the structural features. For instance, the detail reported in the figure 11, clearly shows the position of the riverbed and the thalweg between the two lateral escarpments corresponding to the border of the fluvial terrace.  ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-3-2021 XXIV ISPRS Congress (2021 edition) Figure 11. Edges of the riverbed and terrace borders.
The distinction between various structural features results from the disjunction of the differences between the DEM surface and the reference plane as it has been previously defined. Thus, a textural factor can behave like a structural tool.
It should be noted that if we add the absolute values as shown in the following formula, we obtain a result comparable to that proposed for example by Hani et al. (2012) who estimates that the roughness can be considered as an expression of the vertical deviations by using the original surface and the surface of best fit whether based on all the curved subregions inside the moving window of observation or limited to only flat subregions.
In this case, as shown in figure 12, the representation of roughness is expressed globally and does not reveal the structural elements which determine it.

CONCLUSION
In fact, whether they are results coming from the global calculation or from calculations which differentiate negative values and positive values, areas with high roughness are easily defined. It must be noted, as shown in figure 13, that the most acute crest lines (values between 0 and 110 on the histogram of figure 8B) and the deepest valleys (values between 145 and 255 on the same histogram) are found in the same areas with strong roughness. It is also possible to extract the areas of low roughness by thresholding the values which correspond to the peak of the histogram (123 to 133). In figure 14, these zones appear in grey on which the weak slope zones (0 to 4°) have been plotted in black. In this image, a basaltic volcanic cone is well individualized (low roughness, slopes greater than 4°).

Figure 14.
Light grey: Low roughness regions with a slope greater than 4°: dark grey: Low roughness regions with a slope lower than 4°.
Moreover, the results obtained by using the notion of deviation from the plane passing through the centre of the moving window and the highest point of the perimeter of this window can be compared to those obtained using the three-dimensional fractal measurement (Taud and Parrot, 2005) (Fig. 15) and the calculation of the three-dimensional surface (Parrot, 2007;Parrot and Ramírez-Núñez, 2020) (Fig. 16).  The main contribution of the method presented in this article concerns the characterization of roughness at the local level, differentiating negative irregularities from positive irregularities, which allows to see the real weight of the structures in the relief. In terms of applications, this roughness variable allows detailed classification of landscapes (crest, riverbeds, volcanic domes) as well as its associated processes.