Measuring Hyperscale Topographic Anisotropy as a Continuous Landscape Property

Several landforms are known to exhibit topographic anisotropy, defined as a directional inequality in elevation. The quantitative analysis of topographic anisotropy has largely focused on measurements taken from specific landforms, ignoring the surrounding landscape. Recent research has made progress in measuring topographic anisotropy as a distributed field in natural landscapes. However, current methods are computationally inefficient, as they require specialized hardware and computing environments, or have a limited selection of scales that undermines the feasibility and quality of multiscale analyses by introducing bias. By necessity, current methods operate with a limited set of scales, rather than the full distribution of possible landscapes. Therefore, we present a method for measuring topographic anisotropy in the landscape that has the computational efficiency required for hyperscale analysis by using the integral image filtering approach to compute oriented local topographic position (LTP) measurements, coupled with a root-mean-square deviation (RMSD) model that compares directional samples to an omnidirectional sample. Two tools were developed: One to output a scale signature for a single cell, and the other to output a raster containing the maximum anisotropy value across a range of scales. The performances of both algorithms were tested using two data sets containing repetitive, similarly sized and oriented anisotropic landforms, including a dune field and a drumlin field. The results demonstrated that the method presented has the robustness and sensitivity to identify complex hyperscale anisotropy such as nested features (e.g., a drumlin located within a valley).


Introduction
The term anisotropy describes a directional inequality, occurring when a parameter varies with orientation.Subsequently, topographic anisotropy describes how the elevation field varies with orientation.Topographic anisotropy can be expressed across a wide range of spatial scales; from individual anisotropic landforms, such as drumlins, dunes, valleys and ridges, to entire landscapes in which the arrangement of landforms creates an oriented surface texture.Topographic anisotropy is an important property in several geoscience fields, especially as it relates to geomorphometry, geology [1,2], hydrology [3], and geomorphology [4][5][6][7].Therefore, the ability to measure and characterize topographic anisotropy over a landscape extends beyond specific landforms as a characteristic of the encompassing landscape.
Several methods have been developed to measure topographic anisotropy.Object-based methods are a common way to quantify topographic anisotropy for landforms due to the ease of implementation.For example, Maclachlan and Eyles [6] used topographic data to delineate drumlins, calculating anisotropy in the form of an elongation ratio from the maximal and orthogonal lengths of the drumlin objects.Alternatively, through the patterns observed in cell-based empirical directional variograms, geostatisticians recognized that spatial correlation is dependent on both distance and relative orientation [8,9].As a result, semivariogram attributes, such as the range and the sill, are frequently used to define spatially distributed anisotropy as range-and sill-types of anisotropy.Sill refers to the finite limiting semivariance as the lag distance increases, and range refers is the lag distance at which the sill is reached [10].The more common and intuitive range-type anisotropy measures anisotropy as the ratio of the maximum range with the range of the orthogonal angle, which is assumed to be the minimum range [9][10][11].Therefore, range-type anisotropy measures the maximum effect of orientation on semivariance.However, both types show the correlation length of the topographic surface differs by orientation.
Recent research has continued to develop the directional semi-variogram for hyperscale range-type topographic anisotropy measurements [12].In a hyperscale framework, Roy et al. [1] demonstrated that topographic anisotropy is a distributed and scale-dependent landscape property using the every-direction variogram analysis method.Houser et al. [7] expanded upon this method to sample the standard deviation of elevation across scalable linear transects as a proxy for relief.However, computing fully distributed anisotropy measurements across a range of both scales and rotation axes results in high computational complexity.To circumvent this complexity researchers are forced to either (1) reduce the number of scales or rotations analyzed [7,12], or (2) use specialized hardware or computing environments to reduce the execution time [1].Further discretizing the parameters degrades the quality of the analysis by creating trade-offs between parameter resolution and accuracy.Leveraging more processing power, while offering a powerful option to reduce execution time, is not as practical or scalable compared to algorithmic efficiency.Neither of these methods have the efficiency required to carry out hyperscaled analyses without specific hardware pre-requisites, and therefore introduce barriers to adoption by researchers and practitioners alike.
Where object-based methods omit the surrounding landscape [6,13], local topographic position (LTP) metrics intrinsically relate the elevation of a point to the topography of the surrounding area.Many disciplines such as pedology [14], geomorphometry [15], ecology [16,17], and precision agriculture [18,19], among others, use the relationship with the surrounding landscape to characterize terrain.LTP metrics, including the topographic position index (TPI) and deviation from mean elevation (DEV) provide this information by quantifying how much higher or lower each point is relative to the surrounding area [20].Using a hyperscale analytical framework to evaluate LTP across a range of window sizes allows for the characterization of a point's relief across a range of scales [15,21].Since LTP is fundamentally defined by the scale of analysis, applying hyperscale techniques to LTP analyses not only improves the analysis' accuracy by creating a data-driven scale selection process [22][23][24], but also provides scale information with which to characterize landscape anisotropy [20].
While hyperscale LTP reduces error due to scale selection, it is prone to high computational requirements similar to the semivariogram methods.However, advances in spatial filtering such as the Huang histogram approach [25] and the integral image [26] can be applied to greatly improve computational efficiency.The increased efficiency allows feasible hyperscale analyses to be conducted across a nearly continuous range of scales, eliminating user scale selection bias.Hyperscale sampling eliminates the need for scale selection, instead either analyzing all scales or self-selecting key scales for each grid cell.Newman et al. [27] determined that the DEV metric implemented with the integral image filtering technique was the optimal LTP metric for hyperscale analysis, due to highly robust characterization of LTP and scale-independent time complexity.Moreover, DEV is a relative metric that standardizes the elevation residuals by the elevation variability, and therefore is unaffected by the magnitude of the elevation values.This research presents a practical and generalizable method to characterize hyperscale landscape anisotropy, using a directional DEV sampling strategy and a conservative anisotropy model based on the common root-mean-square deviation (RMSD) that compares all sampled directions simultaneously.

Directional LTP Sampling
Directional LTP sampling was achieved by adapting the integral image method used by Lindsay et al. [20] to create a 3 × 3 matrix of scalable filters around a central cell (z 0 ) (Figure 1).The integral image filtering technique offers constant-time calculations for the number, mean, and variance of elevation samples within the filter, with the trade-off of using strictly rectangular window geometry [20,28,29].The algorithm determines filter geometry by iterating over the range of filter half-lengths specified by the user, where the total length of the filter (L), in cells, is given by: where L is the filter length, and R is the half-length of the filter, both in units of grid cells.The total filter length was partitioned into three sections for each axis (row and column), resulting in a 3 × 3 matrix of sub-filters.The filter lengths of the original filters are used only as an index describing scale, and are not used in the characterization of anisotropy.Directional sampling was achieved by additively combining sub-filters to create four aggregate windows, oriented along the cardinal directions (north-south and east-west oriented windows, as shown in Figure 1a,b respectively), both diagonal directions (northwest-southeast and northeast-southwest oriented windows, as shown in Figure 1c,d respectively), and an isotropic square sample containing all nine sub-windows (Figure 1e).It should be noted that since each oriented window is composed of three sub-filters, the area of all oriented windows is equal.The DEV LTP metric was then calculated for each oriented window with more than two valid elevation samples (variance requires three or more samples), though this limit can be increased to reduce instability in the variance terms at small sample sizes.DEV can be calculated as: DEV = (z 0 − µ)/σ, where z 0 is the central analytical cell, µ is the window mean value, and σ is the window standard deviation.

Directional LTP Sampling
Directional LTP sampling was achieved by adapting the integral image method used by Lindsay et al. [20] to create a 3 × 3 matrix of scalable filters around a central cell (z0) (Figure 1).The integral image filtering technique offers constant-time calculations for the number, mean, and variance of elevation samples within the filter, with the trade-off of using strictly rectangular window geometry [20,28,29].The algorithm determines filter geometry by iterating over the range of filter half-lengths specified by the user, where the total length of the filter (L), in cells, is given by: where L is the filter length, and R is the half-length of the filter, both in units of grid cells.The total filter length was partitioned into three sections for each axis (row and column), resulting in a 3 × 3 matrix of sub-filters.The filter lengths of the original filters are used only as an index describing scale, and are not used in the characterization of anisotropy.Directional sampling was achieved by additively combining sub-filters to create four aggregate windows, oriented along the cardinal directions (northsouth and eastwest oriented windows, as shown in Figures 1a and 1b respectively), both diagonal directions (northwestsoutheast and northeastsouthwest oriented windows, as shown in Figures 1c and 1d respectively), and an isotropic square sample containing all nine sub-windows (Figure 1e).It should be noted that since each oriented window is composed of three sub-filters, the area of all oriented windows is equal.The DEV LTP metric was then calculated for each oriented window with more than two valid elevation samples (variance requires three or more samples), though this limit can be increased to reduce instability in the variance terms at small sample sizes.DEV can be calculated as: where z0 is the central analytical cell, µ is the window mean value, and σ is the window standard deviation.Windows that are not centered over the central analytical cell (z 0 ) measure the difference between the z 0 and the adjacent areas, rather than relating it to the surrounding area.Therefore, windows must be centered to avoid measuring what is effectively an area-weighted gradient.This gradient exists along the spatial limits of the domain with an elevation gradient due to the disproportionate inclusion of 'no data' values in some orientations, skewing the directional LTP value by the elevation gradient.This produces a notable edge effect that results in artificially large anisotropy values.To mitigate this distortion, the edges were buffered by the filter half-length such that no samples were taken from beyond the extent of the data.However, buffering the edges forces z 0 towards the center of the image as the window size increases, limiting the maximum window size along the perimeter.
Topographic anisotropy can be thought of as variation in elevation as a function of orientation, defined between the interval [0, 2π] radians.Since centered windows extend equally in opposite directions about z 0 , all directions are sampled from rotation axes between [0, π] radians.A continuous sampling strategy would compute a value for all rotation axes.However, DEV is calculated over an area and must sample swaths of rotation axes.More specific to the sampling strategy described above, rotation axes zero to π radians are sampled using four windows-providing a DEV value for every π/4 radian (45 • ) swath (Figure 2).Conversely, isotropic sampling integrates all orientations, computing a local average for the entire neighborhood.This sampling strategy achieves the efficiency required for hyperscaled analysis through the integral image filtering technique, which comes at the cost of coarse orientation resolution.Therefore, measures of the orientation of anisotropy have been excluded from the method due to the coarse resolution orientation sampling.
Windows that are not centered over the central analytical cell (z0) measure the difference between the z0 and the adjacent areas, rather than relating it to the surrounding area.Therefore, windows must be centered to avoid measuring what is effectively an area-weighted gradient.This gradient exists along the spatial limits of the domain with an elevation gradient due to the disproportionate inclusion of 'no data' values in some orientations, skewing the directional LTP value by the elevation gradient.This produces a notable edge effect that results in artificially large anisotropy values.To mitigate this distortion, the edges were buffered by the filter half-length such that no samples were taken from beyond the extent of the data.However, buffering the edges forces z0 towards the center of the image as the window size increases, limiting the maximum window size along the perimeter.
Topographic anisotropy can be thought of as variation in elevation as a function of orientation, defined between the interval [0, 2π] radians.Since centered windows extend equally in opposite directions about z0, all directions are sampled from rotation axes between [0, π] radians.A continuous sampling strategy would compute a value for all rotation axes.However, DEV is calculated over an area and must sample swaths of rotation axes.More specific to the sampling strategy described above, rotation axes zero to π radians are sampled using four windowsproviding a DEV value for every π/4 radian (45°) swath (Figure 2).Conversely, isotropic sampling integrates all orientations, computing a local average for the entire neighborhood.This sampling strategy achieves the efficiency required for hyperscaled analysis through the integral image filtering technique, which comes at the cost of coarse orientation resolution.Therefore, measures of the orientation of anisotropy have been excluded from the method due to the coarse resolution orientation sampling.
Figure 2. Schematic demonstrating the theoretical differences between the deviation from mean elevation (DEV) measurements from the π/4 radian swath sampling strategy, an isotropic sampling strategy, and hypothetical continuous sampling strategy using infinitesimally small swaths (effectively infinite line measurements).

Anisotropy Model
The RMSD of the directionally oriented DEV samples (DEVd) was used to characterize anisotropy magnitude (Am), using the isotropically sampled DEV value (DEVi) as the point of comparison in the difference term (Equation (3)).

Anisotropy Model
The RMSD of the directionally oriented DEV samples (DEV d ) was used to characterize anisotropy magnitude (A m ), using the isotropically sampled DEV value (DEV i ) as the point of comparison in the difference term (Equation ( 3)).
where A m is anisotropy magnitude, n is the number of oriented windows, DEV d are DEV values sampled from directionally oriented windows, and DEV i is the DEV value sampled from the isotropic window.Therefore, A m measures variability in DEV across all sampled orientations.A m differs from ratio models used in the semivariogram methods by quantifying the average effect of orientation, as opposed to the maximum difference at an angle, orthogonal or otherwise.It should be noted that A m can be replaced with more sensitive and common metrics such as a ratio or the min-max range using the same sampling strategy.The tools developed for this research are included as plug-in tools in the open-source WhiteboxTools library as part of the Whitebox GAT project [30].The specific details of the algorithm implementation can be viewed in the WhiteboxTools source code.

Analysis
The method described above is designed to measure anisotropy magnitude across both spatial dimensions of a raster digital elevation model (DEM), and across the scale dimension.However, the quantity of information generated threatens to exceed the memory limits of most workstations.Therefore, the interrogation of landscape anisotropy information was implemented in two similar but independent analyses.First, an anisotropy signature was generated for a point location (i.e., a single raster cell), capturing the expression of anisotropy across the full range of scales.This implementation of the method has limited utility beyond preliminary or exploratory analyses for a limited number of locations.Instead, the primary analysis focuses on reducing the scale dimension to a single value that meaningfully characterizes the spatial distribution of hyperscale anisotropy magnitude.While this affords some flexibility in the definition of meaningful, for simplicity this method used the maximum A m value for each cell, similar to the method used by Lindsay et al. (2015) [20] for hyperscale LTP analysis.The A max implementation outputs a second R max raster, which stored the filter half-length at which the A max value occurred.These combined rasters provide a method to feasibly characterize the topographic anisotropy magnitude for a landscape, while also retaining information from the scale dimension.

Data Sets and Study Sites
This study used two fine-resolution data sets containing regularly spaced, similarly sized, and oriented anisotropic landforms to demonstrate the method.Regularly spaced, similarly sized, and oriented anisotropic landforms allow local landscape anisotropy to be spatially associated with landforms known to express anisotropy.The first data set is a 1 m resolution bare-Earth LiDAR DEM from the United States Geological Service's National Elevation Data set and 3D Elevation Program [31], featuring the White Sands National Monument (WSNM) dune field in New Mexico, USA (Figure 3).The desert sand dunes found at WSNM are complex, repeating aeolian dunes that share a similar orientation [32].The WSNM dune field lies to the east of the Alkali Flats and features mainly crescentic dunes in the center, with parabolic dunes around the margin [33].Four individual sites were selected from the WSNM scene as study sites.Site 1 is located in a relatively flat area, absent of dunes or other notable topography.Site 2 is situated on the crest of a small dune, flanked by increasingly large dunes.Due to the crescentic shape of the dunes at Site 2, each dune expresses a large range of orientations.Site 3 is on the crest of a large dune, surrounded by similarly sized and relatively evenly spaced dunes.Site 4 is in between a network of roads, built in the flat areas between similarly sized dunes.The second data set is a 10 m resolution DEM from the Ontario Ministry of Natural Resources' Provincial Digital Elevation Model (v.2) [34] featuring the Peterborough drumlin range (PDR) near Peterborough, Ontario, Canada (Figure 4).Drumlins are anisotropic landforms that are commonly found clustered on postglacial landscapes [35].The PDR contains approximately 3000 drumlins representing a variety of forms, orientations, and elongation ratios [6].The large variation in drumlin scale provides an ideal location to demonstrate the expression of anisotropy magnitude in the landscape texture across scales, in a topographically complex setting.Three individual sites were selected from the PDR scene as study sites.Site 5 is on the crest of a drumlin, surrounded by similarly sized and regularly spaced drumlins.Overall, the drumlins throughout the study site have a northeast-southwest tendency.However, there is substantial local variation in this trend.Along the western edge and the central part of the study site, the orientation of drumlins is more east-west aligned.Towards the eastern-most portions of the PDR the orientation of individual drumlins becomes more northsouth trending.Site 6 is situated on an elevated outcrop that is surrounded on three sides by a floodplain.Site 7 is located in a topographically rough area, with many ridges and valleys, positioned on top of an elevated plateau.The second data set is a 10 m resolution DEM from the Ontario Ministry of Natural Resources' Provincial Digital Elevation Model (v.2) [34] featuring the Peterborough drumlin range (PDR) near Peterborough, Ontario, Canada (Figure 4).Drumlins are anisotropic landforms that are commonly found clustered on postglacial landscapes [35].The PDR contains approximately 3000 drumlins representing a variety of forms, orientations, and elongation ratios [6].The large variation in drumlin scale provides an ideal location to demonstrate the expression of anisotropy magnitude in the landscape texture across scales, in a topographically complex setting.Three individual sites were selected from the PDR scene as study sites.Site 5 is on the crest of a drumlin, surrounded by similarly sized and regularly spaced drumlins.Overall, the drumlins throughout the study site have a northeast-southwest tendency.However, there is substantial local variation in this trend.Along the western edge and the central part of the study site, the orientation of drumlins is more east-west aligned.Towards the eastern-most portions of the PDR the orientation of individual drumlins becomes more north-south trending.Site 6 is situated on an elevated outcrop that is surrounded on three sides by a floodplain.Site 7 is located in a topographically rough area, with many ridges and valleys, positioned on top of an elevated plateau.The second data set is a 10 m resolution DEM from the Ontario Ministry of Natural Resources' Provincial Digital Elevation Model (v.2) [34] featuring the Peterborough drumlin range (PDR) near Peterborough, Ontario, Canada (Figure 4).Drumlins are anisotropic landforms that are commonly found clustered on postglacial landscapes [35].The PDR contains approximately 3000 drumlins representing a variety of forms, orientations, and elongation ratios [6].The large variation in drumlin scale provides an ideal location to demonstrate the expression of anisotropy magnitude in the landscape texture across scales, in a topographically complex setting.Three individual sites were selected from the PDR scene as study sites.Site 5 is on the crest of a drumlin, surrounded by similarly sized and regularly spaced drumlins.Overall, the drumlins throughout the study site have a northeast-southwest tendency.However, there is substantial local variation in this trend.Along the western edge and the central part of the study site, the orientation of drumlins is more east-west aligned.Towards the eastern-most portions of the PDR the orientation of individual drumlins becomes more northsouth trending.Site 6 is situated on an elevated outcrop that is surrounded on three sides by a floodplain.Site 7 is located in a topographically rough area, with many ridges and valleys, positioned on top of an elevated plateau.To improve the visual quality of the outputs, the data sets were smoothed to reduce short-scale roughness in the hillshade images, using the Sun et al. [36] feature-preserving method where edges and corners are unaltered.This also aids the visual interpretation of the results by muting the influence of noise while preserving topographic features, effectively preventing short-scale roughness from dominating the anisotropy signal.In practice, the same effect can be achieved by increasing the minimum filter size.

Topographic Anisotropy-Scale Signatures
Scale signatures are functions of a variable across a range of scales.Therefore, anisotropy-scale signatures characterize how the anisotropy magnitude (A m ) responds to changes in scale, analogous to how spectral signatures characterize surface reflectance response to changes in the wavelength of light.Figure 5 shows scale signatures for the WSNM study sites (Sites 1 through 4) and Figure 6 shows the scale signatures for the PDR study sites (Sites 5 through 7).A m was sampled with filter half-lengths ranging from 7 to 500 cells, incremented by one, for both data sets.This corresponds to filter half-lengths of 7 to 500 m and 70 to 5000 m for the WSNM and PDR data sets, respectively.While a unique scale signature exists for every cell, the seven sites are representative of the varied topographic conditions.Recall that the study areas (as well as study Sites 2, 3, 4, 5, and 6) were chosen for containing similarly sized and oriented, repeating anisotropic landforms.The scale signature for Site 1, a flat area, shows little anisotropy.The signatures for both Sites 3 and 4 also showed minimal anisotropy, despite featuring several linear shaped dunes.The small-scale activity (<90 m) for Site 4 (Figure 6) corresponds with the road network surrounding the location.Conversely, the scale signature for Site 2 captures the proportionally increasing size of the dunes, and the larger dune complex with a multimodal distribution and Sites 5 and 6 show multimodal distributions corresponding to the drumlins and the outcrop respectively, as well as the larger scale topography containing them (i.e., the valley and floodplain).Site 7 featured a topographically rough area, and produced a multimodal signature with relatively small magnitude anisotropy values.
To improve the visual quality of the outputs, the data sets were smoothed to reduce short-scale roughness in the hillshade images, using the Sun et al. [36] feature-preserving method where edges and corners are unaltered.This also aids the visual interpretation of the results by muting the influence of noise while preserving topographic features, effectively preventing short-scale roughness from dominating the anisotropy signal.In practice, the same effect can be achieved by increasing the minimum filter size.

Topographic Anisotropy-Scale Signatures
Scale signatures are functions of a variable across a range of scales.Therefore, anisotropy-scale signatures characterize how the anisotropy magnitude (Am) responds to changes in scale, analogous to how spectral signatures characterize surface reflectance response to changes in the wavelength of light.Figure 5 shows scale signatures for the WSNM study sites (Sites 1 through 4) and Figure 6 shows the scale signatures for the PDR study sites (Sites 5 through 7).Am was sampled with filter half-lengths ranging from 7 to 500 cells, incremented by one, for both data sets.This corresponds to filter half-lengths of 7 to 500 m and 70 to 5000 m for the WSNM and PDR data sets, respectively.While a unique scale signature exists for every cell, the seven sites are representative of the varied topographic conditions.Recall that the study areas (as well as study Sites 2, 3, 4, 5, and 6) were chosen for containing similarly sized and oriented, repeating anisotropic landforms.The scale signature for Site 1, a flat area, shows little anisotropy.The signatures for both Sites 3 and 4 also showed minimal anisotropy, despite featuring several linear shaped dunes.The small-scale activity (<90 m) for Site 4 (Figure 6) corresponds with the road network surrounding the location.Conversely, the scale signature for Site 2 captures the proportionally increasing size of the dunes, and the larger dune complex with a multimodal distribution and Sites 5 and 6 show multimodal distributions corresponding to the drumlins and the outcrop respectively, as well as the larger scale topography containing them (i.e., the valley and floodplain).Site 7 featured a topographically rough area, and produced a multimodal signature with relatively small magnitude anisotropy values.

Spatial Distribution of Topographic Anisotropy
A separate analysis was conducted on each DEM, producing rasters of the maximum anisotropy magnitude value (Amax) for each cell, as well as secondary rasters of the filter half-length at which the maximum anisotropy magnitude value occurred (Rmax).Applying the information from the scale signatures allows the user to avoid the computation of redundant scales that are unlikely to update the Amax value.Combining information from both the Amax and Rmax rasters provides a self-selecting, hyperscale characterization of the spatial distribution of topographic anisotropy in the landscape, with a means to confirm the observed pattern through the Rmax raster.
Amax was calculated for each grid cell of the two DEMs using filter half-lengths ranging from 50 to 300 incremented by 5 cells.This corresponds to filter lengths (in ground units) ranging from 101 to 601 m and 1001 to 6001 m for the WSNM and PDR datasets respectively.For reference, the execution times (excluding data input and output) using an Intel i7-4770 CPU was 6.17 and 1.94 minutes for the WSNM (~116 million cells excluding nodata) and PDR (~40 million cells excluding nodata) data sets respectively.Figure 7 provides examples of the Amax and Rmax raster outputs for both study DEMs, and Figure 8 provides detailed Amax raster examples for the areas surrounding study Sites 2 through 6 to demonstrate the distribution of Amax with respect to local anisotropic textures.The WSNM dune field Rmax raster demonstrates that maximum anisotropy occurred consistently at small scales (<100 m) for the more linear dunes, and occurs at larger scales in the region populated with more crescentic dunes (Figure 7b).The algorithm was also generally able to identify the linear texture created by the dunes, which appear with intermediate and high magnitude anisotropy values (Figures 7a and 8b).This variation in anisotropy magnitude occurred when the dune geometry deviated from a linear pattern to saddle points where dunes merged and smaller scale features with different orientations (e.g., where a dune is transected by a road).
The PDR DEM featured a more complex landscape than WSNM, which was reflected in the highly varied Rmax raster (Figure 7d).Although the PDR Amax and Rmax outputs were more spatially varied compared to the WSNM outputs, the highly anisotropic cells were consistently co-located with anisotropic features such as drumlins, scours, outcrops, valleys and river channels (Figure 8d,e).Relatively larger scale topographic features such as valleys were moderately anisotropic while

Spatial Distribution of Topographic Anisotropy
A separate analysis was conducted on each DEM, producing rasters of the maximum anisotropy magnitude value (A max ) for each cell, as well as secondary rasters of the filter half-length at which the maximum anisotropy magnitude value occurred (R max ).Applying the information from the scale signatures allows the user to avoid the computation of redundant scales that are unlikely to update the A max value.Combining information from both the A max and R max rasters provides a self-selecting, hyperscale characterization of the spatial distribution of topographic anisotropy in the landscape, with a means to confirm the observed pattern through the R max raster.
A max was calculated for each grid cell of the two DEMs using filter half-lengths ranging from 50 to 300 cells, incremented by 5 cells.This corresponds to filter lengths (in ground units) ranging from 101 to 601 m and 1001 to 6001 m for the WSNM and PDR datasets respectively.For reference, the execution times (excluding data input and output) using an Intel i7-4770 CPU was 6.17 and 1.94 minutes for the WSNM (~116 million cells excluding nodata) and PDR (~40 million cells excluding nodata) data sets respectively.Figure 7 provides examples of the A max and R max raster outputs for both study DEMs, and Figure 8 provides detailed A max raster examples for the areas surrounding study Sites 2 through 6 to demonstrate the distribution of A max with respect to local anisotropic textures.The WSNM dune field R max raster demonstrates that maximum anisotropy occurred consistently at small scales (<100 m) for the more linear dunes, and occurs at larger scales in the region populated with more crescentic dunes (Figure 7b).The algorithm was also generally able to identify the linear texture created by the dunes, which appear with intermediate and high magnitude anisotropy values (Figures 7a and 8b).This variation in anisotropy magnitude occurred when the dune geometry deviated from a linear pattern to saddle points where dunes merged and smaller scale features with different orientations (e.g., where a dune is transected by a road).
The PDR DEM featured a more complex landscape than WSNM, which was reflected in the highly varied R max raster (Figure 7d).Although the PDR A max and R max outputs were more spatially varied compared to the WSNM outputs, the highly anisotropic cells were consistently co-located with anisotropic features such as drumlins, scours, outcrops, valleys and river channels (Figure 8d,e).Relatively larger scale topographic features such as valleys were moderately anisotropic while smaller scale features such as drumlins and river channels were more strongly anisotropic (Figure 8d,e).Flat areas were typically isotropic (Figure 7c).Buffering around the edges produced slightly elevated anisotropy values due to the exclusion of larger filter sizes from the analysis.
Geosciences 2018, 8, x FOR PEER REVIEW 9 of 14 smaller scale features such as drumlins and river channels were more strongly anisotropic (Figure 8d,e).Flat areas were typically isotropic (Figure 7c).Buffering around the edges produced slightly elevated anisotropy values due to the exclusion of larger filter sizes from the analysis.

Discussion
The proposed sampling strategy and model implements a computationally efficient hyperscale landscape measure of topographic anisotropy as a scale signature for single cells, and a maximum anisotropy raster.The anisotropy-scale signatures provided a quick method to sample anisotropy values across a range of scales.The scale signatures from all sites produced small anisotropy values at larger scales, corroborating Roy et al.'s finding that landscapes are isotropic at larger spatial scales [1].At smaller scales, the multimodality of the signatures suggests the presence of complex nested anisotropic features (e.g., Figure 5 Site 2, Figure 6 Sites 5 and 6).This is best illustrated by comparing the flat signatures for Sites 1, 3, and 4, with multimodal signatures for Sites 2 and 6.Site 2 featured progressively large dunes as part of an oriented dune complex, and Site 6 featured an outcrop oriented in a different direction from the floodplain.In both cases, nested features appeared as distinct modes in the scale signatures (Figure 5: Site 2, Figure 6: Sites 5 and 6).However, the regularly spaced dunes in Sites 3 and 4 were expected to appear as a series of modes but instead failed to produce any distinct modes over large areas (Figure 5).This is likely an example of mean generalization [37], where the anisotropy signal from the repetitive dunes is muted by the increasing area with which DEV is measured.The multimodal behavior of the scale signatures suggests that the maximum anisotropy analysis can be improved by approximating appropriate scale ranges by each mode, as opposed to using the maximum value across a larger range.
The Amax analysis used the maximum anisotropy values across a range of scales for each cell to create spatially referenced raster images, effectively mapping the highest anisotropy value from every cell's scale signature.Although the multimodal distributions found in the scale signatures suggest that the maximum anisotropy is not the only significant characteristic of signatures (e.g.,

Discussion
The proposed sampling strategy and model implements a computationally efficient hyperscale landscape measure of topographic anisotropy as a scale signature for single cells, and a maximum anisotropy raster.The anisotropy-scale signatures provided a quick method to sample anisotropy values across a range of scales.The scale signatures from all sites produced small anisotropy values at larger scales, corroborating Roy et al.'s finding that landscapes are isotropic at larger spatial scales [1].At smaller scales, the multimodality of the signatures suggests the presence of complex nested anisotropic features (e.g., Figure 5 Site 2, Figure 6 Sites 5 and 6).This is best illustrated by comparing the flat signatures for Sites 1, 3, and 4, with multimodal signatures for Sites 2 and 6.Site 2 featured progressively large dunes as part of an oriented dune complex, and Site 6 featured an outcrop oriented in a different direction from the floodplain.In both cases, nested features appeared as distinct modes in the scale signatures (Figure 5: Site 2, Figure 6: Sites 5 and 6).However, the regularly spaced dunes in Sites 3 and 4 were expected to appear as a series of modes but instead failed to produce any distinct modes over large areas (Figure 5).This is likely an example of mean generalization [37], where the anisotropy signal from the repetitive dunes is muted by the increasing area with which DEV is measured.The multimodal behavior of the scale signatures suggests that the maximum anisotropy analysis can be improved by approximating appropriate scale ranges by each mode, as opposed to using the maximum value across a larger range.
The A max analysis used the maximum anisotropy values across a range of scales for each cell to create spatially referenced raster images, effectively mapping the highest anisotropy value from every cell's scale signature.Although the multimodal distributions found in the scale signatures suggest that the maximum anisotropy is not the only significant characteristic of signatures (e.g., number of peaks indicates signature complexity), high anisotropy values were still consistently co-located with anisotropic landforms.Furthermore, the low spatial variation in the WSNM R max raster (Figure 7b) indicates that the A max values occurred at consistent scales, suggesting that the dune landforms were driving the anisotropy values by influencing the landscape texture.Conversely, the highly varied mosaic of scales in the PDR R max image (Figure 7d) revealed the absence of an optimal single window size with which to measure anisotropy, especially for complex terrain.This also applies to extensive areas where larger scale landscape shaping processes occur simultaneously with smaller scale processes.This emphasizes the value of hyperscale methodological approaches by removing the responsibility of scale selection from the user, and minimizing potential errors as a result.The A max method self-selects the optimal window size for each individual grid cell, resulting in a hyperscale characterization of topographic anisotropy.
This method differs from point [1] or line-based [7] measures of topographic anisotropy by characterizing anisotropy over an area, compared to sampling orientations at discrete intervals.Our approach samples all directions using a small number of oriented sampling areas, at the cost of some information loss.The coarse orientation resolution lacks the ability to sample fine orientation intervals, producing two problems.Firstly, the increased probability that anisotropic features have orientations that exist within multiple windows, which lowers the measured anisotropy by diluting the signal across multiple oriented windows.Close examination of Figure 8a reveals that portions of the crescentic dunes that lie between sampling swaths (e.g., ±22.5 degrees from north) have A max values approximately half as large as other parts of the same dune.This limitation is also expressed as narrow bands of lower A max values at rotation angles that lie in between the swath during the rapid change in orientation of the crest for the crescentic dunes.Secondly, the limited ability to characterize topographic anisotropy for complex geometric patterns.However, the landscape geometries that can be characterized are comparable to the second-order partial derivatives used to define elementary forms [38], morphometric primitives [39], and the Hessian matrices used to detect ridges in the field of computer vision [40,41].The information generated is well-suited as a morphometric characteristic that can augment the translation of continuous landscapes to discrete landforms [38,42,43].Alternatively, this method is also ideal to exploratory analyses, to guide the parameterization (such as scale selection) of targeted measures on known entities.
Texture reveals information about the spatial distribution of variation for a field, and has long been considered a fundamental property of images [44] that also applies to landscapes [12].Trevisani et al. [12] considered surface roughness and anisotropy important elements of textural information for a landscape, revealing the "spatial variability structure" of the area being sampled.Applying a multi-area method to anisotropy sampling extends beyond the local expression of anisotropy by landform objects, and instead characterizing the degree to which the local landscape texture is oriented.The inclusion of anisotropic objects in the search window influences texture for the area.This explains why large anisotropy values often appear in proximity to anisotropic objects.Therefore, this method measures the overall anisotropic texture across areas, not for specific landscape objects.

Conclusions
While measures of anisotropy for specific, delineated objects are well established, continuous measures of topographic anisotropy have only recently been explored.There are a limited number of methodological approaches, all of which feature lower dimension sampling such as points or lines.The method presented here utilizes the computational efficiency of the integral image filtering method to provide a feasible hyperscale, LTP-based measure of anisotropy that relates each point to the landscape containing it.An analysis of the hyperscale behavior of anisotropy through the scale signatures demonstrated that the method can detect complex, hierarchically nested anisotropic topography.The A max and R max rasters demonstrated spatial co-location of large anisotropy magnitude values with anisotropic landforms at scales consistent with the expected scale range of the landforms.Therefore, despite inherently coarse orientation sampling, this method has demonstrated the ability to characterize hyperscale topographic anisotropy in natural landscapes.Future work should compare the RMSD model of anisotropy to the more common ratio models and develop better ways to analyze and visualize voluminous multidimensional spatial data.

Figure 1 .
Figure 1.Graphical representation of the orientated window aggregation scheme.The grey boxes represent the windows from the 3 × 3 matrix that were aggregated to form the (A) northsouth oriented window, (B) eastwest oriented window, (C) northwestsoutheast oriented window, (D) northeastsouthwest oriented window, and (E) the isotropic window.Note that each window from the 3 × 3 matrix is composed of a variable number of grid cells from the input data.

Figure 1 .
Figure 1.Graphical representation of the orientated window aggregation scheme.The grey boxes represent the windows from the 3 × 3 matrix that were aggregated to form the (A) north-south oriented window, (B) east-west oriented window, (C) northwest-southeast oriented window, (D) northeast-southwest oriented window, and (E) the isotropic window.Note that each window from the 3 × 3 matrix is composed of a variable number of grid cells from the input data.

Figure 2 .
Figure 2. Schematic demonstrating the theoretical differences between the deviation from mean elevation (DEV) measurements from the π/4 radian swath sampling strategy, an isotropic sampling strategy, and hypothetical continuous sampling strategy using infinitesimally small swaths (effectively infinite line measurements).

14 Figure 3 .
Figure 3.The White Sands National Monument (WSNM) DEM with the locations of study Sites 1-4.

Figure 4 .
Figure 4.The Peterborough drumlin range (PDR) DEM with the locations of study Sites 5-7.

Figure 3 .
Figure 3.The White Sands National Monument (WSNM) DEM with the locations of study Sites 1-4.

14 Figure 3 .
Figure 3.The White Sands National Monument (WSNM) DEM with the locations of study Sites 1-4.

Figure 4 .
Figure 4.The Peterborough drumlin range (PDR) DEM with the locations of study Sites 5-7.

Figure 4 .
Figure 4.The Peterborough drumlin range (PDR) DEM with the locations of study Sites 5-7.

Figure 5 .
Figure 5. Anisotropy scale signatures for the WSNM data set.The hillshade image for each study site is provided and labelled accordingly.

Figure 5 .
Figure 5. Anisotropy scale signatures for the WSNM data set.The hillshade image for each study site is provided and labelled accordingly.

Figure 6 .
Figure 6.Anisotropy scale signatures for the PDR data set.The hillshade image for each study site is provided and labelled accordingly.

Figure 6 .
Figure 6.Anisotropy scale signatures for the PDR data set.The hillshade image for each study site is provided and labelled accordingly.

Figure 7 .
Figure 7. Spatial distribution of Amax and Rmax rasters for the WNSM (A and B respectively), and PDR (C and D respectively), using filter half-lengths ranging from 50 to 300 grid cells.Both color ramps use blue to represent small magnitude values and red to represent large scale values, relative to the data set.

Figure 7 .
Figure 7. Spatial distribution of A max and R max rasters for the WNSM (A and B respectively), and PDR (C and D respectively), using filter half-lengths ranging from 50 to 300 grid cells.Both color ramps use blue to represent small magnitude values and red to represent large scale values, relative to the data set.

Figure 8 .
Figure 8. High-magnification Amax rasters for the areas surrounding study Sites 2 through 6 using filter size half-lengths ranging from 50 to 300 grid cells.Blue represents isotropy, green and yellow represents intermediate anisotropy, while red represents strong anisotropy, scaled relative to the data set.

Figure 8 .
Figure 8. High-magnification A max rasters for the areas surrounding study Sites 2 through 6 using filter size half-lengths ranging from 50 to 300 grid cells.Blue represents isotropy, green and yellow represents intermediate anisotropy, while red represents strong anisotropy, scaled relative to the data set.