In hrvg/statisticalRoughness: Estimation of roughness exponents across scale for large areas

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8, 
  fig.height = 8,
  warning = FALSE,
  message = FALSE
)

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")

rstr <- raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "submarine_canyon.tif"))
# raster_resolution <- 10
raster_resolution <- 50
# rstr <- raster::aggregate(rstr, fact = 9) # 90 m to match Pastor-Satorras
# raster_resolution <- 90
FT2D <- rstr %>% detrend_dem() %>% raster::as.matrix() %>%
  fft2D(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)
ang_fourier <- 0
rotated_raster <- rotate_raster(rstr, 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.

ind_x <- nrow(rotated_raster) %/% 2
ind_y <- ncol(rotated_raster) %/% 3
.rotated_raster <- rotated_raster
.rotated_raster[ind_x, ] <- NA
.rotated_raster[, ind_y] <- NA
.rotated_raster[, ind_x] <- NA

matrix_for_rayshader(.rotated_raster) %>%
    sphere_shade(texture = "imhof4") %>%
  add_shadow(ray_shade(matrix_for_rayshader(.rotated_raster), sunaltitude = 20), max_darken = 0.3) %>%  
  add_shadow(ambient_shade(matrix_for_rayshader(.rotated_raster)), 0) %>%
  plot_map()

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, raster_resolution, margin = 1)
hhcf_y <- get_hhcf_(rotated_raster, raster_resolution, 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, hhcf_x$autocorr_len[ind_x], do_plot = TRUE)

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

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, hhcf_y$autocorr_len[ind_x], do_plot = TRUE)

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)

alphas <- rbind(
  alpha_x %>% stats::na.omit() %>% dplyr::mutate(direction = 'x'),
  alpha_y %>% stats::na.omit() %>% dplyr::mutate(direction = 'y')
  ) %>% reshape2::melt(measure.vars = c("alpha1", "alpha2"))
p <- ggplot2::ggplot(alphas, ggplot2::aes(x = .data$value, group = .data$direction, color = .data$direction)) +
  ggplot2::geom_density(n = 512) +
  ggplot2::facet_grid(.data$variable ~ .) +
  ggpubr::theme_pubr()
p

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

alpha_x %>% summarise_alpha()
alpha_y %>% summarise_alpha()

x_stats <- alpha_x %>% summarise_alpha()
y_stats <- alpha_y %>% summarise_alpha()

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$ | r x_stats$alpha1_mean %>% signif(2) $\pm$ r x_stats$alpha1_sd %>% signif(2) | r y_stats$alpha1_mean %>% signif(2) $\pm$ r y_stats$alpha1_sd %>% signif(2) | | $\alpha_2$ | r x_stats$alpha2_mean %>% signif(2) $\pm$ r x_stats$alpha2_sd %>% signif(2) | r y_stats$alpha2_mean %>% signif(2) $\pm$ r y_stats$alpha2_sd %>% signif(2) |

These numbers mean that the topography is mainly correlated in a positive way below r mean(c(x_stats$change_point_mean, y_stats$change_point_mean)) %>% signif(3) meters and more correlated in the $x$-direction than in the $y$-direction as $\alpha_{1,x} > \alpha_{1,y} > 0.5$. Above r mean(c(x_stats$change_point_mean, y_stats$change_point_mean)) %>% signif(3) 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$ | r get_zeta_(x_stats$alpha1_mean, y_stats$alpha1_mean, x_stats$alpha1_IQR, y_stats$alpha1_IQR, TRUE) %>% signifNA() | | $\zeta_2$ | r get_zeta_(x_stats$alpha2_mean, y_stats$alpha2_mean, x_stats$alpha2_IQR, y_stats$alpha2_IQR, TRUE) %>% signifNA() |

get_zeta(rstr, raster_resolution, full = TRUE) %>% dplyr::mutate_all(signifNA)

get_zeta(rstr,raster_resolution) %>% dplyr::mutate_all(signifNA)

hrvg/statisticalRoughness documentation built on March 12, 2021, 4:55 p.m.