Scale-Optimized Surface Roughness for Topographic Analysis

Abstract

Surface roughness is a terrain parameter that has been widely applied to the study of geomorphological processes. One of the main challenges in studying roughness is its highly scale-dependent nature. Determining appropriate mapping scales in topographically heterogenous landscapes can be difficult. A method is presented for estimating multiscale surface roughness based on the standard deviation of surface normals. This method utilizes scale partitioning and integral image processing to isolate scales of surface complexity. The computational efficiency of the method enables high scale sampling density and identification of maximum roughness for each grid cell in a digital elevation model (DEM). The approach was applied to a 0.5 m resolution LiDAR DEM of a 210 km² area near Brantford, Canada. The case study demonstrated substantial heterogeneity in roughness properties. At shorter scales, tillage patterns and other micro-topography associated with ground beneath forest cover dominated roughness scale signatures. Extensive agricultural land-use resulted in 35.6% of the site exhibiting maximum roughness at micro-topographic scales. At larger spatial scales, rolling morainal topography and fluvial landforms, including incised channels and meander cut banks, were associated with maximum surface roughness. This method allowed for roughness mapping at spatial scales that are locally adapted to the topographic context of each individual grid cell within a DEM. Furthermore, the analysis revealed significant differences in roughness characteristics among soil texture categories, demonstrating the practical utility of locally adaptive, scale-optimized roughness.

1. Introduction

Surface roughness is an inherent property of topography and is commonly measured using land-surface parameters extracted from digital elevation models (DEMs) [1,2]. A range of DEM-derived surface roughness indices have been widely applied in geoscience and environmental research [3,4]. For example, topographic roughness maps have been used to delineate large-scale geological units and their age [4,5]. Roughness maps have been used to delineate landslides [6,7,8]. Surface roughness has also been widely applied to the study of surface processes in planetary science [9,10,11].

Two related topographic properties are often conflated in common usage of the term roughness [12]. First, roughness can denote local elevation variability, a concept that is referred to in this paper as ruggedness. Landscapes can be characterized along a ruggedness gradient from flat to variable-relief. Elevation range (i.e., local relief), standard deviation in elevation, and standard deviation of topographic residuals are common metrics used to characterize landscape ruggedness [9]. The second dimension of roughness is surface complexity, a measure of topographic texture. Topography can vary from smooth to irregular texture. Roughness metrics that characterize surface complexity either use surface area (e.g., topographic roughness index [1]), or variability in the surface normal vectors (or components of normals, e.g., slope and aspect). Vector dispersion and standard deviation of slope have been used previously to characterize surface complexity [4]. Ruggedness and complexity are orthogonal concepts because areas of complex texture can exhibit relatively low relief (e.g., the rough micro-topography of hummock and hollow peatlands) and vice versa.

Roughness maps are derived by measuring ruggedness or complexity related land-surface parameters within the local neighborhood surrounding each grid cell in a DEM. Therefore, roughness is commonly mapped using the same roving-window approach used for measuring many of the common topographic attributes [4,9]. The size of the local neighborhood dictates the scale at which surface roughness is characterized. The scale-variant nature of roughness is widely recognized in the literature [4,5,10]. Importantly, the scale-dependency of both ruggedness and surface complexity result from the fact that they are defined over areas rather than singular points. While some researchers have perceived the scaling of roughness as problematic [13], others recognize the opportunity scale-variant roughness provides in studying landscapes [4,5]. Ideally, roughness is assessed at a scale that is meaningful with respect to the scale of landforms, geomorphological processes, and the specific application [14]. Ultimately, for heterogeneous landscapes, a single optimal scale to measure surface roughness may not exist. Rather, a range of scales may be more appropriate for capturing the varying complexity of the topographic surface in a region.

The scale-variance of ruggedness is in part a result of the increased likelihood of including more prominent elevation minima and maxima with more extensive neighborhoods. As a result, ruggedness maps tend to highlight elevation minima/maxima and breaks-in-slope. Ruggedness scale signatures, which characterize the relation between roughness and filter kernel size, tend to be simple, monotonically increasing functions. By comparison, the scale signatures for texture-based roughness metrics are generally more complex and are often indicative of variation in the topographic expression of geological units and their age, landform arrangements, and land use/land cover that exhibit micro-topography [15,16].

Grohmann et al. [4] used a combination of filtering and data resampling (i.e., DEM grid coarsening) to study the multiscale nature of roughness for a site in Scotland. They concluded that researchers should derive roughness parameters from higher resolution data, and resample to a coarser resolution as needed. However, while data coarsening significantly improves the computational efficiency of large-scale operations, this approach comes at the cost of potentially losing topographic information. Wood [17] performed multiscale measurements of indices of topographic shape by using first and second derivatives of quadratic surfaces fitted over a range of scales. Lindsay et al. [18] and Newman et al. [19] showed how efficient filtering methods based on integral images [20] can be applied to measures of local topographic position, allowing for efficient ‘hyper-scale’ analysis (analogous to hyper-spectral analysis in the field of remote sensing) of topographic properties without the need for data coarsening. This paper extends these earlier studies by examining the hyper-scale properties of topographic complexity and by developing a method to derive locally adaptive, scale-optimized roughness maps. An approach is presented for measuring a normal vector-based roughness index. This roughness index is applied to the topographic analysis of a fine-resolution LiDAR data set case study.

2. Materials and Methods

2.1. Locally Adaptive Scale-Optimized Surface Roughness Measurement

A surface normal is a 3-D vector that is perpendicular to the surface tangent plane at a point. The three components (x, y, and z) of a unit surface normal are often derived from the slope (S) and aspect (a) of the surface, such that [4,21]:

$$x=\sin \left(a\right)\times \sin \left(S\right)y=\cos \left(a\right)\times \sin \left(S\right)z=\cos \left(S\right),$$

(1)

In practice, however, the three vector components can be derived directly from the elevations within the 3 × 3 window surrounding each grid cell in a DEM, thereby avoiding the calculation of the intermediate S and a terms and negating the need for the computationally expensive trigonometric functions.

Surface complexity can be related to the dispersion of the surface normals within an area; where planar surfaces exhibit no variation in normal orientations and increased irregularity produces greater dispersion (Figure 1).

The spherical standard deviation (σ_s) is a measure of the angular spread among n unit vectors and is defined as [22]:

$${\sigma }_s=\sqrt{-2\ln \left(\frac{R}{n}\right)}·\frac{180}{\pi },$$

(2)

Here, σ_s is measured in degrees and is zero for flat and inclined planes and increases infinitely with greater surface complexity. This formulation of σ_s assumes an underlying wrapped normal distribution. R is the resultant vector length and is derived from the sum of the x, y, and z components of each of the n normals contained within a filter kernel:

$$R=\sqrt{{\left({{\displaystyle \sum }}_i^{n}x\right)}^2+{\left({{\displaystyle \sum }}_i^{n}y\right)}^2+{\left({{\displaystyle \sum }}_i^{n}z\right)}^2},$$

(3)

Since normals are unit vectors, R varies from 0 to n, where larger values indicate less dispersion.

The filter size (D) used to calculate the component summations in Equation (3) determines the scale of roughness characterization and therefore n = D². D must be an odd integer equal to or greater than three. For multiscale analyses, it can be more convenient to refer to kernel sizes using the filter radius (r), where D = 2r + 1 [18]. Thus, r = 1 coincides with the smallest possible 3 × 3 filter, r = 2 is a 5 × 5 filter, etc. Equation (2) can therefore be thought of as defining σ_s(r), a function referred to in this paper as the roughness scale signature.

Grohmann et al. [4] found that vector dispersion, a similar measure of angular spread of surface normals, increases nearly monotonically with increasing scale. This occurred because their measure accumulates the variance of all smaller scales up to the test scale. A more interesting scale relation can therefore be estimated by isolating the amount of surface complexity associated with specific scale ranges. That is, at large spatial scales, σ_s should reflect the texture of large-scale landforms rather than the accumulated complexity at all smaller scales, including high-frequency micro-topographic roughness [14]. As such, the approach adopted in this paper is to suppress the surface complexity of scales that are smaller than the filter size by applying Gaussian blur (with a standard deviation of r/3) to the DEM prior to calculating σ_s(r). In this way, our method is able to isolate and highlight the surface complexity associated with landscape features of a similar scale to that of the filter size. Wu et al. [9] applied a similar approach to scale partitioning in the study of surface roughness, except a mean filter was used to suppress higher-frequency signals in place of a Gaussian filter. Gaussian filters are commonly preferred in computer vision and image processing applications for scale space representation because they do not introduce new spurious structures at coarse scales that do not correspond to simplifications of structures at finer scales [23].

The maximum σ_s (σ_{s max}) for a tested range of spatial scales from lower (r_L) and upper (r_U) bounds can be defined as:

$${\sigma }_{smax}=\max \left\{{\sigma }_s\left(r\right):r=r_L\dots r_U\right\},$$

(4)

Equation (4) represents a locally adaptive measure of surface roughness, in which each grid cell in a DEM is provided a scale-optimized roughness value that is suited to the topographic context in which the cell is situated. If a grid cell in a DEM is situated among the highly complex micro-topography of a tilled agricultural field and another is situated within the much larger-scale undulating surface of a moraine, σ_{s max} can potentially recognize the peak roughness in their local scale signatures to identify their dominant scales.

Calculation of σ_s for larger filter kernels can be very computationally intensive if a naïve filtering approach is used. Fortunately, the component summations in Equation (3) can be computed using integral images, or summed-area tables [20]. An integral image is a raster data structure in which each grid cell’s value represents the sum of an underlying distribution (a normal component in this case) within the rectangular area bounded by the cell and one of the image corners (Figure 2a). Widely applied in the field of computer vision, integral images allow for efficient measurements of the total of the underlying distribution for arbitrary sized rectangular areas in three mathematical operations (Figure 2b). This can facilitate image filtering with constant-time computational efficiency, varying with the number of grid cells in the image, but independent of the size the filter kernel. This is the basis for efficient, hyper-scale geomorphometric analysis proposed by [18].

A plugin, named MultiscaleStdDevNormals (referred to here as MSSDN), was implemented within the open-source geospatial analysis library WhiteboxTools [24] using the Rust programming language. This tool can be used to map the spatial pattern of maximum spherical standard deviation, as well as the scale at which maximum spherical standard deviation occurs (r_max), for each grid cell in a DEM. The general procedure operates as described in

Table 1. Procedure used by the MSSDN tool to calculate σ_{s max} and r_max.

Algorithm Description

Initialize output magnitude (σs max) and scale rasters (rmax), setting cell values to −1.0. For the tested filter size (rL), smooth the input DEM with Gaussian blur. Calculate surface normals for each grid cell in the smoothed DEM and store these data in separate x, y, and z component images. Convert each component image to an integral image. Calculate σs using the tested filter size. For each grid cell in the σs raster, if the value is greater than the current values in the output magnitude raster, update the output and set the scale raster to the tested filter size. Repeat steps 2 to 6 for each filter size in the user-specified scale range up to rU.

.

The MSSDN tool makes extensive use of integral images (i.e., summed-area tables) and parallel processing to ensure computational efficiency. Normal vector calculation (Step 3) is a 3 × 3 neighborhood operation, and is therefore relatively computationally efficient and can be parallelized. Each integral image (Step 4) is calculated sequentially, however, integral-image based filtering operations can be parallelized for further improved efficiency. The DEM smoothing operation (Step 2) utilizes the fast almost-Gaussian filter [25], which applies multiple integral-image-based mean filters of differing sizes to closely approximate Gaussian blur.

2.2. Study Site and Data

Locally adaptive scale-optimized roughness patterns were evaluated for an 840,000,000-grid cell (3.4 GB), 0.5 m resolution LiDAR DEM. The DEM characterizes topography in a 210 km² area east of Brantford Ontario, Canada (Figure 3). The study site is situated at latitude and longitude ranges of 43°5′48” N to 43°13′59” N and 80°2′22” W to 80°12′49” W respectively and overlaps with extensive portions of the Fairchild Creek and Big Creek sub-basins, which are both minor tributaries of the Grand River. The northeastern portion of the study site exhibits a fine drainage texture in the headwaters of Big Creek, where streams are deeply incised. The southern area of the site is dominated by the large incised meanders and broad floodplains of Fairchild Creek and the Grand River. The steep cut banks associated with these features provide significant local relief. The overall relief of the study site is approximately 80 m. The physiography of the site is divided between dissected sand plains and clay plains, with a small area of kame moraine in the northeastern corner of the site [26]. The dominant soil types are silt loam and sandy loam with smaller, isolated areas of silty clay loam and silty clay. The majority of the study site has agricultural land uses; forested areas comprise approximately 14% of the site and is often coincident with the relatively steep slopes of incised streams. There are some smaller urbanized areas within the site, particularly associated with the outlying portions of Brantford.

The DEM was interpolated using a triangulation gridding scheme from last- and only-returns of the source LiDAR point cloud, excluding classified non-ground points (largely vegetation and buildings). While last-returns are generally associated with the ground surface, under heavy forest cover, these points may also include the low-lying portion of tree trunks. This may partially contribute the high roughness of the interpolated DEM in forested areas. The publicly available source data were collected during leaf-off and snow-free ground conditions in the spring of 2018 by a private contractor commissioned by the Ontario Ministry of Agriculture, Food and Rural Affairs (OMAFRA) and the Ministry of Natural Resources and Forestry (MNRF). The LiDAR data were stored in the NAD83 UTM zone 17N (EPSG:2958) projected coordinate system. The average point density of the data set is 8 points⋅m⁻² and the overall vertical accuracy was estimated to be 5 cm. A triangulation gridding scheme was used to interpolate the point data. Hydro-flattening was used to correct the surface of the smaller waterbodies and rivers within the area. WhiteboxTools’ FlattenLakes tool and the publicly available Ontario Hydro Network’s waterbodies data set were used to perform this operation. While many buildings were excluded from the interpolation by removing non-ground point classifications, the DEM was also treated to remove any remaining off-terrain objects using the RemoveOffTerrainObjects tool in WhiteboxTools. Several major roads transect the site, and their embankments are apparent in the DEM. The resolution of the DEM is sufficient to represent the micro-topographic roughness associated with tillage patterns in fields and undulating forest floors. At larger spatial scales, the incised channels provide interesting landscape-scale texture for this study.

3. Results

The MSSDN tool was used to map the spatial distributions of σ_{s max} (Figure 4) and r_max (Figure 5) for the test DEM. Processing completed in 71 min 28.1 s (excluding file input/output time) using an 8-core 3.0 GHz processor with 64 GB of memory. Each tested scale took approximately 85 s to evaluate. The maximum memory requirement during processing was 10.5 times the DEM size (note, the DEM file was stored as 32-bit floating-point values). Fifty scales were sampled, with 1 ≤ r ≤ 748 (in grid cells). Because the DEM grid resolution is 0.5 m, the scale filter radius, in cells, is approximately equivalent to the full filter size (D) in meters. Experience with roughness scale signatures has shown that σ_{s max} is highly variable at shorter scales and changes more gradually at broader scales. Therefore, a nonlinear scale interval was used to ensure that the scale sampling density was higher for short scale ranges and coarser at longer tested scales, such that:

$$r_i=r_L+{\left(i-r_L\right)}^p,$$

(5)

where r_i is the filter radius for step i and p is the nonlinear scaling factor, set to 1.7 in this testing. A constant sampling interval (e.g., r = 1, 3, 5, 7, etc.) was also explored; however, there tends to be more variance at smaller scales, due to micro-topography, than at longer spatial scales. Therefore, it is difficult to select a single sampling interval across the entire scale range without either under-sampling at the shorter scale ranges or needlessly over-sampling at longer ranges. There is nothing significant about the nonlinear scaling factor of 1.7 used in this study except that, for the test DEM, this value provided a good compromise of high sampling density at shorter scales and coarser sampling at the longer scales within the tested scale range. If roughness scale signatures appear smoothly varying, the selected nonlinear scaling value is adequate; if the signature appears more jagged, it is likely that the scaling value should be decreased due to under-sampling.

The highest σ_{s max} values were found to occur along the steep cut banks and terrace scarps of the deeply incised rivers within the site (Figure 4). By comparison, the wide floodplains and scroll-bars of Fairchild Creek and the Grand River are among the least texturally complex locations within the site, being broadly level, depositional topography. Each of the ravines associated with the smaller incised channels within the northeast of the site demonstrated relatively high roughness, typically peaking at intermediate to broader spatial scales (Figure 5). Road embankments, and particularly their associated roadside ditches, frequently exhibited high σ_{s max} values (Figure 4b). The irregular ground associated with forest floors demonstrated elevated roughness relative to their surrounding agricultural areas; an example of this phenomenon is illustrated in the detailed inset provided in Figure 4b. The fine texture associated with tillage patterns within agricultural fields were also evident in the σ_{s max} map, and in fact, the differing tillage direction and furrow depths of neighboring fields, are discernable in the image (Figure 4d).

The spatial distribution of r_max (Figure 5) illustrated several characteristics about σ_{s max}. First, flatter agricultural fields containing prominent tillage patterns and forested areas reached peaks in σ_s at shorter scales, typically with 1 ≤ r ≤ 4. These short scales are very extensive in the study site owing to the prevalence of these land covers (Figure 6). The longest spatial scales of peak σ_s largely coincided with the smooth floodplains on the inside of meander bends (point bars). Fields exhibiting rolling topography possessed more intermediate (20 ≤ r ≤ 70) scales of σ_{s max}. These scales also coincided with a secondary peak in the probability density function (Figure 6) located at r=51 (25.5 m), indicating the prevalence of these roughness scales within the site.

Approximately 35.6% of the grid cells were found to have r_max values less than four grid cells, indicating the dominance of micro-topographic roughness within the 0.5 m resolution study area DEM. The peak in σ_s values was very rarely identified at filter radii greater than approximately 300 cells (150 m) within the study site and only a small proportion of the grid cells were assigned σ_{s max} at the tested r_U value (Figure 6). This would indicate that the r_U value used in this case study (r_U = 748) was sufficient to capture the vast majority of the dominant scales of roughness within the evaluated terrain. The condition of r_max = r_U may be indicative of a truncated peak in the scale signature at the broadest tested scales. Thus, it is advisable to set the r_U value high enough that only a small proportion of grid cells in the DEM are assigned this upper boundary value. Only 0.36% of the grid cells in the study site DEM were found to have r_max = r_U (Figure 6).

Detailed scale signatures were extracted for five sites within the larger study area, with varying topographic and land-cover characteristics (Figure 7). Each of the signatures possessed at least one well-defined peak in roughness and demonstrated decreased σ_s variability for scales greater than approximately 500 m. The diversity in overall signature shape and peak location was found to be high. Site 1 was situated on the inside bend of a meander (Figure 7a) and exhibited generally low σ_s values for r ≤ 350, peaking at r = 508 (Figure 7b). Site 2 was located atop a low hill in an agricultural field typical of the rolling topography found in places throughout the study site. The tillage pattern in the Site 2 field was not as prevalent as in other sites. The site’s signature (Figure 7b) possessed two peaks; the first and highest peak was located at r = 28, while the second much lower peak occurred at r = 344. This first peak was associated with the dimension of the low hill on which the site was situated, while the lower secondary peak was related to a larger-scale topographic undulation within the field. Site 3 was in the middle of a woodlot with a highly irregular forest floor. The first peak in roughness at Site 3 occurred at r = 2 and was the second highest measured roughness peak of the tested signatures. This peak was associated with the micro-topography of the forest floor, while the secondary peak, situated at r = 150, resulted from the larger-scale sloping topography of the woodlot. Site 4 was located along a small incised headwater stream, which passed through an agricultural field. The roughness related to the prominent tillage pattern within the field is apparent in the spike in σ_s value at the smallest scales, but the dominant peak in the scale signature occurs at r = 90 and was caused by the relatively steep hillslopes on either side of the channel. Finally, Site 5 was situated within a flat agricultural field. Unsurprisingly, in the absence of larger-scale topographic features, σ_{s max} occurs at r = 1 and reflects the tillage within the field. Roughness was consistently low across all scales for Site 5, although it did exhibit a slight rise in σ_s for the range 50 ≤ r ≤ 239, before gradually dropping again at larger scales (Figure 7b).

Analysis of the patterns of σ_{s max} and r_max revealed negative exponential relation (Figure 8). This clearly demonstrates the inverse relation between the two indices, such that the grid cells exhibiting the largest σ_{s max} values (> 20°) are typified by r_max < ~50, and similarly, that sites with r_max > 250 exhibit σ_{s max} < ~ 5°.

A comparison of the summary statistics of σ_{s max} (

Table 2. Relation between various soil texture categories within the study site and σ_{s max}.

Soil Texture	Mean	Median	Min	Max	Range	SD
silt loam	3.88	3.05	0.28	92.87	92.59	3.25
sandy loam	3.03	2.36	0.23	82.32	82.08	2.30
silty clay loam	2.59	1.92	0.24	75.45	75.21	2.23
organic	4.83	4.20	1.20	42.34	41.14	2.48
silty clay	2.14	1.59	0.25	72.00	71.75	1.75
loam	3.95	2.72	0.27	80.52	80.25	4.43
fine sandy loam	4.44	3.98	2.50	29.68	27.18	2.01
very fine sandy loam	3.45	2.47	0.64	30.46	29.82	2.56
loamy sand	2.25	1.50	0.60	25.51	24.91	1.95

) and r_max (

Table 3. Relation between various soil texture categories within the study site and r_max.

Soil Texture	Mean	Median	Min	Max	Range	SD
silt loam	62.0	35	1	748	747	92.8
sandy loam	62.4	28	1	748	747	114.2
silty clay loam	59.8	16	1	748	747	101.2
organic	46.1	35	1	344	343	43.8
silty clay	61.1	4	1	748	747	106.0
loam	262.0	325	1	748	747	227.2
fine sandy loam	152.3	178	1	363	362	107.7
very fine sandy loam	43.5	4	1	508	507	63.8
loamy sand	149.4	112	1	599	598	150.6

) distributions was made among the various soil texture categories within the study site. One-way ANOVA testing showed that the soil texture category means were statistically heterogeneous for both σ_{s max} (F_α = 0.05, df₁ = 8, df₂ = 776,113,757 = 4,597,316.144, p < 0.0001) and r_max (F_α = 0.05, df₁ = 8, df₂ = 776,113,757 = 1,560,163.151, p < 0.0001). Of course, given a sufficiently large sample, extremely small and non-notable differences can be found to be statistically significant and statistical significance says nothing about the practical significance of a difference. Nonetheless, examination of Table 2 and Table 3 does reveal noteworthy differences among the distributions of the two roughness indices. For example, silty clay, silty clay loam, and loamy sand each exhibited relatively low average σ_{s max} while organic soils were notable for their relatively high average σ_{s max}. Each of the sandy textures appeared to have substantially less range in σ_{s max} than the other categories. Organic soils exhibited less variability in r_max values while the median r_max for silty clay and very fine sandy loam soils were notably lower than other categories. Loam, fine sandy loam, and loamy sand were discernable for their longer mean and median r_max values. As a whole, the r_max distributions within each soil category were more skewed (see mean relative to median) than the corresponding σ_{s max} values.

4. Discussion

The extensive use of integral images and parallelization within the MSSDN tool makes applying this approach to study the hyperscale characteristics of surface roughness in heterogeneous landscapes practical without the need to adjust DEM resolution using resampling or other data coarsening techniques [4]. The case study shows that it is possible to sample scales densely, even for massive DEMs, allowing for locally optimal scale selection. However, while the approach used in this study was found to be computationally performant, it does require significant memory resources. The large memory requirement of the tool results from the need to store each of the integral images of the three normal vector components with double-precision (64-bits) to avoid rounding errors. Therefore, multi-gigabyte DEMs may not be able to be processed on memory-constrained systems without first splitting the data into smaller pieces. Memory usage issues associated with the application of integral images have been widely recognized in the literature [27,28]. This may become less of a restrictive issue as computing platforms with greater memory become more commonly available and as researchers continue to move toward cloud computing platforms, such as Google Earth Engine, for spatial analysis [29].

Areas of relatively high roughness largely coincided with the smaller incised channels and the cut banks and terrace scarps of larger rivers, as well as road embankments. The impact of road embankments emphasizes the importance of removing these features from fine-resolution LiDAR DEMs in certain applications [30]. The areas of the study site that were under forest cover were each associated with moderate to high values of roughness. This is indicative of how well the source LiDAR data were able to capture the rough ground beneath the forest canopy, a characteristic of LiDAR data sets that has been previously observed in the literature [7,31]. Generally, areas of high roughness were associated with shorter spatial scales and vice versa (Figure 8), with a relation described by a negative exponential distribution. Some of the longest scales of σ_{s max} in the site were situated along the broad floodplains and meander scroll-bars, which were also among the least topographically complex (i.e., smoothest) sites. These highly level sites experience maximum roughness at larger spatial scales, such that search windows overlap with the relatively more complex surrounding topography.

The combination of the maximum roughness (Figure 4) and the scale of maximum roughness (Figure 5) aided in the interpretation of topographic properties of the study site. While roughness is commonly used as an input parameter for modelling applications, we suggest that the additional information of the scale of peak roughness (r_max) may provide additional useful information on landscape structure for these applications. The mosaic of spatial scales associated with maximum surface complexity (Figure 5) revealed substantial information about the topographic character of the test DEM (Figure 3) and was demonstrated to be linked with soil texture (Table 2 and Table 3). The widely ranging r_max values show that there is no single ‘representative’ kernel size that can be chosen to characterize topographic texture in a complex and extensive area. As fine-resolution LiDAR DEMs become increasingly common in terrain modelling applications [32,33], multiscale analysis of land-surface parameters has become more important [12,34,35] and has the potential for improving the accuracy of mapping and classification applications. Future research should explore the extent to which locally-adaptive scale selection, such as the scheme used in this paper, can also impact geomorphometric mapping and classification applications.

By evaluating this locally-adaptive scale selection method for characterizing roughness (i.e., site-specific scale of maximum roughness) it is evident that the information content in the MSSDN map (Figure 4) cannot be replicated by the common approaches of 1) measuring roughness at an, often arbitrary, single ‘representative’ scale, or 2) using a stack of roughness maps derived across a range of spatial scales. This last point may not be immediately evident, particularly given the fact that the MSSDN map is itself derived from the multiscale image stack. However, given the heterogeneity observed in the scale signatures among DEM grid cells (e.g., Figure 7), in diverse terrain, any classification model that is based on applying weights to each image in a roughness scale stack cannot represent the same information as contained in an MSSDN image. Furthermore, there is potential to extend this approach to derive additional information based on other characteristics of scale signatures, e.g., roughness variability, number of peaks, peak height ratios, etc.

The approach to multiscale roughness analysis used in this research is highly adaptive to different DEM data sets and resolutions, as well as applications. For example, while micro-topographic roughness was found to dominate much of the study area, if the larger-scale surface complexity associated with landform features were of specific interest, it would be possible to remove the impacts of smaller-scale roughness by selecting a larger r_L value. If, instead, an application warrants focusing on surface roughness at shorter scales (e.g., Campbell et al.’s [36] recent examination of the impact of roughness on human travel through forested terrain), then it would be appropriate to select a relatively low r_U value. It should be possible to similarly apply this technique to study regional-scale roughness patterns at the continental scale using available global DEM products. That is, the lower and upper bounds of the test scales should be suited to the specific application. The use of Gaussian blur to partition roughness scales ensures that roughness is not accumulated from smaller spatial scales to impact the interpretation of features at a larger scale.

5. Conclusions

To conclude, the scaling of surface roughness (i.e., surface complexity) is shown here to exhibit substantial heterogeneity within a 210 km² study area. This study demonstrated an application of roughness mapping across a broad range of spatial scales with a fine scale resolution. This method allows for the characterization of maximum surface roughness at spatial scales that are optimal for each individual grid cell within a DEM. These data were found to be reflective of the diverse topographic and land-use settings within the study site. Furthermore, the information contained within the scale map can provide additional useful information for landscape interpretation. The analysis revealed significant differences in the maximum roughness and scale of maximum roughness among the various soil texture categories within the study site, demonstrating the practical utility of locally adaptive, scale-optimized roughness characterization for mapping applications.