# Purpose

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

First, we perform all the steps from this vignette.

library("statisticalRoughness")
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