In hrvg/statisticalRoughness: Estimation of roughness exponents across scale for large areas
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 7,
fig.height = 7,
warning = FALSE,
message = FALSE
)
Purpose
This vignette displays the calculation of anisotropy exponent for different landscapes: Gabilan Mesa, Yosemite Valley and Modoc Plateau.
library("statisticalRoughness")
library("rayshader")
view_rotated_raster <- function(rstr, raster_resolution = 9.015, nbin = 50, .Hann = TRUE, .quantile_prob = c(0.9999)){
rstr <- rstr %>% detrend_dem()
FT2D <- fft2D(raster::as.matrix(rstr), dx = raster_resolution, dy = raster_resolution, Hann = .Hann)
binned_power_spectrum <- bin(log10(FT2D$radial_frequency_vector), log10(FT2D$spectral_power_vector), nbin) %>% na.omit()
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 = .quantile_prob)
ang_fourier <- get_fourier_angle(filtered_spectral_power_matrix, FT2D)
rotated_raster <- rotate_raster(raster::as.matrix(rstr), ang_fourier)
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()
}
Gabilan Mesa
Gabilan Mesa landscape is the same landscape that we used in other vignettes and the test bed of @Perron2008 methods.
Using get_zeta() returns a data.frame with the variable of interest for the analysis of topographic roughness and anisotropy.
raster_resolution <- 10
gabilan_mesa <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "gabilan_mesa.tif"))
res_gabi <- get_zeta(gabilan_mesa, raster_resolution) %>% dplyr::mutate_all(signifNA)
Gabilan Mesa landscape has
a minimum elevation of r raster::minValue(gabilan_mesa) m,
an average elevation of r raster::cellStats(gabilan_mesa, stat = "mean") m,
and a maximum elevation of r raster::maxValue(gabilan_mesa) m.
The maximum relief is thus r (raster::maxValue(gabilan_mesa) - raster::minValue(gabilan_mesa)) m.
The figure below showcases the landscape in its original form.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_gabi$theta %% 180$^\circ$ and
$\theta_y =$ r (res_gabi$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_gabi$beta1, and $\beta_2 =$ r res_gabi$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_gabi$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_gabi$beta1) - 2)/2), and
r ((abs(res_gabi$beta2) - 2)/2), respectively.
Before a cut-off lengthscale of
r res_gabi$rc m,
the roughness exponents in the $x$-- and $y$--directions are
$\alpha_{1,x} =$ r res_gabi$alpha1.x, and
$\alpha_{1,y} =$ r res_gabi$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_gabi$zeta1.
After cut-off lengthscale,
the roughness exponents are lower and almost equal with
$\alpha_{2,x} =$ r res_gabi$alpha2.x, and
$\alpha_{2,y} =$ r res_gabi$alpha2.y.
Consequently, the anisotropy exponent is $\zeta_2 =$ r res_gabi$zeta2.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_gabi$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2))
matrix_for_rayshader(raster::as.matrix(gabilan_mesa)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(gabilan_mesa)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(gabilan_mesa))), 0) %>%
plot_map()
gabilan_mesa <- gabilan_mesa %>% spatialEco::raster.gaussian.smooth(
sigma = res_gabi$rc %/% raster_resolution
)
matrix_for_rayshader(raster::as.matrix(gabilan_mesa)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(gabilan_mesa)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(gabilan_mesa))), 0) %>%
plot_map()
Yosemite Valley
raster_resolution <- 10
yosemite <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "yosemite.tif"))
res_yose <- get_zeta(yosemite, raster_resolution) %>% dplyr::mutate_all(signifNA)
Yosemite Valley landscape has
a minimum elevation of r raster::minValue(yosemite) m,
an average elevation of r raster::cellStats(yosemite, stat = "mean") m,
and a maximum elevation of r raster::maxValue(yosemite) m.
The maximum relief is thus r (raster::maxValue(yosemite) - raster::minValue(yosemite)) m.
The figure below showcases the landscape in its original form.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_yose$theta %% 180$^\circ$ and
$\theta_y =$ r (res_yose$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_yose$beta1, and $\beta_2 =$ r res_yose$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_yose$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_yose$beta1) - 2)/2), and
r ((abs(res_yose$beta2) - 2)/2).
Before a cut-off lengthscale of
r res_yose$rc m,
the roughness exponents in the $x$-- and $y$--directions are
$\alpha_{1,x} =$ r res_yose$alpha1.x, and
$\alpha_{1,y} =$ r res_yose$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_yose$zeta1.
After cut-off lengthscale,
the roughness exponents are lower and almost equal with
$\alpha_{2,x} =$ r res_yose$alpha2.x, and
$\alpha_{2,y} =$ r res_yose$alpha2.y.
Consequently, the anisotropy exponent is $\zeta_2 =$ r res_yose$zeta2.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_yose$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2))
matrix_for_rayshader(raster::as.matrix(yosemite)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(yosemite)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(yosemite))), 0) %>%
plot_map()
yosemite <- yosemite %>% spatialEco::raster.gaussian.smooth(
sigma = res_yose$rc %/% raster_resolution
)
matrix_for_rayshader(raster::as.matrix(yosemite)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(yosemite)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(yosemite))), 0) %>%
plot_map()
Modoc Plateau
raster_resolution <- 10
modoc <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "modoc.tif"))
res_modo <- get_zeta(modoc, raster_resolution) %>% dplyr::mutate_all(signifNA)
Modoc Plateau landscape has
a minimum elevation of r raster::minValue(modoc) m,
an average elevation of r raster::cellStats(modoc, stat = "mean") m,
and a maximum elevation of r raster::maxValue(modoc) m.
The maximum relief is thus r (raster::maxValue(modoc) - raster::minValue(modoc)) m.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_modo$theta %% 180$^\circ$ and
$\theta_y =$ r (res_modo$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_modo$beta1, and $\beta_2 =$ r res_modo$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_modo$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_modo$beta1) - 2)/2), and
r ((abs(res_modo$beta2) - 2)/2).
Before a cut-off lengthscale of
r res_modo$rc m,
the roughness exponents in the $x$-- and $y$--directions are almost equal as
$\alpha_{1,x} =$ r res_modo$alpha1.x, and
$\alpha_{1,y} =$ r res_modo$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_modo$zeta1.
After cut-off lengthscale,
the roughness exponents are more distinct as
$\alpha_{2,x} =$ r res_modo$alpha2.x, and
$\alpha_{2,y} =$ r res_modo$alpha2.y.
However, the anisotropy exponent, $\zeta_2$, cannot be derived as at least one ofthe distribution of $\alpha_2$ has a large spread and its mean is an inconclusive description of centrality.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_modo$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2))
matrix_for_rayshader(raster::as.matrix(modoc)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(modoc)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(modoc))), 0) %>%
plot_map()
modoc <- modoc %>% spatialEco::raster.gaussian.smooth(
sigma = res_modo$rc %/% raster_resolution
)
matrix_for_rayshader(raster::as.matrix(modoc)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(modoc)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(modoc))), 0) %>%
plot_map()
Marble Canyon
raster_resolution <- 90
marble_canyon <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "marble_canyon.tif"))
marble_canyon <- raster::aggregate(marble_canyon, fact = 9) # 90 m to match Pastor-Satorras
res_mrbl <- get_zeta(marble_canyon, raster_resolution) %>% dplyr::mutate_all(signifNA)
Marble Canyon landscape has
a minimum elevation of r raster::minValue(marble_canyon) m,
an average elevation of r raster::cellStats(marble_canyon, stat = "mean") m,
and a maximum elevation of r raster::maxValue(marble_canyon) m.
The maximum relief is thus r (raster::maxValue(marble_canyon) - raster::minValue(marble_canyon)) m.
The figure below showcases the landscape in its original form.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_mrbl$theta %% 180$^\circ$ and
$\theta_y =$ r (res_mrbl$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_mrbl$beta1, and $\beta_2 =$ r res_mrbl$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_mrbl$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_mrbl$beta1) - 2)/2), and
r ((abs(res_mrbl$beta2) - 2)/2), respectively.
Before a cut-off lengthscale of
r res_mrbl$rc m,
the roughness exponents in the $x$-- and $y$--directions are
$\alpha_{1,x} =$ r res_mrbl$alpha1.x, and
$\alpha_{1,y} =$ r res_mrbl$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_mrbl$zeta1.
After cut-off lengthscale,
the roughness exponents are lower and almost equal with
$\alpha_{2,x} =$ r res_mrbl$alpha2.x, and
$\alpha_{2,y} =$ r res_mrbl$alpha2.y.
Consequently, the anisotropy exponent is $\zeta_2 =$ r res_mrbl$zeta2.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_mrbl$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2))
matrix_for_rayshader(raster::as.matrix(marble_canyon)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(marble_canyon)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(marble_canyon))), 0) %>%
plot_map()
marble_canyon <- marble_canyon %>% spatialEco::raster.gaussian.smooth(
sigma = res_gabi$rc %/% raster_resolution
)
matrix_for_rayshader(raster::as.matrix(marble_canyon)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(marble_canyon)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(marble_canyon))), 0) %>%
plot_map()
Submarine Canyon
raster_resolution <- 50
submarine_canyon <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "submarine_canyon.tif"))
res_subm <- get_zeta(submarine_canyon, raster_resolution) %>% dplyr::mutate_all(signifNA)
Submarine Canyon landscape has
a minimum elevation of r raster::minValue(submarine_canyon) m,
an average elevation of r raster::cellStats(submarine_canyon, stat = "mean") m,
and a maximum elevation of r raster::maxValue(submarine_canyon) m.
The maximum relief is thus r (raster::maxValue(submarine_canyon) - raster::minValue(submarine_canyon)) m.
The figure below showcases the landscape in its original form.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_subm$theta %% 180$^\circ$ and
$\theta_y =$ r (res_subm$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_subm$beta1, and $\beta_2 =$ r res_subm$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_subm$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_subm$beta1) - 2)/2), and
r ((abs(res_subm$beta2) - 2)/2), respectively.
Before a cut-off lengthscale of
r res_subm$rc m,
the roughness exponents in the $x$-- and $y$--directions are
$\alpha_{1,x} =$ r res_subm$alpha1.x, and
$\alpha_{1,y} =$ r res_subm$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_subm$zeta1.
After cut-off lengthscale,
the roughness exponents are lower and almost equal with
$\alpha_{2,x} =$ r res_subm$alpha2.x, and
$\alpha_{2,y} =$ r res_subm$alpha2.y.
Consequently, the anisotropy exponent is $\zeta_2 =$ r res_subm$zeta2.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_subm$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2))
matrix_for_rayshader(raster::as.matrix(submarine_canyon)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(submarine_canyon)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(submarine_canyon))), 0) %>%
plot_map()
submarine_canyon <- submarine_canyon %>% spatialEco::raster.gaussian.smooth(
sigma = res_gabi$rc %/% raster_resolution
)
matrix_for_rayshader(raster::as.matrix(submarine_canyon)) %>%
sphere_shade(texture = "imhof4") %>%
add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(submarine_canyon)), sunaltitude = 20), max_darken = 0.3) %>%
add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(submarine_canyon))), 0) %>%
plot_map()
Comparing Landscapes
The following table summarizes the metrics for the three landscapes.
res <- rbind(res_gabi, res_yose, res_modo, res_mrbl, res_subm) %>%
dplyr::select(-c("w", "xi", "alpha1", "alpha2"))
rownames(res) <- c("Gabilan Mesa", "Yosemite Valley", "Modoc Plateau", "Marble Canyon", "Submarine Canyon")
cnames <- c(
"$\\beta_1$",
"$\\beta_2$",
"$\\alpha_{1,x}$",
"$\\alpha_{1,y}$",
"$\\zeta_1$",
"$\\alpha_{2,x}$",
"$\\alpha_{2,y}$",
"$\\zeta_2$",
"$\\theta$",
"$1/f_c$",
"$r_c$",
"$\\xi_x$",
"$\\xi_y$",
"$w_x$",
"$w_y$"
)
res %>% knitr::kable("html", col.names = cnames) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),
position = "center",
full_width = FALSE) %>%
kableExtra::scroll_box(width = "7in", height = "3in")
While the estimate of statistical roughness from $\beta$ and $\alpha$ values differ for the three landscapes, their interpretations are coherent.
Gabilan Mesa landscape is weakly anisotropic and mainly correlated below r res_gabi$rc m, and isotropic and mainly anti-correlated above r res_gabi$rc m.
Yosemite Valley landscape is anisotropic and mainly correlated below r res_yose$rc m, and isotropic with no clear correlation above r res_yose$rc m.
Modoc Plateau landscape is isotropic and mainly correlated below r res_modo$rc m, and with undetermined behaviour above r res_modo$rc m.
References
hrvg/statisticalRoughness documentation built on March 12, 2021, 4:55 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 7, warning = FALSE, message = FALSE )
Purpose
This vignette displays the calculation of anisotropy exponent for different landscapes: Gabilan Mesa, Yosemite Valley and Modoc Plateau.
library("statisticalRoughness") library("rayshader") view_rotated_raster <- function(rstr, raster_resolution = 9.015, nbin = 50, .Hann = TRUE, .quantile_prob = c(0.9999)){ rstr <- rstr %>% detrend_dem() FT2D <- fft2D(raster::as.matrix(rstr), dx = raster_resolution, dy = raster_resolution, Hann = .Hann) binned_power_spectrum <- bin(log10(FT2D$radial_frequency_vector), log10(FT2D$spectral_power_vector), nbin) %>% na.omit() 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 = .quantile_prob) ang_fourier <- get_fourier_angle(filtered_spectral_power_matrix, FT2D) rotated_raster <- rotate_raster(raster::as.matrix(rstr), ang_fourier) 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() }
Gabilan Mesa
Gabilan Mesa landscape is the same landscape that we used in other vignettes and the test bed of @Perron2008 methods.
Using get_zeta() returns a data.frame with the variable of interest for the analysis of topographic roughness and anisotropy.
raster_resolution <- 10 gabilan_mesa <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "gabilan_mesa.tif")) res_gabi <- get_zeta(gabilan_mesa, raster_resolution) %>% dplyr::mutate_all(signifNA)
Gabilan Mesa landscape has
a minimum elevation of r raster::minValue(gabilan_mesa) m,
an average elevation of r raster::cellStats(gabilan_mesa, stat = "mean") m,
and a maximum elevation of r raster::maxValue(gabilan_mesa) m.
The maximum relief is thus r (raster::maxValue(gabilan_mesa) - raster::minValue(gabilan_mesa)) m.
The figure below showcases the landscape in its original form.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_gabi$theta %% 180$^\circ$ and
$\theta_y =$ r (res_gabi$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_gabi$beta1, and $\beta_2 =$ r res_gabi$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_gabi$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_gabi$beta1) - 2)/2), and
r ((abs(res_gabi$beta2) - 2)/2), respectively.
Before a cut-off lengthscale of
r res_gabi$rc m,
the roughness exponents in the $x$-- and $y$--directions are
$\alpha_{1,x} =$ r res_gabi$alpha1.x, and
$\alpha_{1,y} =$ r res_gabi$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_gabi$zeta1.
After cut-off lengthscale,
the roughness exponents are lower and almost equal with
$\alpha_{2,x} =$ r res_gabi$alpha2.x, and
$\alpha_{2,y} =$ r res_gabi$alpha2.y.
Consequently, the anisotropy exponent is $\zeta_2 =$ r res_gabi$zeta2.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_gabi$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2)) matrix_for_rayshader(raster::as.matrix(gabilan_mesa)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(gabilan_mesa)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(gabilan_mesa))), 0) %>% plot_map() gabilan_mesa <- gabilan_mesa %>% spatialEco::raster.gaussian.smooth( sigma = res_gabi$rc %/% raster_resolution ) matrix_for_rayshader(raster::as.matrix(gabilan_mesa)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(gabilan_mesa)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(gabilan_mesa))), 0) %>% plot_map()
Yosemite Valley
raster_resolution <- 10 yosemite <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "yosemite.tif")) res_yose <- get_zeta(yosemite, raster_resolution) %>% dplyr::mutate_all(signifNA)
Yosemite Valley landscape has
a minimum elevation of r raster::minValue(yosemite) m,
an average elevation of r raster::cellStats(yosemite, stat = "mean") m,
and a maximum elevation of r raster::maxValue(yosemite) m.
The maximum relief is thus r (raster::maxValue(yosemite) - raster::minValue(yosemite)) m.
The figure below showcases the landscape in its original form.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_yose$theta %% 180$^\circ$ and
$\theta_y =$ r (res_yose$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_yose$beta1, and $\beta_2 =$ r res_yose$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_yose$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_yose$beta1) - 2)/2), and
r ((abs(res_yose$beta2) - 2)/2).
Before a cut-off lengthscale of
r res_yose$rc m,
the roughness exponents in the $x$-- and $y$--directions are
$\alpha_{1,x} =$ r res_yose$alpha1.x, and
$\alpha_{1,y} =$ r res_yose$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_yose$zeta1.
After cut-off lengthscale,
the roughness exponents are lower and almost equal with
$\alpha_{2,x} =$ r res_yose$alpha2.x, and
$\alpha_{2,y} =$ r res_yose$alpha2.y.
Consequently, the anisotropy exponent is $\zeta_2 =$ r res_yose$zeta2.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_yose$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2)) matrix_for_rayshader(raster::as.matrix(yosemite)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(yosemite)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(yosemite))), 0) %>% plot_map() yosemite <- yosemite %>% spatialEco::raster.gaussian.smooth( sigma = res_yose$rc %/% raster_resolution ) matrix_for_rayshader(raster::as.matrix(yosemite)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(yosemite)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(yosemite))), 0) %>% plot_map()
Modoc Plateau
raster_resolution <- 10 modoc <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "modoc.tif")) res_modo <- get_zeta(modoc, raster_resolution) %>% dplyr::mutate_all(signifNA)
Modoc Plateau landscape has
a minimum elevation of r raster::minValue(modoc) m,
an average elevation of r raster::cellStats(modoc, stat = "mean") m,
and a maximum elevation of r raster::maxValue(modoc) m.
The maximum relief is thus r (raster::maxValue(modoc) - raster::minValue(modoc)) m.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_modo$theta %% 180$^\circ$ and
$\theta_y =$ r (res_modo$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_modo$beta1, and $\beta_2 =$ r res_modo$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_modo$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_modo$beta1) - 2)/2), and
r ((abs(res_modo$beta2) - 2)/2).
Before a cut-off lengthscale of
r res_modo$rc m,
the roughness exponents in the $x$-- and $y$--directions are almost equal as
$\alpha_{1,x} =$ r res_modo$alpha1.x, and
$\alpha_{1,y} =$ r res_modo$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_modo$zeta1.
After cut-off lengthscale,
the roughness exponents are more distinct as
$\alpha_{2,x} =$ r res_modo$alpha2.x, and
$\alpha_{2,y} =$ r res_modo$alpha2.y.
However, the anisotropy exponent, $\zeta_2$, cannot be derived as at least one ofthe distribution of $\alpha_2$ has a large spread and its mean is an inconclusive description of centrality.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_modo$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2)) matrix_for_rayshader(raster::as.matrix(modoc)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(modoc)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(modoc))), 0) %>% plot_map() modoc <- modoc %>% spatialEco::raster.gaussian.smooth( sigma = res_modo$rc %/% raster_resolution ) matrix_for_rayshader(raster::as.matrix(modoc)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(modoc)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(modoc))), 0) %>% plot_map()
Marble Canyon
raster_resolution <- 90 marble_canyon <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "marble_canyon.tif")) marble_canyon <- raster::aggregate(marble_canyon, fact = 9) # 90 m to match Pastor-Satorras res_mrbl <- get_zeta(marble_canyon, raster_resolution) %>% dplyr::mutate_all(signifNA)
Marble Canyon landscape has
a minimum elevation of r raster::minValue(marble_canyon) m,
an average elevation of r raster::cellStats(marble_canyon, stat = "mean") m,
and a maximum elevation of r raster::maxValue(marble_canyon) m.
The maximum relief is thus r (raster::maxValue(marble_canyon) - raster::minValue(marble_canyon)) m.
The figure below showcases the landscape in its original form.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_mrbl$theta %% 180$^\circ$ and
$\theta_y =$ r (res_mrbl$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_mrbl$beta1, and $\beta_2 =$ r res_mrbl$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_mrbl$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_mrbl$beta1) - 2)/2), and
r ((abs(res_mrbl$beta2) - 2)/2), respectively.
Before a cut-off lengthscale of
r res_mrbl$rc m,
the roughness exponents in the $x$-- and $y$--directions are
$\alpha_{1,x} =$ r res_mrbl$alpha1.x, and
$\alpha_{1,y} =$ r res_mrbl$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_mrbl$zeta1.
After cut-off lengthscale,
the roughness exponents are lower and almost equal with
$\alpha_{2,x} =$ r res_mrbl$alpha2.x, and
$\alpha_{2,y} =$ r res_mrbl$alpha2.y.
Consequently, the anisotropy exponent is $\zeta_2 =$ r res_mrbl$zeta2.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_mrbl$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2)) matrix_for_rayshader(raster::as.matrix(marble_canyon)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(marble_canyon)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(marble_canyon))), 0) %>% plot_map() marble_canyon <- marble_canyon %>% spatialEco::raster.gaussian.smooth( sigma = res_gabi$rc %/% raster_resolution ) matrix_for_rayshader(raster::as.matrix(marble_canyon)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(marble_canyon)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(marble_canyon))), 0) %>% plot_map()
Submarine Canyon
raster_resolution <- 50 submarine_canyon <- raster::raster(file.path(system.file("extdata/rasters/", package = "statisticalRoughness"), "submarine_canyon.tif")) res_subm <- get_zeta(submarine_canyon, raster_resolution) %>% dplyr::mutate_all(signifNA)
Submarine Canyon landscape has
a minimum elevation of r raster::minValue(submarine_canyon) m,
an average elevation of r raster::cellStats(submarine_canyon, stat = "mean") m,
and a maximum elevation of r raster::maxValue(submarine_canyon) m.
The maximum relief is thus r (raster::maxValue(submarine_canyon) - raster::minValue(submarine_canyon)) m.
The figure below showcases the landscape in its original form.
The landscape was rotated so that in this original frame of reference,
$\theta_x =$ r res_subm$theta %% 180$^\circ$ and
$\theta_y =$ r (res_subm$theta %% 180 + 90)$^\circ$.
(In the rotated frame of reference, $\theta_x' = 0^\circ$ and $\theta_y' = 90^\circ$.)
The segmented regression on the power spectrum identifies scalings with slope $\beta_1 =$ r res_subm$beta1, and $\beta_2 =$ r res_subm$beta2, before and after a cut-off lengthscale of $1/f_c =$ r (res_subm$inv.fc) m, respectively.
As $\beta = 2H + d$, and $d = 2$, the Hurst coefficients $H$ corresponding to $\beta_1$ and $\beta_2$ are
r ((abs(res_subm$beta1) - 2)/2), and
r ((abs(res_subm$beta2) - 2)/2), respectively.
Before a cut-off lengthscale of
r res_subm$rc m,
the roughness exponents in the $x$-- and $y$--directions are
$\alpha_{1,x} =$ r res_subm$alpha1.x, and
$\alpha_{1,y} =$ r res_subm$alpha1.y.
Consequently, the anisotropy exponent is $\zeta_1 =$ r res_subm$zeta1.
After cut-off lengthscale,
the roughness exponents are lower and almost equal with
$\alpha_{2,x} =$ r res_subm$alpha2.x, and
$\alpha_{2,y} =$ r res_subm$alpha2.y.
Consequently, the anisotropy exponent is $\zeta_2 =$ r res_subm$zeta2.
The figure below showcases the landscape in its original form (left) and smoothed with a Gaussian kernel of r res_subm$rc %/% raster_resolution standard deviations to match the value of $r_c$ (right).
par(mfrow = c(1,2)) matrix_for_rayshader(raster::as.matrix(submarine_canyon)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(submarine_canyon)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(submarine_canyon))), 0) %>% plot_map() submarine_canyon <- submarine_canyon %>% spatialEco::raster.gaussian.smooth( sigma = res_gabi$rc %/% raster_resolution ) matrix_for_rayshader(raster::as.matrix(submarine_canyon)) %>% sphere_shade(texture = "imhof4") %>% add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(submarine_canyon)), sunaltitude = 20), max_darken = 0.3) %>% add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(submarine_canyon))), 0) %>% plot_map()
Comparing Landscapes
The following table summarizes the metrics for the three landscapes.
res <- rbind(res_gabi, res_yose, res_modo, res_mrbl, res_subm) %>% dplyr::select(-c("w", "xi", "alpha1", "alpha2")) rownames(res) <- c("Gabilan Mesa", "Yosemite Valley", "Modoc Plateau", "Marble Canyon", "Submarine Canyon")
cnames <- c( "$\\beta_1$", "$\\beta_2$", "$\\alpha_{1,x}$", "$\\alpha_{1,y}$", "$\\zeta_1$", "$\\alpha_{2,x}$", "$\\alpha_{2,y}$", "$\\zeta_2$", "$\\theta$", "$1/f_c$", "$r_c$", "$\\xi_x$", "$\\xi_y$", "$w_x$", "$w_y$" ) res %>% knitr::kable("html", col.names = cnames) %>% kableExtra::kable_styling(bootstrap_options = c("striped", "hover"), position = "center", full_width = FALSE) %>% kableExtra::scroll_box(width = "7in", height = "3in")
While the estimate of statistical roughness from $\beta$ and $\alpha$ values differ for the three landscapes, their interpretations are coherent.
Gabilan Mesa landscape is weakly anisotropic and mainly correlated below r res_gabi$rc m, and isotropic and mainly anti-correlated above r res_gabi$rc m.
Yosemite Valley landscape is anisotropic and mainly correlated below r res_yose$rc m, and isotropic with no clear correlation above r res_yose$rc m.
Modoc Plateau landscape is isotropic and mainly correlated below r res_modo$rc m, and with undetermined behaviour above r res_modo$rc m.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.