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
)

r pretty_region_name(region)

results_directory <- file.path(system.file(paste0("extdata/zeta_results/", region), package = "statisticalRoughness"))
zeta_results <- read_zeta_raster(out_dir = results_directory, Lmax = 1E3, raster_resolution = 10, .len = 10, proxy = FALSE)
rstr <- raster::raster(file.path(system.file("extdata/big_rasters/", package = "statisticalRoughness"), paste0(region, ".tif")))

Lmax <- min(head(dim(rstr), 2))
spatial_scales <- get_all_R_L(Lmax, 5, len = 13)$allL %>% head(-1)
spatial_scales <- spatial_scales[spatial_scales >= 30]
zeta_results$spatial_scales <- tail(spatial_scales, length(zeta_results$raster_list)) * 10
zeta_results$spatial_scales

band_names <- read.csv(file.path(results_directory, "band_names.csv")) %>% names()

if ("theta" %in% .selected) zeta_results$raster_list <- zeta_results$raster_list %>% make_angular(match("theta", band_names))

DEM

rstr <- raster::aggregate(rstr, fact = 10)
library(rayshader)
matrix_for_rayshader(raster::as.matrix(rstr)) %>%
  sphere_shade(colorintensity = 0.5, texture = "imhof4") %>%
    add_shadow(ray_shade(matrix_for_rayshader(raster::as.matrix(rstr)), sunaltitude = 20), max_darken = 0.1) %>%  
    # add_shadow(ambient_shade(matrix_for_rayshader(raster::as.matrix(rstr))), 0) %>%
    plot_map()

Maps - clamped - denoised {.tabset}

sliced_clamped_raster <- slice_clamp(
  raster_list = zeta_results$raster_list,
  att_names = band_names,
  selected = .selected,
  clamp_raster = TRUE, 
  clamp_params = cparams)

sliced_clamped_raster <- reduce_spatial_noise(sliced_clamped_raster)

#' @param raster_list a `list` of `stars` objects
#' @param target_id the id of the layer to be masked
#' @return a `list` of `stars` 
logtransform_layer <- function(raster_list, target_id){
  transformed_raster_list <- lapply(seq_along(raster_list), function(n){
    s <- raster_list[[n]] %>% as("Raster")
    transformed_values <- raster::getValues(s[[target_id]])
    transformed_values <- log10(transformed_values)
    transformed_values[!is.finite(transformed_values)] <- NA
    s <- raster::setValues(s, transformed_values, layer = target_id)
    return(stars::st_as_stars(s))
  })
  return(transformed_raster_list)
}

.sliced_clamped_raster <- sliced_clamped_raster
for (x in c("alpha1.y", "alpha1.x", "alpha2.y", "alpha2.x", "w.y", "w.x", "xi.y", "xi.x", "zeta1", "zeta2")){
    .sliced_clamped_raster <- logtransform_layer(.sliced_clamped_raster, match(x, .selected))
}

res <- lapply(.selected, function(x) {
  knitr::knit_child(text = c(
    '### `r x`',
    '',
    '```r',
    'raster_select(.sliced_clamped_raster, band_id = match(x, .selected)) %>% 
  make_leaflet_map(ttl = x, groups = zeta_results$spatial_scales, style = "continuous")',
    '```',
    ''
  ), envir = environment(), quiet = TRUE)
})
cat(unlist(res), sep = '\n')

Distribution plots {.tabset}

ind <- four_values_check(sliced_clamped_raster)
sliced_clamped_raster <- sliced_clamped_raster[ind]
spatial_scales <- zeta_results$spatial_scales[ind]

res <- lapply(.selected, function(x) {
  knitr::knit_child(text = c(
    '### `r x`',
    '',
    '```r',
    'gb <- make_all_plots(sliced_clamped_raster, spatial_scales, match(x, .selected), x, begin = 0.1, end = 0.85, direction = 1, option = "viridis")',
    'layout_mat <- matrix(c(1, 1, 1, 1, 1, 2, 2, 2, 3, 3), byrow = FALSE, nrow = 5, ncol = 2)',
    'gridExtra::grid.arrange(grobs = gb, layout_matrix = layout_mat)',
    '```',
    ''
  ), envir = environment(), quiet = TRUE)
})
cat(unlist(res), sep = '\n')

Comparing $x-$ and $y-$ variables

is_band1_gt_band2 <- function(band1, band2, rasters, .selected, spatial_scales){
    allx <- raster_select(rasters, band_id = match(band1, .selected))$rasters
    ally <- raster_select(rasters, band_id = match(band2, .selected))$rasters
    xy <- foreach(x = allx, y = ally, .combine = rbind) %do% {
        x <- as(x, "Raster")
        y <- as(y, "Raster")
        s <- raster::stack(x, y)
        raster::calc(s, fun = function(x){x[1] > x[2]}) %>% raster::getValues() %>% rstatix::freq_table() %>% dplyr::filter(group == "TRUE") %>% dplyr::select(-group)
    }
    # rownames(xy) <- spatial_scales
  mean_xy <- xy %>% dplyr::summarize(n = mean(n), prop = mean(prop))
  # rownames(mean_xy) <- "mean"
  xy <- rbind(xy, mean_xy)
    xy  %>% 
    knitr::kable("html", caption = paste0("Is ", band1, " greater than ", band2, "?")) %>% 
    kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),
        position = "center",
        full_width = TRUE)
}

vars <- .selected[grepl(".x", .selected)]
res <- foreach(var = vars) %do% {
    is_band1_gt_band2(var, gsub(".x", ".y", var), sliced_clamped_raster, .selected, spatial_scales)
}
cat(unlist(res), sep = '\n')

Cross-scale correlations

crossscale_correlations(
    raster_list = sliced_clamped_raster,
    att_names = band_names,
    selected = .selected,
    spatial_scales = spatial_scales,
    corr_type = "pearson")

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