zeta.Rmd
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("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.
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()
.
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
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.
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 |