Purpose

This vignette showcases the functions used to derive the roughness and anisotropy exponents from a rotated landscape.

Loading data

First, we perform all the steps from this vignette.

library("statisticalRoughness")
library("rayshader")
library("raster")

gabilan_mesa <- raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "gabilan_mesa.tif"))
gabilan_mesa <- gabilan_mesa %>% detrend_dem() 
raster_resolution <- 9.015
FT2D <- fft2D(raster::as.matrix(gabilan_mesa), dx = raster_resolution, dy = raster_resolution, Hann = TRUE)
nbin <- 20
binned_power_spectrum <- bin(log10(FT2D$radial_frequency_vector), log10(FT2D$spectral_power_vector), nbin)
binned_power_spectrum <- na.omit(binned_power_spectrum)
normalized_spectral_power_matrix <- get_normalized_spectral_power_matrix(binned_power_spectrum, FT2D)
filtered_spectral_power_matrix <- filter_spectral_power_matrix(normalized_spectral_power_matrix, FT2D, quantile_prob = c(0.9999))
ang_fourier <- get_fourier_angle(filtered_spectral_power_matrix, FT2D)
rotated_raster <- rotate_raster(raster::as.matrix(gabilan_mesa), ang_fourier)

On the figure below, the horizontal line is the middle row, the leftmost vertical line is the first third column and the rightmost vertical line is the middle column.

Height-height correlation functions

The height-height correlation function (HHCF) is defined as:

\[C(r) = \left\langle \sqrt{[h(x+r,t)-h(x,t)]^2} \right\rangle\]

with \(h\) the elevation, \(x\) the spatial coordinate, \(t\) the time coordinate, \(r\) a spatial increment and where brackets indicate averaging. With this definition, we can now calculate the HHCF for the rotated_raster in the row-wise (margin = 1) and column-wise directions (margin = 2) using get_hhcf().

hhcf_x <- get_hhcf(rotated_raster, margin = 1)
hhcf_y <- get_hhcf(rotated_raster, margin = 2)

Roughness exponent, \(\alpha\)

Definition and examples

The HHCF scales with increasing \(r\) as:

\[ C(r) \sim r^\alpha \]

with \(\alpha\) the roughness exponent. We now use get_alpha() to derive the roughness exponents \(\alpha_x\) and \(\alpha_y\) in each direction. The option do_plot = TRUE provides a graphical display of the step-wise binned regression process and identification of the cross-over lengthscale. As the regression is a piecewise regression, for each direction, two slopes are returned, for example: \(\alpha_{1,x}\) and \(\alpha_{2,x}\).

get_alpha(hhcf_x$hhcf[ind_x, ], raster_resolution, do_plot = TRUE)

#>        rc    alpha1    alpha2    rmax  alpha.r2
#> x 112.492 0.8514611 0.3053454 243.405 0.9991539

get_alpha(hhcf_y$hhcf[ind_y, ], raster_resolution, do_plot = TRUE)

#>         rc    alpha1    alpha2    rmax  alpha.r2
#> x 69.71757 0.8454504 0.2776742 243.405 0.9993926

get_alpha() includes a number of safeguard and the function will return NA. Here is an example the roughly follow a ridgeline, the middle column, with no statistically consistent scaling.

get_alpha(hhcf_y$hhcf[ind_x, ], raster_resolution, do_plot = TRUE)

#>         rc    alpha1    alpha2   rmax  alpha.r2
#> x 69.15329 0.7871953 0.2112622 126.21 0.9987383

Computing all roughness exponents

Provided a list of HHCF from get_hhcf(), the function get_all_alpha() computes all step-wise regressions.

alpha_x <- get_all_alpha(hhcf_x, raster_resolution)
alpha_y <- get_all_alpha(hhcf_y, raster_resolution)

The resulting data.frames can easily be summarized with summarise_alpha(). For example:

alpha_x %>% summarise_alpha()
#>   alpha1_min alpha1_mean alpha1_max  alpha1_sd alpha1_IQR alpha2_min
#> 1  0.7525691   0.8576394  0.9194875 0.02333023 0.02668279  0.1089438
#>   alpha2_mean alpha2_max  alpha2_sd alpha2_IQR   rc_min  rc_mean   rc_max
#> 1   0.3469332   0.648869 0.08089364  0.1061726 42.56599 117.3809 229.5888
#>      rc_sd   rc_IQR rmax_min rmax_mean rmax_max  rmax_sd rmax_IQR alpha.r2_min
#> 1 25.04692 26.14729    90.15  404.4445   973.62 216.4101    180.3    0.9932127
#>   alpha.r2_mean alpha.r2_max  alpha.r2_sd alpha.r2_IQR
#> 1     0.9989099      0.99965 0.0006989915 0.0004849006
alpha_y %>% summarise_alpha()
#>   alpha1_min alpha1_mean alpha1_max  alpha1_sd alpha1_IQR alpha2_min
#> 1  0.4445179   0.8035647  0.9064168 0.06828347 0.05439905  0.1101542
#>   alpha2_mean alpha2_max alpha2_sd alpha2_IQR   rc_min  rc_mean   rc_max
#> 1   0.3291653   1.002387 0.1200056  0.0899577 49.48651 120.3672 1587.951
#>      rc_sd   rc_IQR rmax_min rmax_mean rmax_max  rmax_sd rmax_IQR alpha.r2_min
#> 1 171.9672 28.32148  117.195  369.5543 2677.455 391.4812  207.345    0.9760342
#>   alpha.r2_mean alpha.r2_max alpha.r2_sd alpha.r2_IQR
#> 1     0.9981882    0.9997217 0.002859781 0.0006231902

In a more human-readable format, we can summarize the entire landscape on rotated_raster by:

	row-wise direction, \(x\)	column-wise direction, \(y\)
\(\alpha_1\)	0.86 \(\pm\) 0.023	0.8 \(\pm\) 0.068
\(\alpha_2\)	0.35 \(\pm\) 0.081	0.33 \(\pm\) 0.12

These numbers mean that the topography is mainly correlated in a positive way below NA meters and more correlated in the \(x\)-direction than in the \(y\)-direction as \(\alpha_{1,x} > \alpha_{1,y} > 0.5\). Above NA meters, the topography is negatively correlated and more so in the y-direction than in the \(x\)-direction \(\alpha_{2,y} < \alpha_{2,x} < 0.5\). These results also mean that the \(x\)-direction corresponds to the down-slope direction, \(\mathbf{e}_\parallel\), while the \(y\)-direction corresponds to the across-slope direction, \(\mathbf{e}_\perp\). Further, the dominant slope in the \(x\)-direction results in increased smoothing and higher values of roughness exponents. This is what we could have expected from a visual inspection of the landscape.

Anisotropy exponent, \(\zeta\)

The anisotropy exponent is the ratio between the down-slope and across-slope roughness exponents:

\[ \zeta = \frac{\alpha_\perp}{\alpha_\parallel}\]

Because of the increased smoothing due to gravity in the down-slope direction, \(\alpha_\perp \geq \alpha_\parallel\) and \(\zeta \geq 1\). From the central values of \(\alpha_{1,x}\) and \(\alpha_{1,y}\), the value of the anisotropy exponent \(\zeta_1\) can be derived as:

\[ \zeta_1 = \frac{\max{\alpha_{1,x}, \alpha_{1,x}}}{\min{\alpha_{1,x}, \alpha_{1,x}}}\]

This estimate of \(\zeta\) is handled by an internal function which includes a check for the spread of the distribution of \(\alpha\). If the distribution is not constrained enough, NA values are returned. For the current landscape:

	anisotropy exponent, \(\zeta\)
\(\zeta_1\)	1.07
\(\zeta_2\)	1.05

Deriving the anisotropy exponent

Hervé Guillon