diff --git a/geohabnet/DESCRIPTION b/geohabnet/DESCRIPTION index ae3c459..15d21fa 100644 --- a/geohabnet/DESCRIPTION +++ b/geohabnet/DESCRIPTION @@ -1,7 +1,7 @@ Package: geohabnet Title: Geographical Risk Analysis Based on Habitat Connectivity -Version: 2.0.0 -Date: 2024-02-27 +Version: 2.1.1 +Date: 2024-04-05 Authors@R: c( person("Krishna", "Keshav", , "kkeshav@ufl.edu", role = c("aut", "cre")), person("Aaron", "Plex", , "plexaaron@ufl.edu", role = "aut", diff --git a/geohabnet/NEWS.md b/geohabnet/NEWS.md index 4988a68..ea552cb 100644 --- a/geohabnet/NEWS.md +++ b/geohabnet/NEWS.md @@ -1,3 +1,10 @@ +# geohabnet 2.1.1 +- Removed dependency with Amazon S3 + +# geohabnet 2.1.0 +- Dispersal kernels can now either be Inverse Power Law or Exponential or both +- Renamed CCRI to HCI + # geohabnet 2.0.0 - Adoption of habitat connectivity - Corrected difference map diff --git a/geohabnet/R/ccri_helper.R b/geohabnet/R/ccri_helper.R deleted file mode 100644 index 2e6f10e..0000000 --- a/geohabnet/R/ccri_helper.R +++ /dev/null @@ -1,416 +0,0 @@ -library(yaml) - -# Global constants -------------------------------------------------------- - -.kparameters_file_type <- "parameters" -.kzeroraster_file_type <- "zero" -.kmapgreybackground_file_type <- "map_grey" - -# utility functions for CCRI ---------------------------------------------- - -.is_packed_rast <- function(x) { - if (tolower(class(x)) == "packedspatraster") { - TRUE - } else { - FALSE - } -} -.unpack_rast_ifnot <- function(x) { - if (.is_packed_rast(x)) { - terra::unwrap(x) - } else { - x - } -} - -.stopifnot_sprast <- function(x) { - stopifnot("Require argument of type SpatRaster" = methods::isClass(x, "SpatRaster")) -} - -.host_map <- function(...) { - if (length(...) == 1) { - return(...[[1]]) - } - # assumes only 2 elements, since only 2 agg methods are supported, - # i.e. sum and mean. Take mean if both are present, else take whatever. - sum(...[[1]], ...[[2]], na.rm = TRUE) / 2 -} - -# Meta-programming approach with eval_tidy -.cal_dist <- function(latilongimatr, method) { - - #---- use Geosphere package, fun distVincentyEllipsoid() is used to calculate the distance, default distance is meter - # reference of standard distance in meter for one degree - - method <- tolower(method) - supported <- dist_methods() - stopifnot("Distance strategy not supported. See dist_methods()\n" = method %in% supported) - - n <- nrow(latilongimatr) - temp_matrix <- matrix(-999, n, n) - - f <- switch(method, - "geodesic" = geosphere::distGeo, - "vincentyellipsoid" = geosphere::distVincentyEllipsoid) - - dvse <- f(c(0, 0), cbind(1, 0)) - - # Calculate the distances - for (i in seq_len(n)) { - temp_matrix[i, ] <- f(round(latilongimatr[i, ], 5), latilongimatr) / dvse - } - - return(temp_matrix) -} - -#' Get risk indices -#' -#' @description Get risk indices from GeoRasters object. -#' @param ri GeoRasters object -#' @return List of risk indices. If the `ri` is global, the list will contain two elements, -#' one for each hemisphere. e.g. `list(east = list(), west = list())`. If the `ri` is not global, -#' the list will contain a single element, e.g. `list()`. -#' @details -#' This function will unpack SpatRasters from GeoMadel and thus is [future::future()] safe. -#' -#' @export -risk_indices <- function(ri) { - stopifnot("Object is not of type GeoRasters" = class(ri) == "GeoRasters") - .ew_split <- function() { - ew_indices <- list(list(), list()) - names(ew_indices) <- c(STR_EAST, STR_WEST) - - for (grast in ri$global_rast) { - for (mod in grast$east) { - ew_indices[[STR_EAST]] <- c(ew_indices[[STR_EAST]], terra::rast(mod$index)) - } - for (mod in grast$west) { - ew_indices[[STR_WEST]] <- c(ew_indices[[STR_WEST]], terra::rast(mod$index)) - } - } - return(ew_indices) - } - - if (ri$global) { - #east-west split - .ew_split() - } else { - unlist(lapply( - ri$rasters, - FUN = function(x) { - terra::rast(x$index) - } - ), recursive = FALSE) - } -} - -.to_ext <- function(geoscale) { - return(terra::ext(geoscale)) -} - -.showmsg <- function(...) { - if (getOption("verbose")) { - message(...) - } -} - -.valid_vector_input <- function(vector_to_check) { - if (!is.vector(vector_to_check) || length(vector_to_check) == 0) { - return(FALSE) - } - return(TRUE) -} - -.utilrast_locs <- function() { - - base_uri <- "https://geohabnet.s3.us-east-2.amazonaws.com/util-rasters/" - full_uri <- function(name) { - paste(base_uri, name, sep = "") - } - - locs <- matrix(c(.kzeroraster_file_type, full_uri("ZeroRaster.tif"), - .kmapgreybackground_file_type, full_uri("map_grey_background.tif"), - "avocado", full_uri("avocado_HarvestedAreaFraction.tif")), # only for testing purpose - ncol = 2, byrow = TRUE) - colnames(locs) <- c("raster", "uri") - return(locs) -} - -.utilrast_uri <- function(typ) { - locs <- .utilrast_locs() - stopifnot("internal error - wrong type" = typ %in% locs[, 1]) - idx <- which(locs[, 1] == typ) - return(locs[idx, ][[2]]) -} - -.download <- function(uri) { - f <- paste(tempfile(), ".tif", sep = "") - stopifnot("download failed " = utils::download.file(uri, - destfile = f, - method = "auto", - mode = "wb", - quiet = getOption("verbose")) == 0) - return(f) -} - -.utilrast <- function(typ) { - return(.download(.utilrast_uri(typ))) -} - -.apply_agg <- function(rast, - reso, - method) { - - return(terra::aggregate(rast, - fact = reso, - fun = method, - na.rm = TRUE, - na.action = stats::na.omit)) -} - -.onLoad <- function(libname, pkgname) { - - # this will cleanup previously created config and replace it with new one - # in case, if it's a first installation, it will simply copy the config. - reset_params() - - .utilrast <<- memoise::memoise(.utilrast) - .cal_mgb <<- memoise::memoise(.cal_mgb) -} - -.get_helper_filepath <- function(file_type) { - file_path <- if (file_type == "parameters") { - system.file("parameters.yaml", - package = "geohabnet", - mustWork = TRUE) - } else { - .utilrast(file_type) - } - return(file_path) -} - -# Recursively check the structure of nested sections -.check_nested_structure <- function(existing_section, provided_section, - section_name) { - if (!identical(names(existing_section), names(provided_section))) { - stop(paste("The", section_name, "section in the provided YAML file does not - match the structure of the existing YAML file.")) - } - - # Check the structure of nested sections - for (key in names(existing_section)) { - if (is.list(existing_section[[key]]) && is.list(provided_section[[key]])) { - .check_nested_structure( - existing_section[[key]], provided_section[[key]], - paste(section_name, key, sep = " > ") - ) - } - } -} - -.check_yaml_structure <- function(existing_yaml_file, provided_yaml_file) { - # Read the existing YAML file - existing_yaml <- yaml::yaml.load_file(existing_yaml_file) - - # Read the provided YAML file - provided_yaml <- yaml::yaml.load_file(provided_yaml_file) - - # Compare the structure of the YAML files - if (!identical(names(existing_yaml), names(provided_yaml))) { - stop("The provided YAML file does not match the structure of the existing - YAML file.") - } - # Check the structure of nested sections - for (key in names(existing_yaml)) { - if (is.list(existing_yaml[[key]]) && is.list(provided_yaml[[key]])) { - .check_nested_structure(existing_yaml[[key]], provided_yaml[[key]], key) - } - } - # If the function reaches this point, the YAML structures match - return(TRUE) -} - -.param_hosts <- function(param_config = load_parameters()) { - return(paste(param_config$`CCRI parameters`$Hosts, collapse = ", ")) -} - -.cal_mgb <- function(geoscale, isglobal) { - # calculate map grey background - map_grey_background <- terra::rast(.get_helper_filepath(.kmapgreybackground_file_type)) - map_grey_background_ext <- if (isglobal == FALSE) { - terra::crop(map_grey_background, .to_ext(geoscale)) - } else { - map_grey_background - } - return(map_grey_background_ext) -} - -.cal_zerorast <- function(in_rast, reso) { - - # Create zero_rast with the same dimensions as in_rast - zero_rast <- terra::rast(.get_helper_filepath(.kzeroraster_file_type)) - # Set extent of zero_rast to match in_rast - zero_rast <- terra::resample(zero_rast, in_rast, threads = TRUE) - - return(zero_rast) -} - - -.get_palette <- function() { - palette1 <- c( - "#F4E156FF", "#F6D746FF", "#F8CD37FF", "#FAC329FF", "#FBB91EFF", - "#FCAF13FF", "#FCA50BFF", "#FB9C06FF", "#FA9207FF", "#F8890CFF", - "#F68013FF", "#F37819FF", "#F06F20FF", "#EC6727FF", "#E85F2EFF", - "#E25834FF", "#DD5139FF", "#D74B3FFF", "#D04545FF", "#CA404AFF", - "#C33B4FFF", "#BC3754FF", "#B43359FF", "#AC305EFF", "#A42C60FF", - "#9B2964FF", "#932667FF", "#922568FF", "#902568FF", "#8F2469FF", - "#8D2369FF", "#8C2369FF", "#8A226AFF", "#88226AFF", "#87216BFF", - "#85216BFF", "#84206BFF", "#82206CFF", "#801F6CFF", "#7F1E6CFF", - "#7D1E6DFF", "#7C1D6DFF", "#7A1D6DFF", "#781C6DFF", "#771C6DFF", - "#751B6EFF", "#741A6EFF", "#721A6EFF", "#71196EFF", "#6E196EFF", - "#6D186EFF", "#6B186EFF", "#6A176EFF", "#68166EFF", "#66166EFF", - "#65156EFF", "#63156EFF", "#61136EFF", "#60136EFF", "#5E126EFF", - "#5C126EFF", "#5B126EFF", "#59106EFF", "#58106EFF", "#560F6DFF", - "#540F6DFF", "#530E6DFF", "#510E6CFF", "#500D6CFF", "#4D0D6CFF", - "#4C0C6BFF", "#4A0C6BFF", "#490B6AFF", "#470B6AFF", "#450A69FF", - "#440A68FF", "#420A68FF", "#400A67FF", "#3E0966FF", "#3D0965FF", - "#3B0964FF", "#390963FF", "#380962FF", "#360961FF", "#340A5FFF", - "#320A5EFF", "#310A5CFF", "#2F0A5BFF", "#2D0B59FF", "#2B0B57FF", - "#290B55FF", "#280B53FF", "#250C51FF", "#240C4EFF", "#230C4BFF", - "#200C49FF", "#1F0C47FF", "#1D0C44FF", "#1C0C42FF", "#1A0C40FF", - "#190C3DFF", "#170C3BFF", "#150B38FF", "#150B36FF", "#130A33FF", - "#110A31FF", "#11092EFF", "#0F092CFF", "#0D082AFF", "#0C0827FF", - "#0B0725FF", "#0A0723FF", "#090620FF", "#08051EFF", "#07051CFF", - "#060419FF", "#050418FF", "#040315FF", "#040312FF", "#030210FF", - "#02020EFF", "#02020CFF", "#02010AFF", "#010108FF", "#010106FF", - "#010005FF", "#000004FF", "#000004FF", "#000004FF" - ) - return(palette1) -} - -.get_palette_for_diffmap <- function() { - - # ```{r ,fig.width=6, fig.height=7, dpi=150} - paldif <- viridisLite::viridis(80, direction = 1, alpha = 0.9) - return(paldif) -} - -.write_yaml <- function(yaml_obj, file_path) { - # Validate YAML object - if (is.null(yaml_obj) || !is.list(yaml_obj)) { - stop("Invalid YAML object. Please provide a non-null list as the YAML object.") - } - - # Validate file path and type - if (!is.character(file_path) || !grepl("\\.yaml$|\\.yml$", file_path, ignore.case = TRUE)) { - stop("Invalid file path. Please provide a valid YAML file path with '.yaml' or '.yml' extension.") - } - - # Write YAML to file - tryCatch( - yaml::write_yaml(yaml_obj, file_path), - error = function(e) { - stop("Error writing YAML to file:", conditionMessage(e)) - } - ) - - .showmsg("YAML object successfully written to file: ", file_path) -} - -.get_cropharvest_raster_helper <- function(crop_name, data_source) { - if (data_source == "monfreda") { - geodata::crop_monfreda(crop = crop_name, path = tempdir(), var = "area_f") - } else if (data_source %in% c("mapspam2010", "mapspam2017Africa")) { - x <- if (data_source == "mapspam2010") { - sp_rast(crp = crop_name) - } else { - sp_rast(crp = crop_name, africa = TRUE) - } - x * 0.0001 - } else { - stop(paste("unsupported source: ", data_source)) - } -} - -#' Get supported sources of crops -#' -#' When provided, [cropharvest_rast()] will -#' look for cropland data in this specific source. -#' @returns Vector of supported sources. -#' Also used as a lookup to find get raster object. -#' @export -#' @examples -#' # Get currently supported sources -#' supported_sources() -supported_sources <- function() { - return(c(monfreda(), mapspam())) -} - -#' Supported sources for monfreda -#' -#' @export -monfreda <- function() { - return(c("monfreda")) -} - -#' Supported sources for Mapspam -#' -#' @export -mapspam <- function() { - return(c("mapspam2010", "mapspam2017Africa")) -} - -#' Search for crop -#' -#' It returns the dataset sources in which crop data is available. -#' It's a wrapper around [geodata::spamCrops()] and [geodata::monfredaCrops()] -#' @param name name of crop -#' @return Logical. Sources iin crop data is available. -#' @export -#' @examples -#' search_crop("coffee") -#' search_crop("wheat") -#' \donttest{ -#' search_crop("jackfruit") -#' } -#' -#' @seealso [supported_sources()] -search_crop <- function(name) { - crp <- tolower(trimws(name)) - - funs <- c("monfreda", "spam") - srcs <- character(0) - - for (src in funs) { - f <- paste0("geodata::", src, "Crops()") - res <- rlang::eval_tidy(rlang::parse_expr(f)) - if (src == "monfreda") { - res <- res$name - } - if (crp %in% res) { - srcs <- c(srcs, src) - } - } - - srcs <- if (is.null(srcs) || length(srcs) < 1) { - "Crop not present in supported sources." - } else { - srcs - } - - return(srcs) -} - -#' Distance methods supported -#' -#' Contains supported strategies to calculate distance between two points. -#' Use of one the methods in [sean()] or [sensitivity_analysis()]. -#' @return vector -#' @export -#' -#' @examples -#' dist_methods() -#' -dist_methods <- function() { - return(c("geodesic", "vincentyellipsoid")) -} diff --git a/geohabnet/R/hci_helper.R b/geohabnet/R/hci_helper.R index b2943d6..bb6378b 100644 --- a/geohabnet/R/hci_helper.R +++ b/geohabnet/R/hci_helper.R @@ -3,8 +3,8 @@ library(yaml) # Global constants -------------------------------------------------------- .kparameters_file_type <- "parameters" -.kzeroraster_file_type <- "zero" -.kmapgreybackground_file_type <- "map_grey" +.kzeroraster_fname <- "ZeroRaster.tif" +.kmapgreybackground_fname <- "map_grey_background.tif" # utility functions for CCRI ---------------------------------------------- @@ -121,28 +121,6 @@ risk_indices <- function(ri) { return(TRUE) } -.utilrast_locs <- function() { - - base_uri <- "https://geohabnet.s3.us-east-2.amazonaws.com/util-rasters/" - full_uri <- function(name) { - paste(base_uri, name, sep = "") - } - - locs <- matrix(c(.kzeroraster_file_type, full_uri("ZeroRaster.tif"), - .kmapgreybackground_file_type, full_uri("map_grey_background.tif"), - "avocado", full_uri("avocado_HarvestedAreaFraction.tif")), # only for testing purpose - ncol = 2, byrow = TRUE) - colnames(locs) <- c("raster", "uri") - return(locs) -} - -.utilrast_uri <- function(typ) { - locs <- .utilrast_locs() - stopifnot("internal error - wrong type" = typ %in% locs[, 1]) - idx <- which(locs[, 1] == typ) - return(locs[idx, ][[2]]) -} - .download <- function(uri) { f <- paste(tempfile(), ".tif", sep = "") stopifnot("download failed " = utils::download.file(uri, @@ -153,8 +131,13 @@ risk_indices <- function(ri) { return(f) } -.utilrast <- function(typ) { - return(.download(.utilrast_uri(typ))) +.utilrast <- function(fname) { + + stopifnot("Internal error. TIFF file not available." = + fname %in% c(.kzeroraster_fname, .kmapgreybackground_fname)) + return(system.file(fname, + package = utils::packageName(), + mustWork = TRUE)) } .apply_agg <- function(rast, @@ -176,7 +159,6 @@ risk_indices <- function(ri) { .utilrast <<- memoise::memoise(.utilrast) .cal_mgb <<- memoise::memoise(.cal_mgb) - #.apply_agg <<- memoise::memoise(.apply_agg) } .get_helper_filepath <- function(file_type) { @@ -237,7 +219,7 @@ risk_indices <- function(ri) { .cal_mgb <- function(geoscale, isglobal) { # calculate map grey background - map_grey_background <- terra::rast(.get_helper_filepath(.kmapgreybackground_file_type)) + map_grey_background <- terra::rast(.get_helper_filepath(.kmapgreybackground_fname)) map_grey_background_ext <- if (isglobal == FALSE) { terra::crop(map_grey_background, .to_ext(geoscale)) } else { @@ -246,17 +228,16 @@ risk_indices <- function(ri) { return(map_grey_background_ext) } -.cal_zerorast <- function(in_rast, reso) { +.cal_zerorast <- function(in_rast) { # Create zero_rast with the same dimensions as in_rast - zero_rast <- terra::rast(.get_helper_filepath(.kzeroraster_file_type)) + zero_rast <- terra::rast(.get_helper_filepath(.kzeroraster_fname)) # Set extent of zero_rast to match in_rast zero_rast <- terra::resample(zero_rast, in_rast, threads = TRUE) return(zero_rast) } - .get_palette <- function() { palette1 <- c( "#F4E156FF", "#F6D746FF", "#F8CD37FF", "#FAC329FF", "#FBB91EFF", diff --git a/geohabnet/R/maps.R b/geohabnet/R/maps.R index 42b9961..a6394c0 100644 --- a/geohabnet/R/maps.R +++ b/geohabnet/R/maps.R @@ -49,7 +49,7 @@ hci_mean <- function(indices, } dis_mean_id <- terra::disagg(mean_index, fact = c(res, res)) - toplot <- dis_mean_id + .cal_zerorast(dis_mean_id, res) + toplot <- dis_mean_id + .cal_zerorast(dis_mean_id) plt_ret <- if (plt == TRUE) { .plot(toplot, @@ -93,7 +93,7 @@ hci_variance <- function(indices, var_disag_rast[var_disag_rast[] == 0] <- NA - var_disag_rast + .cal_zerorast(var_disag_rast, res) + var_disag_rast + .cal_zerorast(var_disag_rast) } geoext <- geoscale @@ -185,7 +185,7 @@ hci_diff <- function(x, diagg_rast <- terra::disagg(scaled_rast, fact = c(res, res)) - diagg_rast + .cal_zerorast(diagg_rast, res) + diagg_rast + .cal_zerorast(diagg_rast) } diff_out <- .cal_diff(x, y, geoext) diff --git a/geohabnet/R/metrics.R b/geohabnet/R/metrics.R index ec653bb..9385d2a 100644 --- a/geohabnet/R/metrics.R +++ b/geohabnet/R/metrics.R @@ -98,7 +98,7 @@ betweeness <- function(crop_dm, ...) { #' @rdname nn_sum ev <- function(crop_dm, ...) { - eigenvectorvalues <- igraph::evcent(crop_dm, ...) + eigenvectorvalues <- igraph::eigen_centrality(crop_dm, ...) evv <- eigenvectorvalues$vector evv[is.na(evv)] <- 0 evp <- if (max(evv) == 0) { diff --git a/geohabnet/geohabnet-report.html b/geohabnet/geohabnet-report.html deleted file mode 100644 index 2c6f62a..0000000 --- a/geohabnet/geohabnet-report.html +++ /dev/null @@ -1,13751 +0,0 @@ - - -
- - - - - - - - - - - - - - - - - - - - - - -1 | -- |
- library(yaml)- |
-
2 | -- | - - | -
3 | -- |
- # Global constants --------------------------------------------------------- |
-
4 | -- | - - | -
5 | -- |
- .kparameters_file_type <- "parameters"- |
-
6 | -- |
- .kzeroraster_file_type <- "zero"- |
-
7 | -- |
- .kmapgreybackground_file_type <- "map_grey"- |
-
8 | -- | - - | -
9 | -- |
- # utility functions for CCRI ----------------------------------------------- |
-
10 | -- |
- # Meta-programming approach with eval_tidy- |
-
11 | -- |
- .cal_dist <- function(latilongimatr, method) {- |
-
12 | -- | - - | -
13 | -8x | -
- method <- tolower(method)- |
-
14 | -8x | -
- supported <- dist_methods()- |
-
15 | -8x | -
- stopifnot("Distance strategy not supported. See dist_methods()\n" = method %in% supported)- |
-
16 | -- | - - | -
17 | -8x | -
- n <- nrow(latilongimatr)- |
-
18 | -8x | -
- temp_matrix <- matrix(-999, n, n)- |
-
19 | -- | - - | -
20 | -8x | -
- f <- switch(method,- |
-
21 | -8x | -
- "geodesic" = geosphere::distGeo,- |
-
22 | -8x | -
- "vincentyellipsoid" = geosphere::distVincentyEllipsoid)- |
-
23 | -- | - - | -
24 | -8x | -
- dvse <- f(c(0, 0), cbind(1, 0))- |
-
25 | -- | - - | -
26 | -- |
- # Calculate the distances- |
-
27 | -8x | -
- for (i in seq_len(n)) {- |
-
28 | -2116x | -
- temp_matrix[i, ] <- f(round(latilongimatr[i, ], 5), latilongimatr) / dvse- |
-
29 | -- |
- }- |
-
30 | -- | - - | -
31 | -8x | -
- return(temp_matrix)- |
-
32 | -- |
- }- |
-
33 | -- | - - | -
34 | -- |
- .flatten_ri <- function(isglobal, ri) {- |
-
35 | -3x | -
- .ew_split <- function() {- |
-
36 | -3x | -
- ew_indices <- list(list(), list())- |
-
37 | -3x | -
- names(ew_indices) <- c(STR_EAST, STR_WEST)- |
-
38 | -3x | -
- for (indices in ri) {- |
-
39 | -4x | -
- ew_indices[[STR_EAST]] <- c(ew_indices[[STR_EAST]], indices[[STR_EAST]])- |
-
40 | -4x | -
- ew_indices[[STR_WEST]] <- c(ew_indices[[STR_WEST]], indices[[STR_WEST]])- |
-
41 | -- |
- }- |
-
42 | -3x | -
- return(ew_indices)- |
-
43 | -- |
- }- |
-
44 | -- | - - | -
45 | -3x | -
- if (isglobal) {- |
-
46 | -- |
- #east-west split- |
-
47 | -3x | -
- .ew_split()- |
-
48 | -- |
- } else {- |
-
49 | -! | -
- unlist(ri, recursive = FALSE)- |
-
50 | -- |
- }- |
-
51 | -- |
- }- |
-
52 | -- | - - | -
53 | -- |
- .to_ext <- function(geoscale) {- |
-
54 | -12x | -
- return(terra::ext(geoscale))- |
-
55 | -- |
- }- |
-
56 | -- | - - | -
57 | -- |
- .valid_vector_input <- function(vector_to_check) {- |
-
58 | -! | -
- if (!is.vector(vector_to_check) || length(vector_to_check) == 0) {- |
-
59 | -! | -
- return(FALSE)- |
-
60 | -- |
- }- |
-
61 | -! | -
- return(TRUE)- |
-
62 | -- |
- }- |
-
63 | -- | - - | -
64 | -- |
- .utilrast_locs <- function() {- |
-
65 | -- | - - | -
66 | -3x | -
- base_uri <- "https://geohabnet.s3.us-east-2.amazonaws.com/util-rasters/"- |
-
67 | -3x | -
- full_uri <- function(name) {- |
-
68 | -9x | -
- paste(base_uri, name, sep = "")- |
-
69 | -- |
- }- |
-
70 | -- | - - | -
71 | -3x | -
- locs <- matrix(c(.kzeroraster_file_type, full_uri("ZeroRaster.tif"),- |
-
72 | -3x | -
- .kmapgreybackground_file_type, full_uri("map_grey_background.tif"),- |
-
73 | -3x | -
- "avocado", full_uri("avocado_HarvestedAreaFraction.tif")), # only for testing purpose- |
-
74 | -3x | -
- ncol = 2, byrow = TRUE)- |
-
75 | -3x | -
- colnames(locs) <- c("raster", "uri")- |
-
76 | -3x | -
- return(locs)- |
-
77 | -- |
- }- |
-
78 | -- | - - | -
79 | -- |
- .utilrast_uri <- function(typ) {- |
-
80 | -3x | -
- locs <- .utilrast_locs()- |
-
81 | -3x | -
- stopifnot("internal error - wrong type" = typ %in% locs[, 1])- |
-
82 | -3x | -
- idx <- which(locs[, 1] == typ)- |
-
83 | -3x | -
- return(locs[idx, ][[2]])- |
-
84 | -- |
- }- |
-
85 | -- | - - | -
86 | -- |
- .download <- function(uri) {- |
-
87 | -3x | -
- f <- paste(tempfile(), ".tif", sep = "")- |
-
88 | -3x | -
- stopifnot("dowload failed " = utils::download.file(uri, destfile = f, method = "auto", mode = "wb") == 0)- |
-
89 | -3x | -
- return(f)- |
-
90 | -- |
- }- |
-
91 | -- | - - | -
92 | -- |
- .utilrast <- function(typ) {- |
-
93 | -3x | -
- return(.download(.utilrast_uri(typ)))- |
-
94 | -- |
- }- |
-
95 | -- | - - | -
96 | -- |
- .apply_agg <- function(rast,- |
-
97 | -- |
- reso,- |
-
98 | -- |
- method) {- |
-
99 | -- | - - | -
100 | -4x | -
- return(terra::aggregate(rast,- |
-
101 | -4x | -
- fact = reso,- |
-
102 | -4x | -
- fun = method,- |
-
103 | -4x | -
- na.action = stats::na.omit))- |
-
104 | -- |
- }- |
-
105 | -- | - - | -
106 | -- |
- .onLoad <- function(libname, pkgname) {- |
-
107 | -! | -
- .utilrast <<- memoise::memoise(.utilrast)- |
-
108 | -! | -
- .cal_mgb <<- memoise::memoise(.cal_mgb)- |
-
109 | -! | -
- .apply_agg <<- memoise::memoise(.apply_agg)- |
-
110 | -- |
- #.crop_rast <<- memoise::memoise(.crop_rast)- |
-
111 | -- | - - | -
112 | -! | -
- metric_funs <<- stats::setNames(metric_funs,- |
-
113 | -! | -
- c(STR_NEAREST_NEIGHBORS_SUM,- |
-
114 | -! | -
- STR_NODE_STRENGTH,- |
-
115 | -! | -
- STR_BETWEENNESS,- |
-
116 | -! | -
- STR_EIGEN_VECTOR_CENTRALITY,- |
-
117 | -! | -
- STR_CLOSENESS_CENTRALITY,- |
-
118 | -! | -
- STR_PAGE_RANK,- |
-
119 | -! | -
- STR_DEGREE))- |
-
120 | -- |
- }- |
-
121 | -- | - - | -
122 | -- |
- .get_helper_filepath <- function(file_type) {- |
-
123 | -12x | -
- file_path <- if (file_type == "parameters") {- |
-
124 | -5x | -
- system.file("parameters.yaml",- |
-
125 | -5x | -
- package = "geohabnet",- |
-
126 | -5x | -
- mustWork = TRUE- |
-
127 | -- |
- )- |
-
128 | -- |
- } else {- |
-
129 | -7x | -
- .utilrast(file_type)- |
-
130 | -- |
- }- |
-
131 | -12x | -
- return(file_path)- |
-
132 | -- |
- }- |
-
133 | -- | - - | -
134 | -- |
- # Recursively check the structure of nested sections- |
-
135 | -- |
- .check_nested_structure <- function(existing_section, provided_section,- |
-
136 | -- |
- section_name) {- |
-
137 | -13x | -
- if (!identical(names(existing_section), names(provided_section))) {- |
-
138 | -! | -
- stop(paste("The", section_name, "section in the provided YAML file does not- |
-
139 | -! | -
- match the structure of the existing YAML file."))- |
-
140 | -- |
- }- |
-
141 | -- | - - | -
142 | -- |
- # Check the structure of nested sections- |
-
143 | -13x | -
- for (key in names(existing_section)) {- |
-
144 | -29x | -
- if (is.list(existing_section[[key]]) && is.list(provided_section[[key]])) {- |
-
145 | -12x | -
- .check_nested_structure(- |
-
146 | -12x | -
- existing_section[[key]], provided_section[[key]],- |
-
147 | -12x | -
- paste(section_name, key, sep = " > ")- |
-
148 | -- |
- )- |
-
149 | -- |
- }- |
-
150 | -- |
- }- |
-
151 | -- |
- }- |
-
152 | -- | - - | -
153 | -- |
- .check_yaml_structure <- function(existing_yaml_file, provided_yaml_file) {- |
-
154 | -- |
- # Read the existing YAML file- |
-
155 | -1x | -
- existing_yaml <- yaml::yaml.load_file(existing_yaml_file)- |
-
156 | -- | - - | -
157 | -- |
- # Read the provided YAML file- |
-
158 | -1x | -
- provided_yaml <- yaml::yaml.load_file(provided_yaml_file)- |
-
159 | -- | - - | -
160 | -- |
- # Compare the structure of the YAML files- |
-
161 | -1x | -
- if (!identical(names(existing_yaml), names(provided_yaml))) {- |
-
162 | -! | -
- stop("The provided YAML file does not match the structure of the existing- |
-
163 | -! | -
- YAML file.")- |
-
164 | -- |
- }- |
-
165 | -- |
- # Check the structure of nested sections- |
-
166 | -1x | -
- for (key in names(existing_yaml)) {- |
-
167 | -1x | -
- if (is.list(existing_yaml[[key]]) && is.list(provided_yaml[[key]])) {- |
-
168 | -1x | -
- .check_nested_structure(existing_yaml[[key]], provided_yaml[[key]], key)- |
-
169 | -- |
- }- |
-
170 | -- |
- }- |
-
171 | -- |
- # If the function reaches this point, the YAML structures match- |
-
172 | -1x | -
- return(TRUE)- |
-
173 | -- |
- }- |
-
174 | -- | - - | -
175 | -- |
- .param_hosts <- function(param_config = load_parameters()) {- |
-
176 | -! | -
- return(paste(param_config$`CCRI parameters`$Hosts, collapse = ", "))- |
-
177 | -- |
- }- |
-
178 | -- | - - | -
179 | -- |
- .cal_mgb <- function(geoscale, isglobal) {- |
-
180 | -- |
- # calculate map grey background- |
-
181 | -! | -
- map_grey_background <- terra::rast(.get_helper_filepath(.kmapgreybackground_file_type))- |
-
182 | -! | -
- map_grey_background_ext <- if (isglobal == FALSE) {- |
-
183 | -! | -
- terra::crop(map_grey_background, .to_ext(geoscale))- |
-
184 | -- |
- } else {- |
-
185 | -! | -
- map_grey_background- |
-
186 | -- |
- }- |
-
187 | -! | -
- return(map_grey_background_ext)- |
-
188 | -- |
- }- |
-
189 | -- | - - | -
190 | -- |
- .cal_zerorast <- function(in_rast, reso) {- |
-
191 | -- | - - | -
192 | -- |
- # Create zero_rast with the same dimensions as in_rast- |
-
193 | -4x | -
- zero_rast <- terra::rast(.get_helper_filepath(.kzeroraster_file_type))- |
-
194 | -- |
- # Set extent of zero_rast to match in_rast- |
-
195 | -4x | -
- zero_rast <- terra::resample(zero_rast, in_rast, threads = TRUE)- |
-
196 | -- | - - | -
197 | -4x | -
- return(zero_rast)- |
-
198 | -- |
- }- |
-
199 | -- | - - | -
200 | -- | - - | -
201 | -- |
- .get_palette <- function() {- |
-
202 | -! | -
- palette1 <- c(- |
-
203 | -! | -
- "#F4E156FF", "#F6D746FF", "#F8CD37FF", "#FAC329FF", "#FBB91EFF",- |
-
204 | -! | -
- "#FCAF13FF", "#FCA50BFF", "#FB9C06FF", "#FA9207FF", "#F8890CFF",- |
-
205 | -! | -
- "#F68013FF", "#F37819FF", "#F06F20FF", "#EC6727FF", "#E85F2EFF",- |
-
206 | -! | -
- "#E25834FF", "#DD5139FF", "#D74B3FFF", "#D04545FF", "#CA404AFF",- |
-
207 | -! | -
- "#C33B4FFF", "#BC3754FF", "#B43359FF", "#AC305EFF", "#A42C60FF",- |
-
208 | -! | -
- "#9B2964FF", "#932667FF", "#922568FF", "#902568FF", "#8F2469FF",- |
-
209 | -! | -
- "#8D2369FF", "#8C2369FF", "#8A226AFF", "#88226AFF", "#87216BFF",- |
-
210 | -! | -
- "#85216BFF", "#84206BFF", "#82206CFF", "#801F6CFF", "#7F1E6CFF",- |
-
211 | -! | -
- "#7D1E6DFF", "#7C1D6DFF", "#7A1D6DFF", "#781C6DFF", "#771C6DFF",- |
-
212 | -! | -
- "#751B6EFF", "#741A6EFF", "#721A6EFF", "#71196EFF", "#6E196EFF",- |
-
213 | -! | -
- "#6D186EFF", "#6B186EFF", "#6A176EFF", "#68166EFF", "#66166EFF",- |
-
214 | -! | -
- "#65156EFF", "#63156EFF", "#61136EFF", "#60136EFF", "#5E126EFF",- |
-
215 | -! | -
- "#5C126EFF", "#5B126EFF", "#59106EFF", "#58106EFF", "#560F6DFF",- |
-
216 | -! | -
- "#540F6DFF", "#530E6DFF", "#510E6CFF", "#500D6CFF", "#4D0D6CFF",- |
-
217 | -! | -
- "#4C0C6BFF", "#4A0C6BFF", "#490B6AFF", "#470B6AFF", "#450A69FF",- |
-
218 | -! | -
- "#440A68FF", "#420A68FF", "#400A67FF", "#3E0966FF", "#3D0965FF",- |
-
219 | -! | -
- "#3B0964FF", "#390963FF", "#380962FF", "#360961FF", "#340A5FFF",- |
-
220 | -! | -
- "#320A5EFF", "#310A5CFF", "#2F0A5BFF", "#2D0B59FF", "#2B0B57FF",- |
-
221 | -! | -
- "#290B55FF", "#280B53FF", "#250C51FF", "#240C4EFF", "#230C4BFF",- |
-
222 | -! | -
- "#200C49FF", "#1F0C47FF", "#1D0C44FF", "#1C0C42FF", "#1A0C40FF",- |
-
223 | -! | -
- "#190C3DFF", "#170C3BFF", "#150B38FF", "#150B36FF", "#130A33FF",- |
-
224 | -! | -
- "#110A31FF", "#11092EFF", "#0F092CFF", "#0D082AFF", "#0C0827FF",- |
-
225 | -! | -
- "#0B0725FF", "#0A0723FF", "#090620FF", "#08051EFF", "#07051CFF",- |
-
226 | -! | -
- "#060419FF", "#050418FF", "#040315FF", "#040312FF", "#030210FF",- |
-
227 | -! | -
- "#02020EFF", "#02020CFF", "#02010AFF", "#010108FF", "#010106FF",- |
-
228 | -! | -
- "#010005FF", "#000004FF", "#000004FF", "#000004FF"- |
-
229 | -- |
- )- |
-
230 | -! | -
- return(palette1)- |
-
231 | -- |
- }- |
-
232 | -- | - - | -
233 | -- |
- .get_palette_for_diffmap <- function() {- |
-
234 | -- | - - | -
235 | -- |
- # ```{r ,fig.width=6, fig.height=7, dpi=150}- |
-
236 | -! | -
- paldif <- viridisLite::viridis(80, direction = 1, alpha = 0.9)- |
-
237 | -! | -
- return(paldif)- |
-
238 | -- |
- }- |
-
239 | -- | - - | -
240 | -- |
- .write_yaml <- function(yaml_obj, file_path) {- |
-
241 | -- |
- # Validate YAML object- |
-
242 | -! | -
- if (is.null(yaml_obj) || !is.list(yaml_obj)) {- |
-
243 | -! | -
- stop("Invalid YAML object. Please provide a non-null list as the YAML object.")- |
-
244 | -- |
- }- |
-
245 | -- | - - | -
246 | -- |
- # Validate file path and type- |
-
247 | -! | -
- if (!is.character(file_path) || !grepl("\\.yaml$|\\.yml$", file_path, ignore.case = TRUE)) {- |
-
248 | -! | -
- stop("Invalid file path. Please provide a valid YAML file path with '.yaml' or '.yml' extension.")- |
-
249 | -- |
- }- |
-
250 | -- | - - | -
251 | -- |
- # Write YAML to file- |
-
252 | -! | -
- tryCatch(- |
-
253 | -! | -
- yaml::write_yaml(yaml_obj, file_path),- |
-
254 | -! | -
- error = function(e) {- |
-
255 | -! | -
- stop("Error writing YAML to file:", conditionMessage(e))- |
-
256 | -- |
- }- |
-
257 | -- |
- )- |
-
258 | -- | - - | -
259 | -! | -
- message("YAML object successfully written to file:", file_path)- |
-
260 | -- |
- }- |
-
261 | -- | - - | -
262 | -- |
- .get_cropharvest_raster_helper <- function(crop_name, data_source) {- |
-
263 | -2x | -
- if (data_source == "monfreda") {- |
-
264 | -2x | -
- geodata::crop_monfreda(crop = crop_name, path = tempdir(), var = "area_f")- |
-
265 | -! | -
- } else if (data_source == "mapspam") {- |
-
266 | -! | -
- sp_rast(crp = crop_name) / 10000- |
-
267 | -- |
- } else {- |
-
268 | -! | -
- stop(paste("Encountered unsupported source: ", data_source))- |
-
269 | -- |
- }- |
-
270 | -- |
- }- |
-
271 | -- | - - | -
272 | -- |
- #' Get supported sources of crops- |
-
273 | -- |
- #'- |
-
274 | -- |
- #' When provided, [get_cropharvest_raster()] will- |
-
275 | -- |
- #' look for cropland data in this specific source.- |
-
276 | -- |
- #' @returns return vector of supported sources.- |
-
277 | -- |
- #' Also used as a lookup to find get raster object.- |
-
278 | -- |
- #' @export- |
-
279 | -- |
- #' @examples- |
-
280 | -- |
- #' # Get currently supported sources- |
-
281 | -- |
- #' get_supported_sources()- |
-
282 | -- |
- get_supported_sources <- function() {- |
-
283 | -4x | -
- return(c("monfreda", "mapspam"))- |
-
284 | -- |
- }- |
-
285 | -- | - - | -
286 | -- |
- #' Search for crop- |
-
287 | -- |
- #'- |
-
288 | -- |
- #' It returns the dataset sources in which crop data is available.- |
-
289 | -- |
- #' It's a wrapper around [geodata::spamCrops()] and [geodata::monfredaCrops()]- |
-
290 | -- |
- #' @param name name of crop- |
-
291 | -- |
- #' @export- |
-
292 | -- |
- #' @examples- |
-
293 | -- |
- #' search_crop("coffee")- |
-
294 | -- |
- #' search_crop("wheat")- |
-
295 | -- |
- #' \dontrun{- |
-
296 | -- |
- #' search_crop("jackfruit")- |
-
297 | -- |
- #' }- |
-
298 | -- |
- #' @seealso [get_supported_sources()]- |
-
299 | -- |
- search_crop <- function(name) {- |
-
300 | -2x | -
- crp <- tolower(trimws(name))- |
-
301 | -- | - - | -
302 | -2x | -
- funs <- c("monfreda", "spam")- |
-
303 | -2x | -
- srcs <- character(0)- |
-
304 | -- | - - | -
305 | -2x | -
- for (src in funs) {- |
-
306 | -4x | -
- f <- paste0("geodata::", src, "Crops()")- |
-
307 | -4x | -
- res <- rlang::eval_tidy(rlang::parse_expr(f))- |
-
308 | -4x | -
- if (src == "monfreda") {- |
-
309 | -2x | -
- res <- res$name- |
-
310 | -- |
- }- |
-
311 | -4x | -
- if (crp %in% res) {- |
-
312 | -4x | -
- srcs <- c(srcs, src)- |
-
313 | -- |
- }- |
-
314 | -- |
- }- |
-
315 | -- | - - | -
316 | -2x | -
- stopifnot("Crop not present in supported sources." = length(srcs) > 0)- |
-
317 | -- | - - | -
318 | -2x | -
- return(srcs)- |
-
319 | -- |
- }- |
-
320 | -- | - - | -
321 | -- |
- #' Distance methods supported- |
-
322 | -- |
- #'- |
-
323 | -- |
- #' Contains supported strategies to calculate distance between two points.- |
-
324 | -- |
- #' Use of one the methods in [sean()] or [sensitivity_analysis()].- |
-
325 | -- |
- #' @return vector- |
-
326 | -- |
- #' @export- |
-
327 | -- |
- dist_methods <- function() {- |
-
328 | -8x | -
- return(c("geodesic", "vincentyellipsoid"))- |
-
329 | -- |
- }- |
-
1 | -- |
- #' raster for mapspam crop.- |
-
2 | -- |
- #'- |
-
3 | -- |
- #' get raster for crop in mapspam dataset- |
-
4 | -- |
- #' @param crp character. name of a crop. Case-insensitive.- |
-
5 | -- |
- #' @return spatRaster- |
-
6 | -- |
- #' @details- |
-
7 | -- |
- #' See [geodata::spamCrops()] for supported crops.- |
-
8 | -- |
- #' @export- |
-
9 | -- |
- #' @examples- |
-
10 | -- |
- #' \dontrun{- |
-
11 | -- |
- #' sp_rast("rice")- |
-
12 | -- |
- #' }- |
-
13 | -- |
- #' @seealso- |
-
14 | -- |
- #' [geodata::spamCrops()]- |
-
15 | -- |
- #' [search_crop()]- |
-
16 | -- |
- #' @references- |
-
17 | -- |
- #' www.mapspam.com/data- |
-
18 | -- |
- #' International Food Policy Research Institute, 2020.- |
-
19 | -- |
- #' Spatially-Disaggregated Crop Production Statistics Data in Africa South of the Sahara for 2017.- |
-
20 | -- |
- #' https://doi.org/10.7910/DVN/FSSKBW, Harvard Dataverse, V2- |
-
21 | -- |
- sp_rast <- function(crp) {- |
-
22 | -! | -
- return(terra::rast(.gen_url(crp)))- |
-
23 | -- |
- }- |
-
24 | -- | - - | -
25 | -- |
- .gen_url <- function(crop) {- |
-
26 | -- | - - | -
27 | -! | -
- crp <- tolower(trimws(crop))- |
-
28 | -! | -
- crops <- geodata::spamCrops()- |
-
29 | -! | -
- if (!(crp %in% crops)) {- |
-
30 | -! | -
- stop("crop not in mapspam; see spamCrops()")- |
-
31 | -- |
- }- |
-
32 | -! | -
- i <- which(crp == crops)[1]- |
-
33 | -! | -
- if (i > nrow(crops))- |
-
34 | -! | -
- i <- i - nrow(crops)- |
-
35 | -! | -
- crp <- toupper(crops[i, 2])- |
-
36 | -! | -
- uri <- paste("https://geohabnet.s3.us-east-2.amazonaws.com/spamcrops/spam2010V2r0_global_H",- |
-
37 | -! | -
- crp,- |
-
38 | -! | -
- "A.tif",- |
-
39 | -! | -
- sep = "_")- |
-
40 | -- | - - | -
41 | -! | -
- return(.download(uri))- |
-
42 | -- |
- }- |
-
1 | -- |
- #' @exportPattern ^[^\\.].*- |
-
2 | -- | - - | -
3 | -- |
- # Utility functions -------------------------------------------------------- |
-
4 | -- | - - | -
5 | -- |
- .loadparam_ifnull <- function() {- |
-
6 | -28x | -
- if (is.null(the$parameters_config)) {- |
-
7 | -! | -
- the$parameters_config <- load_parameters()- |
-
8 | -- |
- }- |
-
9 | -- |
- }- |
-
10 | -- | - - | -
11 | -- |
- #----------- Extract cropland density data ------------------------ |
-
12 | -- |
- .extract_cropland_density <- function(cropharvest_agg_crop, host_density_threshold) {- |
-
13 | -8x | -
- crop_values <- terra::values(cropharvest_agg_crop)- |
-
14 | -8x | -
- max_val <- max(crop_values, na.rm = TRUE)- |
-
15 | -8x | -
- if (max_val <= host_density_threshold) {- |
-
16 | -! | -
- stop(paste("host density threshold: ", host_density_threshold,- |
-
17 | -! | -
- " is greater than the max value: ", max_val, " of aggregate raster"))- |
-
18 | -- |
- }- |
-
19 | -8x | -
- crop_cells_above_threshold <- which(crop_values > host_density_threshold)- |
-
20 | -8x | -
- thresholded_crop_values <- crop_values[crop_cells_above_threshold]- |
-
21 | -8x | -
- return(list(agg_crop = cropharvest_agg_crop,- |
-
22 | -8x | -
- crop_value = thresholded_crop_values,- |
-
23 | -8x | -
- crop_values_at = crop_cells_above_threshold))- |
-
24 | -- |
- }- |
-
25 | -- | - - | -
26 | -- |
- # Initialize --------------------------------------------------- |
-
27 | -- |
- #' Only meant to global variables- |
-
28 | -- |
- #' @keywords internal- |
-
29 | -- |
- the <- new.env(parent = emptyenv())- |
-
30 | -- |
- the$is_initialized <- FALSE- |
-
31 | -- |
- the$parameters_config <- NULL- |
-
32 | -- |
- the$distance_matrix <- NULL- |
-
33 | -- |
- the$cropharvest_aggtm <- NULL- |
-
34 | -- |
- the$cropharvest_agglm_crop <- NULL- |
-
35 | -- |
- the$cropharvest_aggtm_crop <- NULL- |
-
36 | -- |
- the$gan <- list(sum = list("east" = NULL, "west" = NULL),- |
-
37 | -- |
- mean = list("east" = NULL, "west" = NULL))- |
-
38 | -- | - - | -
39 | -- |
- .resetgan <- function() {- |
-
40 | -3x | -
- the$cropharvest_aggtm <- NULL- |
-
41 | -3x | -
- the$cropharvest_agglm_crop <- NULL- |
-
42 | -3x | -
- the$cropharvest_aggtm_crop <- NULL- |
-
43 | -3x | -
- the$gan <- list(sum = list("east" = NULL, "west" = NULL),- |
-
44 | -3x | -
- mean = list("east" = NULL, "west" = NULL))- |
-
45 | -3x | -
- invisible()- |
-
46 | -- |
- }- |
-
47 | -- | - - | -
48 | -- |
- .gan_table <- function(row, col, val) {- |
-
49 | -- |
- #stopifnot("Not a spatRaster" = tolower(class(val)) == "spatraster")- |
-
50 | -8x | -
- the$gan[[row]][[col]] <- val- |
-
51 | -8x | -
- invisible()- |
-
52 | -- |
- }- |
-
53 | -- | - - | -
54 | -- |
- .crop_rast <- function(agg_method, cropharvest_agg, resolution, geo_scale) {- |
-
55 | -8x | -
- postagg_rast <- if (agg_method == "sum") {- |
-
56 | -- | - - | -
57 | -4x | -
- the$cropharvest_aggtm <- cropharvest_agg / resolution / resolution # TOTAL MEAN- |
-
58 | -- |
- # crop cropland area for the given extent- |
-
59 | -4x | -
- the$cropharvest_aggtm_crop <- terra::crop(the$cropharvest_aggtm, .to_ext(geo_scale))- |
-
60 | -4x | -
- the$cropharvest_aggtm_crop- |
-
61 | -8x | -
- } else if (agg_method == "mean") {- |
-
62 | -- | - - | -
63 | -4x | -
- the$cropharvest_agglm <- cropharvest_agg- |
-
64 | -- |
- # crop cropland area for the given extent- |
-
65 | -4x | -
- the$cropharvest_agglm_crop <- terra::crop(the$cropharvest_agglm, .to_ext(geo_scale))- |
-
66 | -4x | -
- the$cropharvest_agglm_crop- |
-
67 | -- |
- } else {- |
-
68 | -! | -
- stop("aggregation strategy is not supported")- |
-
69 | -- |
- }- |
-
70 | -8x | -
- return(postagg_rast)- |
-
71 | -- |
- }- |
-
72 | -- | - - | -
73 | -- |
- .init_cd <- function(cropharvest_raster,- |
-
74 | -- |
- resolution = 12,- |
-
75 | -- |
- geo_scale,- |
-
76 | -- |
- host_density_threshold = 0,- |
-
77 | -- |
- agg_method = "sum",- |
-
78 | -- |
- dist_method = "geodesic") {- |
-
79 | -- | - - | -
80 | -- |
- # aggregation will be cached- |
-
81 | -8x | -
- cropharvest_agg <- .apply_agg(cropharvest_raster,- |
-
82 | -8x | -
- resolution,- |
-
83 | -8x | -
- agg_method)- |
-
84 | -- | - - | -
85 | -8x | -
- temp_rast <- .crop_rast(agg_method,- |
-
86 | -8x | -
- cropharvest_agg,- |
-
87 | -8x | -
- resolution,- |
-
88 | -8x | -
- geo_scale)- |
-
89 | -8x | -
- density_data <- .extract_cropland_density(- |
-
90 | -8x | -
- temp_rast,- |
-
91 | -8x | -
- host_density_threshold)- |
-
92 | -- | - - | -
93 | -8x | -
- if (is.null(density_data) || (!is.list(density_data))) {- |
-
94 | -! | -
- stop("unable to extract density data, longitude/latitude")- |
-
95 | -- |
- }- |
-
96 | -- | - - | -
97 | -- |
- # Prepare arguments elements values for the CCRI functions- |
-
98 | -- | - - | -
99 | -- |
- # save the latitude and longitude as new matrix- |
-
100 | -8x | -
- latilongimatr <- terra::xyFromCell(density_data$agg_crop, cell = density_data$crop_values_at)- |
-
101 | -- |
- #---- use Geosphere package, fun distVincentyEllipsoid() is used to calculate the distance, default distance is meter- |
-
102 | -- |
- # reference of standard distance in meter for one degree- |
-
103 | -- |
- #dvse <- geosphere::distVincentyEllipsoid(c(0, 0), cbind(1, 0))- |
-
104 | -8x | -
- latilongimatr <- as.matrix(latilongimatr)- |
-
105 | -8x | -
- temp_matrix <- .cal_dist(latilongimatr, dist_method)- |
-
106 | -- | - - | -
107 | -8x | -
- the$distance_matrix <- temp_matrix- |
-
108 | -- | - - | -
109 | -8x | -
- the$is_initialized <- TRUE- |
-
110 | -8x | -
- return(density_data)- |
-
111 | -- |
- }- |
-
112 | -- | - - | -
113 | -- |
- # Aggregate ------------------------------------------------------ |
-
114 | -- | - - | -
115 | -- |
- .ccri_powerlaw <- function(betas,- |
-
116 | -- |
- link_threshold = 0,- |
-
117 | -- |
- metrics = the$parameters_config$`CCRI parameters`$NetworkMetrics$InversePowerLaw,- |
-
118 | -- |
- rast,- |
-
119 | -- |
- crop_cells_above_threshold = NULL,- |
-
120 | -- |
- thresholded_crop_values = NULL) {- |
-
121 | -- | - - | -
122 | -8x | -
- if (!.validate_index_cal(betas)) {- |
-
123 | -! | -
- return(0)- |
-
124 | -- |
- }- |
-
125 | -- | - - | -
126 | -8x | -
- .loadparam_ifnull()- |
-
127 | -- | - - | -
128 | -8x | -
- index_list <- lapply(betas, model_powerlaw,- |
-
129 | -8x | -
- link_threshold = link_threshold,- |
-
130 | -8x | -
- the$distance_matrix,- |
-
131 | -8x | -
- thresholded_crop_values,- |
-
132 | -8x | -
- adj_mat = NULL,- |
-
133 | -8x | -
- rast,- |
-
134 | -8x | -
- crop_cells_above_threshold,- |
-
135 | -8x | -
- metrics = metrics- |
-
136 | -- |
- )- |
-
137 | -8x | -
- return(index_list)- |
-
138 | -- | - - | -
139 | -- |
- #the$result_index_list <- c(the$result_index_list, index_list)- |
-
140 | -- |
- #invisible()- |
-
141 | -- |
- }- |
-
142 | -- | - - | -
143 | -- |
- .ccri_negative_exp <- function(gammas,- |
-
144 | -- |
- link_threshold = 0,- |
-
145 | -- |
- metrics = the$parameters_config$`CCRI parameters`$NetworkMetrics$InversePowerLaw,- |
-
146 | -- |
- rast,- |
-
147 | -- |
- crop_cells_above_threshold = NULL,- |
-
148 | -- |
- thresholded_crop_values = NULL) {- |
-
149 | -- | - - | -
150 | -8x | -
- if (!.validate_index_cal(gammas)) {- |
-
151 | -! | -
- return(0)- |
-
152 | -- |
- }- |
-
153 | -- | - - | -
154 | -8x | -
- .loadparam_ifnull()- |
-
155 | -- | - - | -
156 | -8x | -
- index_list <- lapply(gammas,- |
-
157 | -8x | -
- model_neg_exp,- |
-
158 | -8x | -
- link_threshold = link_threshold,- |
-
159 | -8x | -
- the$distance_matrix,- |
-
160 | -8x | -
- thresholded_crop_values,- |
-
161 | -8x | -
- adj_mat = NULL,- |
-
162 | -8x | -
- rast,- |
-
163 | -8x | -
- crop_cells_above_threshold,- |
-
164 | -8x | -
- metrics = metrics- |
-
165 | -- |
- )- |
-
166 | -- | - - | -
167 | -8x | -
- return(index_list)- |
-
168 | -- |
- #the$result_index_list <- c(the$result_index_list, index_list)- |
-
169 | -- |
- #invisible()- |
-
170 | -- |
- }- |
-
171 | -- | - - | -
172 | -- | - - | -
173 | -- |
- # Utility functions -------------------------------------------------------- |
-
174 | -- | - - | -
175 | -- |
- .validate_index_cal <- function(vals) {- |
-
176 | -16x | -
- ready <- TRUE- |
-
177 | -16x | -
- if (!the$is_initialized) {- |
-
178 | -! | -
- stop("Not initialized. Call initializeCroplandData()")- |
-
179 | -- |
- }- |
-
180 | -16x | -
- stopifnot("dispersal values missing" = length(vals) > 0)- |
-
181 | -16x | -
- if (!is.vector(vals)) {- |
-
182 | -! | -
- warning("argument is not a vector")- |
-
183 | -! | -
- ready <- FALSE- |
-
184 | -- |
- }- |
-
185 | -16x | -
- return(ready)- |
-
186 | -- |
- }- |
-
187 | -- | - - | -
188 | -- |
- # CCRI functions ----------------------------------------------------------- |
-
189 | -- | - - | -
190 | -- |
- .ccri <- function(- |
-
191 | -- |
- link_threshold = 0,- |
-
192 | -- |
- power_law_metrics = the$parameters_config$`CCRI parameters`$NetworkMetrics$InversePowerLaw,- |
-
193 | -- |
- negative_exponential_metrics = the$parameters_config$`CCRI parameters`$NetworkMetrics$NegativeExponential,- |
-
194 | -- |
- rast,- |
-
195 | -- |
- crop_cells_above_threshold,- |
-
196 | -- |
- thresholded_crop_values) {- |
-
197 | -- | - - | -
198 | -8x | -
- .loadparam_ifnull()- |
-
199 | -8x | -
- risk_indexes <- list()- |
-
200 | -- | - - | -
201 | -- |
- # TODO: parallelize them- |
-
202 | -8x | -
- betas <- as.numeric(the$parameters_config$`CCRI parameters`$DispersalKernelModels$InversePowerLaw$beta)- |
-
203 | -- | - - | -
204 | -8x | -
- if (length(betas) > 0) {- |
-
205 | -8x | -
- stopifnot("beta values are not valid" = is.numeric(betas) == TRUE, is.vector(betas) == TRUE)- |
-
206 | -8x | -
- risk_indexes <- c(risk_indexes,- |
-
207 | -8x | -
- .ccri_powerlaw(betas,- |
-
208 | -8x | -
- link_threshold,- |
-
209 | -8x | -
- metrics = power_law_metrics,- |
-
210 | -8x | -
- rast,- |
-
211 | -8x | -
- crop_cells_above_threshold = crop_cells_above_threshold,- |
-
212 | -8x | -
- thresholded_crop_values = thresholded_crop_values- |
-
213 | -- |
- ))- |
-
214 | -- |
- }- |
-
215 | -- | - - | -
216 | -8x | -
- gammas <- as.numeric(the$parameters_config$`CCRI parameters`$DispersalKernelModels$NegativeExponential$gamma)- |
-
217 | -8x | -
- if (length(gammas) > 0) {- |
-
218 | -8x | -
- stopifnot("gamma values are not valid" = is.numeric(gammas) == TRUE, is.vector(gammas) == TRUE)- |
-
219 | -8x | -
- risk_indexes <- c(risk_indexes,- |
-
220 | -8x | -
- .ccri_negative_exp(gammas,- |
-
221 | -8x | -
- link_threshold,- |
-
222 | -8x | -
- metrics = negative_exponential_metrics,- |
-
223 | -8x | -
- rast,- |
-
224 | -8x | -
- crop_cells_above_threshold = crop_cells_above_threshold,- |
-
225 | -8x | -
- thresholded_crop_values = thresholded_crop_values- |
-
226 | -- |
- ))- |
-
227 | -- |
- }- |
-
228 | -- | - - | -
229 | -8x | -
- return(risk_indexes)- |
-
230 | -- |
- }- |
-
231 | -- | - - | -
232 | -- |
- # Sensitivity analysis ----------------------------------------------------- |
-
233 | -- | - - | -
234 | -- |
- #' Calculate sensitivity analysis on cropland harvested area fraction- |
-
235 | -- |
- #'- |
-
236 | -- |
- #' This function calculates sensitivity analysis on cropland harvested area fraction based on provided parameters.- |
-
237 | -- |
- #' It can be used as an entry point for sensitivity analysis.- |
-
238 | -- |
- #' @param link_threshold numeric. A threshold value for link- |
-
239 | -- |
- #' @param host_density_threshold A host density threshold value- |
-
240 | -- |
- #' @inheritParams sa_onrasters- |
-
241 | -- |
- #' @return A list of calculated CCRI values using negative exponential- |
-
242 | -- |
- #' @export- |
-
243 | -- |
- #' @details- |
-
244 | -- |
- #' When `global = TRUE`, `geoscale` is ignored and [global_scales()] is used- |
-
245 | -- |
- #'- |
-
246 | -- |
- #' @seealso Uses [connectivity()]- |
-
247 | -- |
- sean <- function(rast,- |
-
248 | -- |
- global = TRUE,- |
-
249 | -- |
- geoscale,- |
-
250 | -- |
- agg_methods = c("sum", "mean"),- |
-
251 | -- |
- dist_method = "geodesic",- |
-
252 | -- |
- link_threshold = 0,- |
-
253 | -- |
- host_density_threshold = 0,- |
-
254 | -- |
- res = reso(),- |
-
255 | -- |
- maps = TRUE) {- |
-
256 | -- | - - | -
257 | -2x | -
- .resetgan()- |
-
258 | -2x | -
- .loadparam_ifnull()- |
-
259 | -- | - - | -
260 | -2x | -
- mets <- get_param_metrics(the$parameters_config)- |
-
261 | -- | - - | -
262 | -2x | -
- sean_geo <- function(geoext) {- |
-
263 | -4x | -
- message(- |
-
264 | -4x | -
- paste(- |
-
265 | -4x | -
- "\nRunning sensitivity analysis for the extent: [",- |
-
266 | -4x | -
- paste(geoext, collapse = ", "),- |
-
267 | -4x | -
- "],\n",- |
-
268 | -4x | -
- "Link threshold: ",- |
-
269 | -4x | -
- link_threshold,- |
-
270 | -4x | -
- "\n",- |
-
271 | -4x | -
- "Host density threshold: ",- |
-
272 | -4x | -
- host_density_threshold,- |
-
273 | -4x | -
- "\n"- |
-
274 | -- |
- )- |
-
275 | -- |
- )- |
-
276 | -- | - - | -
277 | -4x | -
- lrisk_indexes <- list()- |
-
278 | -- | - - | -
279 | -4x | -
- for (agg_method in agg_methods) {- |
-
280 | -8x | -
- density_data <- .init_cd(rast,- |
-
281 | -8x | -
- res,- |
-
282 | -8x | -
- geoext,- |
-
283 | -8x | -
- host_density_threshold = host_density_threshold,- |
-
284 | -8x | -
- agg_method,- |
-
285 | -8x | -
- dist_method)- |
-
286 | -- | - - | -
287 | -8x | -
- lrisk_indexes <- c(lrisk_indexes,- |
-
288 | -8x | -
- .ccri(link_threshold,- |
-
289 | -8x | -
- power_law_metrics = mets$pl,- |
-
290 | -8x | -
- negative_exponential_metrics = mets$ne,- |
-
291 | -8x | -
- rast = density_data$agg_crop,- |
-
292 | -8x | -
- crop_cells_above_threshold =- |
-
293 | -8x | -
- density_data$crop_values_at,- |
-
294 | -8x | -
- thresholded_crop_values =- |
-
295 | -8x | -
- density_data$crop_value- |
-
296 | -- |
- ))- |
-
297 | -- |
- }- |
-
298 | -4x | -
- return(lrisk_indexes)- |
-
299 | -- |
- }- |
-
300 | -- | - - | -
301 | -2x | -
- .addto_tab <- function(hemi) {- |
-
302 | -4x | -
- .gan_table("sum", hemi, the$cropharvest_aggtm_crop)- |
-
303 | -4x | -
- .gan_table("mean", hemi, the$cropharvest_agglm_crop)- |
-
304 | -4x | -
- invisible()- |
-
305 | -- |
- }- |
-
306 | -- | - - | -
307 | -2x | -
- risk_indexes <- if (global) {- |
-
308 | -- | - - | -
309 | -2x | -
- global_exts <- global_scales()- |
-
310 | -- | - - | -
311 | -2x | -
- east_indexes <- sean_geo(global_exts[[STR_EAST]])- |
-
312 | -2x | -
- .addto_tab(STR_EAST)- |
-
313 | -- | - - | -
314 | -2x | -
- west_indexes <- sean_geo(global_exts[[STR_WEST]])- |
-
315 | -2x | -
- .addto_tab(STR_WEST)- |
-
316 | -2x | -
- list(east = east_indexes, west = west_indexes)- |
-
317 | -- |
- } else {- |
-
318 | -! | -
- sean_geo(geoscale)- |
-
319 | -- |
- }- |
-
320 | -- | - - | -
321 | -2x | -
- if (maps == TRUE) {- |
-
322 | -! | -
- connectivity(risk_indexes,- |
-
323 | -! | -
- global,- |
-
324 | -! | -
- geoscale,- |
-
325 | -! | -
- res,- |
-
326 | -! | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$MeanCC),- |
-
327 | -! | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$Variance),- |
-
328 | -! | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$Difference)- |
-
329 | -- |
- )- |
-
330 | -- |
- }- |
-
331 | -2x | -
- the$is_initialized <- FALSE- |
-
332 | -2x | -
- return(risk_indexes)- |
-
333 | -- |
- }- |
-
334 | -- | - - | -
335 | -- |
- .sean_linkweights <- function(link_threshold = 0,- |
-
336 | -- |
- host_density_thresholds,- |
-
337 | -- |
- global = TRUE,- |
-
338 | -- |
- geoscale,- |
-
339 | -- |
- agg_methods,- |
-
340 | -- |
- rast,- |
-
341 | -- |
- res,- |
-
342 | -- |
- dist_method = "geodesic",- |
-
343 | -- |
- maps = TRUE) {- |
-
344 | -- | - - | -
345 | -1x | -
- risk_indexes <- lapply(host_density_thresholds,- |
-
346 | -1x | -
- function(threshold) {- |
-
347 | -2x | -
- invisible(- |
-
348 | -2x | -
- sean(- |
-
349 | -2x | -
- link_threshold = link_threshold,- |
-
350 | -2x | -
- host_density_threshold = threshold,- |
-
351 | -2x | -
- global = global,- |
-
352 | -2x | -
- geoscale = geoscale,- |
-
353 | -2x | -
- agg_methods = agg_methods,- |
-
354 | -2x | -
- dist_method = dist_method,- |
-
355 | -2x | -
- rast = rast,- |
-
356 | -2x | -
- res = res,- |
-
357 | -2x | -
- maps = FALSE- |
-
358 | -- |
- )- |
-
359 | -- |
- )- |
-
360 | -- |
- })- |
-
361 | -- | - - | -
362 | -1x | -
- risk_indexes <- .flatten_ri(global, risk_indexes)- |
-
363 | -1x | -
- if (maps == TRUE) {- |
-
364 | -- | - - | -
365 | -! | -
- connectivity(risk_indexes,- |
-
366 | -! | -
- global = global,- |
-
367 | -! | -
- geoscale,- |
-
368 | -! | -
- res,- |
-
369 | -! | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$MeanCC),- |
-
370 | -! | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$Variance),- |
-
371 | -! | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$Difference)- |
-
372 | -- |
- )- |
-
373 | -- |
- }- |
-
374 | -1x | -
- return(risk_indexes)- |
-
375 | -- |
- }- |
-
376 | -- | - - | -
377 | -- |
- #' Run senstivity analysis- |
-
378 | -- |
- #'- |
-
379 | -- |
- #' Same as [sensitivity_analysis()] but it takes raster object and other parameters as an input.- |
-
380 | -- |
- #' @param rast Raster object which will be used in analysis.- |
-
381 | -- |
- #' @seealso Use [get_rasters()] to obtain raster object.- |
-
382 | -- |
- #' @param global Logical. `TRUE` if global analysis, `FALSE` otherwise.- |
-
383 | -- |
- #' Default is `TRUE`- |
-
384 | -- |
- #' @param geoscale Vector. Geographical coordinates- |
-
385 | -- |
- #' in the form of c(Xmin, Xmax, Ymin, Ymax)- |
-
386 | -- |
- #' @param link_thresholds vector. link threshold values- |
-
387 | -- |
- #' @param host_density_thresholds vector. host density threshold values- |
-
388 | -- |
- #' @param agg_methods vector. Aggregation methods- |
-
389 | -- |
- #' @param dist_method character. One of the values from [dist_methods()]- |
-
390 | -- |
- #' @param res numeric.- |
-
391 | -- |
- #' resolution at which operations will run.- |
-
392 | -- |
- #' Default is [reso()]- |
-
393 | -- |
- #' @param maps logical. `TRUE` if maps are to be plotted, `FALSE` otherwise- |
-
394 | -- |
- #' @return A list of calculated CCRI indices after operations.- |
-
395 | -- |
- #' An index is generated for each combination of paramters.- |
-
396 | -- |
- #' One combination is equivalent to [sean()] function.- |
-
397 | -- |
- #' @export- |
-
398 | -- |
- #' @details- |
-
399 | -- |
- #' When `global = TRUE`, `geo_scale` is ignored.- |
-
400 | -- |
- #' Instead uses scales from [global_scales()].- |
-
401 | -- |
- #'- |
-
402 | -- |
- #' @examples- |
-
403 | -- |
- #' \dontrun{- |
-
404 | -- |
- #' rr <- get_rasters(list(monfreda = c("avocado")))- |
-
405 | -- |
- #' sa_onrasters(rr[[1]],- |
-
406 | -- |
- #' global = FALSE,- |
-
407 | -- |
- #' geoscale = c(-115, -75, 5, 32),- |
-
408 | -- |
- #' c(0.0001, 0.00004),- |
-
409 | -- |
- #' c(0.0001, 0.00005),- |
-
410 | -- |
- #' c("sum", "mean"),- |
-
411 | -- |
- #' res = 24)- |
-
412 | -- |
- #' sa_onrasters(rr[[1]],- |
-
413 | -- |
- #' global = TRUE,- |
-
414 | -- |
- #' link_thresholds = c(0.000001),- |
-
415 | -- |
- #' host_density_thresholds = c(0.00015),- |
-
416 | -- |
- #' agg_methods = c("sum"),- |
-
417 | -- |
- #' res = 24)- |
-
418 | -- |
- #'}- |
-
419 | -- |
- #' @inherit sensitivity_analysis seealso- |
-
420 | -- |
- sa_onrasters <- function(rast,- |
-
421 | -- |
- global = TRUE,- |
-
422 | -- |
- geoscale,- |
-
423 | -- |
- link_thresholds,- |
-
424 | -- |
- host_density_thresholds,- |
-
425 | -- |
- agg_methods = c("sum", "mean"),- |
-
426 | -- |
- dist_method = "geodesic",- |
-
427 | -- |
- res = reso(),- |
-
428 | -- |
- maps = TRUE) {- |
-
429 | -- | - - | -
430 | -1x | -
- cat("New analysis started for given raster")- |
-
431 | -- | - - | -
432 | -1x | -
- .loadparam_ifnull()- |
-
433 | -- | - - | -
434 | -1x | -
- risk_indexes <- lapply(link_thresholds,- |
-
435 | -1x | -
- function(lthreshold) {- |
-
436 | -1x | -
- invisible(- |
-
437 | -1x | -
- .sean_linkweights(- |
-
438 | -1x | -
- link_threshold = lthreshold,- |
-
439 | -1x | -
- host_density_thresholds = host_density_thresholds,- |
-
440 | -1x | -
- global = global,- |
-
441 | -1x | -
- geoscale = geoscale,- |
-
442 | -1x | -
- agg_methods = agg_methods,- |
-
443 | -1x | -
- dist_method = dist_method,- |
-
444 | -1x | -
- rast = rast,- |
-
445 | -1x | -
- res = res,- |
-
446 | -1x | -
- maps = FALSE- |
-
447 | -- |
- )- |
-
448 | -- |
- )- |
-
449 | -- |
- })- |
-
450 | -- | - - | -
451 | -- | - - | -
452 | -1x | -
- risk_indexes <- .flatten_ri(global, risk_indexes)- |
-
453 | -- | - - | -
454 | -1x | -
- if (maps == TRUE) {- |
-
455 | -! | -
- connectivity(risk_indexes,- |
-
456 | -! | -
- global,- |
-
457 | -! | -
- geoscale,- |
-
458 | -! | -
- res,- |
-
459 | -! | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$MeanCC),- |
-
460 | -! | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$Variance),- |
-
461 | -! | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$Difference)- |
-
462 | -- |
- )- |
-
463 | -- |
- }- |
-
464 | -1x | -
- return(risk_indexes)- |
-
465 | -- |
- }- |
-
466 | -- | - - | -
467 | -- |
- #' @title Calculate sensitivity analysis on parameters- |
-
468 | -- |
- #'- |
-
469 | -- |
- #' @description- |
-
470 | -- |
- #' This function runs sensitivity analysis on parameters based on- |
-
471 | -- |
- #' parameters provided through [set_parameters()].- |
-
472 | -- |
- #' It can be used as an entry point for CCRI.- |
-
473 | -- |
- #' By default, it runs analysis on global sclaes[global_scales()].- |
-
474 | -- |
- #' After analysis is complete,- |
-
475 | -- |
- #' it will suppress maps for outcomes if `maps = FALSE` or- |
-
476 | -- |
- #' [interactive()] is `FALSE`.- |
-
477 | -- |
- #' @param maps logical. `TRUE` if maps are to be plotted, `FALSE` otherwise- |
-
478 | -- |
- #' @param alert logical. `TRUE` if beep sound is to be played, `FALSE` otherwise- |
-
479 | -- |
- #' @return logical. `TRUE` if analysis is completed, `FALSE` otherwise.- |
-
480 | -- |
- #' Errors are not handled.- |
-
481 | -- |
- #' @export- |
-
482 | -- |
- #' @examples- |
-
483 | -- |
- #' \dontrun{- |
-
484 | -- |
- #' # Run analysis on specified parameters.yaml- |
-
485 | -- |
- #' sensitivity_analysis()- |
-
486 | -- |
- #' sensitivity_analysis(FALSE, FALSE)- |
-
487 | -- |
- #' sensitivity_analysis(TRUE, FALSE)- |
-
488 | -- |
- #' }- |
-
489 | -- |
- #' @details- |
-
490 | -- |
- #'- |
-
491 | -- |
- #' \code{vignette("global_analysis", package = "geohabnet")}- |
-
492 | -- |
- #'- |
-
493 | -- |
- #' @seealso- |
-
494 | -- |
- #' [sa_onrasters()]- |
-
495 | -- |
- #' [sean()]- |
-
496 | -- |
- #' [global_scales()]- |
-
497 | -- |
- #' [get_parameters()]- |
-
498 | -- |
- #' [set_parameters()]- |
-
499 | -- |
- #' [connectivity()]- |
-
500 | -- |
- sensitivity_analysis <- function(maps = TRUE, alert = TRUE) {- |
-
501 | -- | - - | -
502 | -- |
- #.resetglobals()- |
-
503 | -1x | -
- .resetgan()- |
-
504 | -1x | -
- the$is_initialized <- FALSE- |
-
505 | -1x | -
- the$parameters_config <- load_parameters()- |
-
506 | -- | - - | -
507 | -- |
- # cutoff adjacency matrix- |
-
508 | -1x | -
- cropland_thresholds <- the$parameters_config$`CCRI parameters`$HostDensityThreshold- |
-
509 | -- | - - | -
510 | -- |
- # crop data- |
-
511 | -1x | -
- crop_rasters <- get_rasters(the$parameters_config$`CCRI parameters`$Hosts)- |
-
512 | -1x | -
- agg_methods <- the$parameters_config$`CCRI parameters`$AggregationStrategy # list- |
-
513 | -- | - - | -
514 | -- |
- # resolution- |
-
515 | -1x | -
- resolution <- the$parameters_config$`CCRI parameters`$Resolution- |
-
516 | -- | - - | -
517 | -- |
- # global analysis- |
-
518 | -1x | -
- isglobal <- the$parameters_config$`CCRI parameters`$GeoExtent$global- |
-
519 | -1x | -
- geoscale <- geoscale_param()- |
-
520 | -- | - - | -
521 | -1x | -
- risk_indexes <- lapply(crop_rasters,- |
-
522 | -1x | -
- invisible(function(rast) {- |
-
523 | -1x | -
- sa_onrasters(- |
-
524 | -1x | -
- rast = rast,- |
-
525 | -1x | -
- global = isglobal,- |
-
526 | -1x | -
- geoscale = geoscale,- |
-
527 | -1x | -
- link_thresholds = the$parameters_config$`CCRI parameters`$LinkThreshold,- |
-
528 | -1x | -
- host_density_thresholds = cropland_thresholds,- |
-
529 | -1x | -
- agg_methods = agg_methods,- |
-
530 | -1x | -
- dist_method = the$parameters_config$`CCRI parameters`$DistanceStrategy,- |
-
531 | -1x | -
- res = resolution,- |
-
532 | -1x | -
- maps = FALSE- |
-
533 | -- |
- )- |
-
534 | -- |
- }))- |
-
535 | -- | - - | -
536 | -1x | -
- risk_indexes <- .flatten_ri(isglobal, risk_indexes)- |
-
537 | -- | - - | -
538 | -1x | -
- if (maps == TRUE) {- |
-
539 | -1x | -
- connectivity(risk_indexes,- |
-
540 | -1x | -
- isglobal,- |
-
541 | -1x | -
- geoscale,- |
-
542 | -1x | -
- resolution,- |
-
543 | -1x | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$MeanCC),- |
-
544 | -1x | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$Variance),- |
-
545 | -1x | -
- as.logical(the$parameters_config$`CCRI parameters`$PriorityMaps$Difference)- |
-
546 | -- |
- )- |
-
547 | -- |
- }- |
-
548 | -- | - - | -
549 | -1x | -
- message("sensitivity analysis completed. Refer to maps for results.")- |
-
550 | -1x | -
- if (alert == TRUE) {- |
-
551 | -1x | -
- beepr::beep(2)- |
-
552 | -- |
- }- |
-
553 | -- | - - | -
554 | -1x | -
- return(TRUE)- |
-
555 | -- |
- }- |
-
1 | -- |
- #' Calculate and plot maps- |
-
2 | -- |
- #'- |
-
3 | -- |
- #' Calculate mean, variance and difference. The result is produced in form of maps plotted with predefined settings.- |
-
4 | -- |
- #' Currently, the settings for plot cannot be customized.- |
-
5 | -- |
- #' Default value is `TRUE` for all logical arguments- |
-
6 | -- |
- #' @param indexes list of rasters. See details.- |
-
7 | -- |
- #' @param global logical. `TRUE` if global analysis is required, `FALSE` otherwise.- |
-
8 | -- |
- #' When `TRUE`, `geoscale` is ignored. Default is `TRUE`.- |
-
9 | -- |
- #' @param geoscale vector. geographical scale- |
-
10 | -- |
- #' @param res numeric. map resolution- |
-
11 | -- |
- #' @param pmean `TRUE` if map of mean should be plotted, `FALSE` otherwise.- |
-
12 | -- |
- #' @param pvar `TRUE` if variance map should be plotted, `FALSE` otherwise- |
-
13 | -- |
- #' @param pdiff `TRUE` if difference map should be plotted, `FALSE` otherwise- |
-
14 | -- |
- #' @details- |
-
15 | -- |
- #' `indexes` are actually risk resulting from operations on crop's raster and- |
-
16 | -- |
- #' parameters provided in either `parameters.yaml` or [sean()].- |
-
17 | -- |
- #'- |
-
18 | -- |
- #' It will save all the opted plots using - `pmean`, `pvar` and `pdiff`.- |
-
19 | -- |
- #' File will be saved in [getwd()].If [interactive()] is `TRUE`,- |
-
20 | -- |
- #' then plots can be seen in active plot window. E.g. Rstudio- |
-
21 | -- |
- #'- |
-
22 | -- |
- #' @export- |
-
23 | -- |
- connectivity <- function(indexes,- |
-
24 | -- |
- global = TRUE,- |
-
25 | -- |
- geoscale,- |
-
26 | -- |
- res = reso(),- |
-
27 | -- |
- pmean = TRUE,- |
-
28 | -- |
- pvar = TRUE,- |
-
29 | -- |
- pdiff = TRUE) {- |
-
30 | -- | - - | -
31 | -1x | -
- mean_rast <- ccri_mean(indexes, global, geoscale, pmean)- |
-
32 | -- | - - | -
33 | -1x | -
- if (pvar == TRUE) {- |
-
34 | -1x | -
- ccri_variance(- |
-
35 | -1x | -
- indexes,- |
-
36 | -1x | -
- mean_rast,- |
-
37 | -1x | -
- global,- |
-
38 | -1x | -
- geoscale,- |
-
39 | -1x | -
- res)- |
-
40 | -- |
- }- |
-
41 | -- | - - | -
42 | -1x | -
- if (pdiff == TRUE) {- |
-
43 | -1x | -
- if (global) {- |
-
44 | -1x | -
- geoscale <- .global_ext()- |
-
45 | -- |
- }- |
-
46 | -- | - - | -
47 | -1x | -
- ccri_diff(- |
-
48 | -1x | -
- mean_rast,- |
-
49 | -1x | -
- the$cropharvest_aggtm_crop,- |
-
50 | -1x | -
- the$cropharvest_agglm_crop,- |
-
51 | -1x | -
- global,- |
-
52 | -1x | -
- geoscale,- |
-
53 | -1x | -
- res)- |
-
54 | -- |
- }- |
-
55 | -1x | -
- invisible()- |
-
56 | -- |
- }- |
-
57 | -- | - - | -
58 | -- |
- # maps --------------------------------------------------------------------- |
-
59 | -- | - - | -
60 | -- |
- #' Calculate mean of raster objects- |
-
61 | -- |
- #'- |
-
62 | -- |
- #' Wrapper for [terra::mean()]. Calculates mean of list of rasters.- |
-
63 | -- |
- #' @inheritParams connectivity- |
-
64 | -- |
- #' @param plt `TRUE` if need to plot mean map, `FALSE` otherwise and `geoscale`.- |
-
65 | -- |
- #' @export- |
-
66 | -- |
- ccri_mean <- function(indexes,- |
-
67 | -- |
- global = TRUE,- |
-
68 | -- |
- geoscale = NULL,- |
-
69 | -- |
- plt = TRUE) {- |
-
70 | -- | - - | -
71 | -1x | -
- .cal_mean <- function(ext_indices) {- |
-
72 | -2x | -
- mean_idx <- terra::app(terra::rast(ext_indices), sum, na.rm = TRUE) / length(ext_indices)- |
-
73 | -2x | -
- mean_index_vals <- terra::values(mean_idx)- |
-
74 | -2x | -
- zeroid <- which(mean_index_vals == 0)- |
-
75 | -2x | -
- mean_idx[zeroid] <- NaN- |
-
76 | -2x | -
- mean_idx- |
-
77 | -- |
- }- |
-
78 | -- | - - | -
79 | -1x | -
- mean_index <- if (global == TRUE) {- |
-
80 | -- | - - | -
81 | -1x | -
- .gan_paramok(indexes)- |
-
82 | -1x | -
- east <- .cal_mean(indexes[[STR_EAST]])- |
-
83 | -1x | -
- west <- .cal_mean(indexes[[STR_WEST]])- |
-
84 | -1x | -
- geoscale <- .global_ext()- |
-
85 | -1x | -
- terra::merge(east, west)- |
-
86 | -- | - - | -
87 | -- |
- } else {- |
-
88 | -! | -
- .cal_mean(indexes)- |
-
89 | -- |
- }- |
-
90 | -- | - - | -
91 | -1x | -
- if (plt == TRUE) {- |
-
92 | -1x | -
- .plot(mean_index,- |
-
93 | -1x | -
- paste("Mean cropland connectivity risk index\n"),- |
-
94 | -1x | -
- global,- |
-
95 | -1x | -
- geoscale,- |
-
96 | -1x | -
- zlim = c(0, 1),- |
-
97 | -1x | -
- typ = "mean")- |
-
98 | -- |
- }- |
-
99 | -- | - - | -
100 | -1x | -
- return(mean_index)- |
-
101 | -- |
- }- |
-
102 | -- | - - | -
103 | -- |
- #' Calculate variance of CCRI- |
-
104 | -- |
- #'- |
-
105 | -- |
- #' This function produces a map of variance of CCRI based on input parameters- |
-
106 | -- |
- #' @inheritParams connectivity- |
-
107 | -- |
- #' @param rast A raster object. It will be used in calculating variance.- |
-
108 | -- |
- #' @export- |
-
109 | -- |
- ccri_variance <- function(indexes,- |
-
110 | -- |
- rast,- |
-
111 | -- |
- global,- |
-
112 | -- |
- geoscale,- |
-
113 | -- |
- res = reso()) {- |
-
114 | -- | - - | -
115 | -1x | -
- .cal_var <- function(ext_indices, scale) {- |
-
116 | -2x | -
- var_rastvect <-- |
-
117 | -2x | -
- apply(do.call(cbind, lapply(ext_indices, terra::values)), 1, stats::var, na.rm = TRUE)- |
-
118 | -- | - - | -
119 | -2x | -
- scaled_rast <- terra::crop(rast, .to_ext(scale))- |
-
120 | -2x | -
- scaled_rast[] <- var_rastvect- |
-
121 | -- | - - | -
122 | -2x | -
- var_disag_rast <-- |
-
123 | -2x | -
- terra::disagg(scaled_rast, fact = c(res, res))- |
-
124 | -- | - - | -
125 | -2x | -
- var_disag_rast[var_disag_rast[] == 0] <- NA- |
-
126 | -- | - - | -
127 | -2x | -
- var_disag_rast + .cal_zerorast(var_disag_rast, res)- |
-
128 | -- |
- }- |
-
129 | -- | - - | -
130 | -1x | -
- var_out <- if (global == TRUE) {- |
-
131 | -- | - - | -
132 | -1x | -
- .gan_paramok(indexes)- |
-
133 | -1x | -
- exts <- global_scales()- |
-
134 | -1x | -
- east <-- |
-
135 | -1x | -
- .cal_var(indexes[[STR_EAST]], exts[[STR_EAST]])- |
-
136 | -1x | -
- west <-- |
-
137 | -1x | -
- .cal_var(indexes[[STR_WEST]], exts[[STR_WEST]])- |
-
138 | -- | - - | -
139 | -1x | -
- geoscale <- .global_ext()- |
-
140 | -1x | -
- terra::merge(east, west)- |
-
141 | -- | - - | -
142 | -- |
- } else {- |
-
143 | -! | -
- .cal_var(indexes, geoscale)- |
-
144 | -- |
- }- |
-
145 | -- | - - | -
146 | -1x | -
- z_var_w <- range(var_out[which(var_out[] > 0)])- |
-
147 | -1x | -
- .plot(var_out,- |
-
148 | -1x | -
- "Variance in cropland connectivity",- |
-
149 | -1x | -
- global,- |
-
150 | -1x | -
- geoscale,- |
-
151 | -1x | -
- zlim = z_var_w,- |
-
152 | -1x | -
- typ = "variance")- |
-
153 | -- | - - | -
154 | -1x | -
- invisible(1)- |
-
155 | -- |
- }- |
-
156 | -- | - - | -
157 | -- |
- #' Calculate difference map- |
-
158 | -- |
- #'- |
-
159 | -- |
- #' This function produces a map of difference b/w mean and sum indexes in rank of cropland harvested area fraction.- |
-
160 | -- |
- #' @param rast A raster object for mean index raster difference- |
-
161 | -- |
- #' @param x A raster object for cropland harvest- |
-
162 | -- |
- #' @param y A raster object for cropland harvest- |
-
163 | -- |
- #' @inheritParams connectivity- |
-
164 | -- |
- #' @export- |
-
165 | -- |
- ccri_diff <- function(rast,- |
-
166 | -- |
- x,- |
-
167 | -- |
- y,- |
-
168 | -- |
- global,- |
-
169 | -- |
- geoscale,- |
-
170 | -- |
- res = reso()) {- |
-
171 | -- |
- # difference map- |
-
172 | -- |
- # Function to check for missing or null values- |
-
173 | -1x | -
- .params_ok <- function(...) {- |
-
174 | -1x | -
- !any(sapply(list(...), function(r) missing(r) || is.null(r)))- |
-
175 | -- |
- }- |
-
176 | -- | - - | -
177 | -1x | -
- zr2 <- c(0, 0)- |
-
178 | -- | - - | -
179 | -1x | -
- .cal_diff <- function(r1, r2, scale) {- |
-
180 | -- | - - | -
181 | -2x | -
- scaled_rast <- terra::crop(rast, .to_ext(scale))- |
-
182 | -2x | -
- ccri_id <- which(scaled_rast[] > 0)- |
-
183 | -2x | -
- meantotalland_w <- sum(r1, r2, na.rm = TRUE) / 2- |
-
184 | -- | - - | -
185 | -2x | -
- meanindexcell_w <- scaled_rast[][ccri_id]- |
-
186 | -2x | -
- meantotallandcell_w <- meantotalland_w[][ccri_id]- |
-
187 | -- | - - | -
188 | -- |
- # mean cropland minus mean index,- |
-
189 | -- |
- # negative value means importance of cropland reduce,- |
-
190 | -- |
- # positive value means importance increase,- |
-
191 | -- |
- # zero means the importance of cropland doesn't change.- |
-
192 | -2x | -
- rankdifferent_w <-- |
-
193 | -2x | -
- rank(meantotallandcell_w * (-1)) - rank(meanindexcell_w * (-1))- |
-
194 | -2x | -
- scaled_rast[] <- NaN- |
-
195 | -2x | -
- scaled_rast[][ccri_id] <- rankdifferent_w- |
-
196 | -- | - - | -
197 | -2x | -
- maxrank_w <- max(abs(rankdifferent_w))- |
-
198 | -2x | -
- zr2 <- max(zr2, range(-maxrank_w, maxrank_w))- |
-
199 | -- | - - | -
200 | -2x | -
- diagg_rast <- terra::disagg(scaled_rast,- |
-
201 | -2x | -
- fact = c(res, res))- |
-
202 | -2x | -
- diagg_rast + .cal_zerorast(diagg_rast, res)- |
-
203 | -- |
- }- |
-
204 | -- | - - | -
205 | -1x | -
- diff_out <- if (global == TRUE) {- |
-
206 | -- | - - | -
207 | -1x | -
- if (!.params_ok(the$gan[["sum"]][[STR_EAST]],- |
-
208 | -1x | -
- the$gan[["mean"]][[STR_EAST]],- |
-
209 | -1x | -
- the$gan[["sum"]][[STR_WEST]],- |
-
210 | -1x | -
- the$gan[["mean"]][[STR_WEST]])) {- |
-
211 | -! | -
- message("Either sum or mean aggregate is missing. Aborting difference calculation")- |
-
212 | -! | -
- return(NULL)- |
-
213 | -- |
- }- |
-
214 | -- | - - | -
215 | -1x | -
- exts <- global_scales()- |
-
216 | -- | - - | -
217 | -1x | -
- east_var <-- |
-
218 | -1x | -
- .cal_diff(the$gan[["sum"]][[STR_EAST]], the$gan[["mean"]][[STR_EAST]], exts[[STR_EAST]])- |
-
219 | -1x | -
- west_var <-- |
-
220 | -1x | -
- .cal_diff(the$gan[["sum"]][[STR_WEST]], the$gan[["mean"]][[STR_WEST]], exts[[STR_WEST]])- |
-
221 | -- | - - | -
222 | -1x | -
- geoscale <- .global_ext(exts)- |
-
223 | -1x | -
- terra::merge(east_var, west_var)- |
-
224 | -- |
- } else {- |
-
225 | -! | -
- if (!.params_ok(x, y)) {- |
-
226 | -! | -
- message("Either sum or mean aggregate is missing. Aborting difference calculation")- |
-
227 | -! | -
- return(NULL)- |
-
228 | -- |
- }- |
-
229 | -! | -
- .cal_diff(x, y, geoscale)- |
-
230 | -- |
- }- |
-
231 | -- | - - | -
232 | -1x | -
- .plot(diff_out,- |
-
233 | -1x | -
- "Difference in rank of host connectivity and host density",- |
-
234 | -1x | -
- global,- |
-
235 | -1x | -
- geoscale,- |
-
236 | -1x | -
- .get_palette_for_diffmap(),- |
-
237 | -1x | -
- zr2,- |
-
238 | -1x | -
- typ = "difference")- |
-
239 | -- | - - | -
240 | -1x | -
- invisible()- |
-
241 | -- |
- }- |
-
242 | -- | - - | -
243 | -- |
- # private functions -------------------------------------------------------- |
-
244 | -- | - - | -
245 | -- |
- .plot <- function(rast,- |
-
246 | -- |
- label,- |
-
247 | -- |
- isglobal,- |
-
248 | -- |
- geoscale,- |
-
249 | -- |
- colorss = .get_palette(),- |
-
250 | -- |
- zlim,- |
-
251 | -- |
- typ = "plot",- |
-
252 | -- |
- plotf = .plotmap) {- |
-
253 | -- | - - | -
254 | -3x | -
- .saverast(typ, rast)- |
-
255 | -- | - - | -
256 | -3x | -
- if (is.null(plotf)) {- |
-
257 | -! | -
- .plotmap(rast, geoscale, isglobal, label, colorss, zlim)- |
-
258 | -- |
- } else {- |
-
259 | -3x | -
- plotf(rast = rast,- |
-
260 | -3x | -
- geoscale = geoscale,- |
-
261 | -3x | -
- isglobal = isglobal,- |
-
262 | -3x | -
- label = label,- |
-
263 | -3x | -
- col_pal = colorss,- |
-
264 | -3x | -
- zlim = zlim)- |
-
265 | -- |
- }- |
-
266 | -- | - - | -
267 | -3x | -
- invisible()- |
-
268 | -- |
- }- |
-
269 | -- | - - | -
270 | -- |
- .saverast <- function(typ, rast) {- |
-
271 | -- |
- # Create the "plots" directory if it doesn't exist- |
-
272 | -3x | -
- if (!dir.exists("plots")) {- |
-
273 | -1x | -
- dir.create("plots")- |
-
274 | -- |
- }- |
-
275 | -- | - - | -
276 | -- |
- # Save the plot as a raster file- |
-
277 | -3x | -
- fname <- paste("plots", "/", typ,- |
-
278 | -- |
- "_",- |
-
279 | -3x | -
- stringi::stri_rand_strings(1, 5),- |
-
280 | -3x | -
- ".tif", sep = "")- |
-
281 | -3x | -
- terra::writeRaster(rast, overwrite = TRUE,- |
-
282 | -3x | -
- filename = fname,- |
-
283 | -3x | -
- gdal = c("COMPRESS=NONE"))- |
-
284 | -3x | -
- message(paste("raster created", fname, sep = ": "), "\n")- |
-
285 | -- |
- }- |
-
286 | -- | - - | -
287 | -- |
- .plotmap <- function(rast, geoscale, isglobal, label, col_pal, zlim) {- |
-
288 | -3x | -
- if (interactive()) {- |
-
289 | -- | - - | -
290 | -- |
- # Set the plot parameters- |
-
291 | -! | -
- graphics::par(bg = "aliceblue")- |
-
292 | -- |
- # Plot the base map- |
-
293 | -! | -
- terra::plot(.cal_mgb(geoscale, isglobal),- |
-
294 | -! | -
- col = "grey85",- |
-
295 | -! | -
- xaxt = "n",- |
-
296 | -! | -
- yaxt = "n",- |
-
297 | -! | -
- legend = FALSE,- |
-
298 | -! | -
- main = label,- |
-
299 | -! | -
- cex.main = 0.9)- |
-
300 | -- |
- # Plot the raster- |
-
301 | -! | -
- if (isglobal == TRUE) {- |
-
302 | -! | -
- gs <- terra::ext(rast)- |
-
303 | -! | -
- terra::plot(rast,- |
-
304 | -! | -
- col = col_pal,- |
-
305 | -! | -
- xaxt = "n",- |
-
306 | -! | -
- yaxt = "n",- |
-
307 | -! | -
- zlim = zlim,- |
-
308 | -! | -
- add = TRUE,- |
-
309 | -! | -
- lwd = 0.7,- |
-
310 | -! | -
- legend = TRUE,- |
-
311 | -! | -
- plg = list(loc = "bottom",- |
-
312 | -! | -
- ext = c(gs[1] + 30, gs[2] - 30, gs[3] - 30, gs[3] - 20),- |
-
313 | -! | -
- horizontal = TRUE)- |
-
314 | -- |
- )- |
-
315 | -- |
- } else {- |
-
316 | -! | -
- terra::plot(rast,- |
-
317 | -! | -
- col = col_pal,- |
-
318 | -! | -
- xaxt = "n",- |
-
319 | -! | -
- yaxt = "n",- |
-
320 | -! | -
- zlim = zlim,- |
-
321 | -! | -
- add = TRUE,- |
-
322 | -! | -
- lwd = 0.7,- |
-
323 | -! | -
- legend = TRUE- |
-
324 | -- |
- )- |
-
325 | -- |
- }- |
-
326 | -- | - - | -
327 | -- |
- # plg = list(loc = "bottom", horizontal = TRUE)- |
-
328 | -- |
- # Plot the country boundaries- |
-
329 | -! | -
- world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")- |
-
330 | -! | -
- world <- world[which(world$continent != "Antarctica"), ]["geometry"]- |
-
331 | -! | -
- world <- terra::vect(world)- |
-
332 | -- | - - | -
333 | -! | -
- if (isglobal == FALSE) {- |
-
334 | -! | -
- world <- terra::crop(world, terra::ext(rast))- |
-
335 | -- |
- }- |
-
336 | -- | - - | -
337 | -! | -
- terra::plot(world, col = NA, border = "black", add = TRUE)- |
-
338 | -- |
- }- |
-
339 | -3x | -
- invisible()- |
-
340 | -- |
- }- |
-
341 | -- | - - | -
342 | -- |
- .gan_paramok <- function(indices) {- |
-
343 | -2x | -
- stopifnot("Require list of east and west indices" = all(c(STR_EAST, STR_WEST) %in% names(indices)))- |
-
344 | -2x | -
- return(TRUE)- |
-
345 | -- |
- }- |
-
1 | -- |
- #' Calculate risk index using inbuilt models.- |
-
2 | -- |
- #'- |
-
3 | -- |
- #' @description- |
-
4 | -- |
- #' - [`model_powerlaw()`] calculates risk index using power law.- |
-
5 | -- |
- #' - [`model_neg_exp()`] calculates risk index using negative exponential.- |
-
6 | -- |
- #' @param beta A list of beta values. `DispersalParameterBeta` in `parameters.yaml`.- |
-
7 | -- |
- #' @param gamma_val A list of beta values. `DispersalParameterGamma` in `parameters.yaml`.- |
-
8 | -- |
- #' @param link_threshold A threshold value for link.- |
-
9 | -- |
- #' @param distance_matrix distance matrix, generated during [sean()].- |
-
10 | -- |
- #' @param thresholded_crop_values crop values above threshold.- |
-
11 | -- |
- #' @param adj_mat Adjacency matrix(optional) representing un-directed graph network.- |
-
12 | -- |
- #' If this is provided, then gamma_val, distance_matrix, link_threshold and thresholded_crop_values are ignored.- |
-
13 | -- |
- #' These ignored parameters are used to generate adjacency matrix internally.- |
-
14 | -- |
- #' This is the only way to use custom adjacency matrix.- |
-
15 | -- |
- #' @param crop_raster A raster object for cropland harvest.- |
-
16 | -- |
- #' @param crop_cells_above_threshold crop cells above threshold. Only contains cells and not the the values.- |
-
17 | -- |
- #' @param metrics A list 2 vectors - metrics and weights.- |
-
18 | -- |
- #' @return risk index- |
-
19 | -- |
- #' @details- |
-
20 | -- |
- #' Network metrics should be passed as a list of vectors e.g. `list(metrics = c("betweeness"), weights = c(100))`.- |
-
21 | -- |
- #' Default values are fetched from `parameters.yaml` and arguments uses the same structure.- |
-
22 | -- |
- #'- |
-
23 | -- |
- #' @export- |
-
24 | -- |
- model_powerlaw <- function(beta,- |
-
25 | -- |
- link_threshold,- |
-
26 | -- |
- distance_matrix = the$distance_matrix,- |
-
27 | -- |
- thresholded_crop_values,- |
-
28 | -- |
- adj_mat = NULL,- |
-
29 | -- |
- crop_raster,- |
-
30 | -- |
- crop_cells_above_threshold,- |
-
31 | -- |
- metrics = the$parameters_config$`CCRI parameters`$NetworkMetrics$InversePowerLaw) {- |
-
32 | -- | - - | -
33 | -24x | -
- metrics <- .refined_mets(metrics)- |
-
34 | -- |
- #### create adjacency matrix- |
-
35 | -- | - - | -
36 | -24x | -
- stan <- if (is.null(adj_mat)) {- |
-
37 | -24x | -
- distancematr <- distance_matrix # pairwise distance matrix- |
-
38 | -- |
- #---- end of code- |
-
39 | -- |
- # use function C=AX^(-beta), here A=1, X=distancematr- |
-
40 | -24x | -
- distancematrexp <- distancematr^(-beta)- |
-
41 | -24x | -
- cropmatr <- thresholded_crop_values # complete gravity model with crop data- |
-
42 | -24x | -
- cropmatr1 <- matrix(cropmatr, , 1)- |
-
43 | -24x | -
- cropmatr2 <- matrix(cropmatr, 1, )- |
-
44 | -- | - - | -
45 | -24x | -
- cropmatrix <- cropmatr1 %*% cropmatr2- |
-
46 | -24x | -
- cropmatrix <- as.matrix(cropmatrix)- |
-
47 | -- |
- # adjacency matrix- |
-
48 | -24x | -
- cropdistancematr <- distancematrexp * cropmatrix # make available to users- |
-
49 | -- |
- # adjacency matrix after threshold- |
-
50 | -24x | -
- logicalmatr <- cropdistancematr > link_threshold- |
-
51 | -24x | -
- adjmat <- cropdistancematr * logicalmatr- |
-
52 | -- |
- # use round() because betweeness() may have problem when do the calculation- |
-
53 | -24x | -
- round(adjmat, 6)- |
-
54 | -- |
- } else {- |
-
55 | -! | -
- adj_mat- |
-
56 | -- |
- }- |
-
57 | -- | - - | -
58 | -- |
- # create graph object using adjacency matrix- |
-
59 | -24x | -
- cropdistancematrix <- igraph::graph.adjacency(stan,- |
-
60 | -24x | -
- mode = c("undirected"),- |
-
61 | -24x | -
- diag = FALSE, weighted = TRUE)- |
-
62 | -- | - - | -
63 | -24x | -
- indexpre <- crop_raster- |
-
64 | -24x | -
- indexpre[] <- 0- |
-
65 | -24x | -
- indexpre[crop_cells_above_threshold] <- .apply_met(metrics, cropdistancematrix)- |
-
66 | -24x | -
- indexv <- indexpre- |
-
67 | -24x | -
- return(indexv)- |
-
68 | -- |
- }- |
-
69 | -- | - - | -
70 | -- |
- #' @rdname model_powerlaw- |
-
71 | -- |
- model_neg_exp <- function(gamma_val,- |
-
72 | -- |
- link_threshold,- |
-
73 | -- |
- distance_matrix = the$distance_matrix,- |
-
74 | -- |
- thresholded_crop_values,- |
-
75 | -- |
- adj_mat = NULL,- |
-
76 | -- |
- crop_raster,- |
-
77 | -- |
- crop_cells_above_threshold,- |
-
78 | -- |
- metrics = the$parameters_config$`CCRI parameters`$NetworkMetrics$InversePowerLaw) {- |
-
79 | -32x | -
- metrics <- .refined_mets(metrics)- |
-
80 | -- | - - | -
81 | -- |
- #### create adjacency matrix- |
-
82 | -- |
- ####- |
-
83 | -32x | -
- stan <- if (is.null(adj_mat)) {- |
-
84 | -32x | -
- distancematr <- distance_matrix- |
-
85 | -- | - - | -
86 | -32x | -
- eulernumber <- exp(1)- |
-
87 | -- |
- # exponential model- |
-
88 | -32x | -
- distancematrexponential <-- |
-
89 | -32x | -
- eulernumber^(-gamma_val * distancematr)- |
-
90 | -32x | -
- cropmatr <- thresholded_crop_values # complete gravity model with crop data- |
-
91 | -32x | -
- cropmatr1 <- matrix(cropmatr, , 1) # complete gravity model with crop data- |
-
92 | -32x | -
- cropmatr2 <- matrix(cropmatr, 1, )- |
-
93 | -32x | -
- cropmatrix <- cropmatr1 %*% cropmatr2- |
-
94 | -32x | -
- cropmatrix <- as.matrix(cropmatrix)- |
-
95 | -32x | -
- cropdistancematr <- distancematrexponential * cropmatrix- |
-
96 | -32x | -
- logicalmatr <- cropdistancematr > link_threshold- |
-
97 | -32x | -
- adjmat <- cropdistancematr * logicalmatr- |
-
98 | -- |
- # use round() because betweeness() may have problem when do the calculation- |
-
99 | -32x | -
- round(adjmat, 6)- |
-
100 | -- |
- } else {- |
-
101 | -! | -
- adj_mat- |
-
102 | -- |
- }- |
-
103 | -- | - - | -
104 | -- |
- # create graph object from adjacency matrix- |
-
105 | -32x | -
- cropdistancematrix <- igraph::graph.adjacency(stan,- |
-
106 | -32x | -
- mode = c("undirected"),- |
-
107 | -32x | -
- diag = FALSE, weighted = TRUE- |
-
108 | -- |
- )- |
-
109 | -- | - - | -
110 | -32x | -
- indexpre <- crop_raster- |
-
111 | -32x | -
- indexpre[] <- 0- |
-
112 | -32x | -
- indexpre[crop_cells_above_threshold] <- .apply_met(metrics, cropdistancematrix)- |
-
113 | -32x | -
- indexv <- indexpre- |
-
114 | -32x | -
- return(indexv)- |
-
115 | -- |
- }- |
-
116 | -- | - - | -
117 | -- |
- # private methods ---------------------------------------------------------- |
-
118 | -- | - - | -
119 | -- |
- .apply_met <- function(mets, adj_graph) {- |
-
120 | -- | - - | -
121 | -56x | -
- mets <- Map(c, mets[[1]], mets[[2]])- |
-
122 | -56x | -
- index <- 0- |
-
123 | -- | - - | -
124 | -- |
- # Iterate over the metric names and values in 'mets'- |
-
125 | -56x | -
- for (mname in names(mets)) {- |
-
126 | -224x | -
- if (mname %in% names(metric_funs)) {- |
-
127 | -224x | -
- val <- mets[[mname]][[2]]- |
-
128 | -224x | -
- mfun <- metric_funs[[mname]]- |
-
129 | -224x | -
- index <- index + mfun(adj_graph, val)- |
-
130 | -- |
- }- |
-
131 | -- |
- }- |
-
132 | -- | - - | -
133 | -56x | -
- replace(index, index == 0, NA)- |
-
134 | -- | - - | -
135 | -56x | -
- return(index)- |
-
136 | -- |
- }- |
-
1 | -- |
- .get_param_file_path <- function() {- |
-
2 | -3x | -
- return(system.file("parameters.yaml", package = "geohabnet", mustWork = TRUE))- |
-
3 | -- |
- }- |
-
4 | -- | - - | -
5 | -- |
- .get_directoryfromuser <- function() {- |
-
6 | -! | -
- return(easycsv::choose_dir())- |
-
7 | -- |
- }- |
-
8 | -- | - - | -
9 | -- |
- .open_file_selection_prompt <- function() {- |
-
10 | -! | -
- return(file.choose())- |
-
11 | -- |
- }- |
-
12 | -- | - - | -
13 | -- |
- .copy_file <- function(from, to) {- |
-
14 | -3x | -
- file.copy(from = from, to = to, copy.mode = FALSE, overwrite = TRUE)- |
-
15 | -- |
- }- |
-
16 | -- | - - | -
17 | -- |
- .extract_hosts <- function(params = load_parameters()) {- |
-
18 | -! | -
- monfredas <- params$`CCRI parameters`$Hosts$monfreda- |
-
19 | -! | -
- spams <- params$`CCRI parameters`$Hosts$mapspam- |
-
20 | -! | -
- crops <- list()- |
-
21 | -! | -
- if (!is.null(monfredas) && !is.list(monfredas)) {- |
-
22 | -! | -
- crops[["monfreda"]] <- monfredas- |
-
23 | -- |
- }- |
-
24 | -! | -
- if (!is.null(spams) && !is.list(spams)) {- |
-
25 | -! | -
- crops[["mapspam"]] <- spams- |
-
26 | -- |
- }- |
-
27 | -! | -
- return(crops)- |
-
28 | -- |
- }- |
-
29 | -- | - - | -
30 | -- |
- #' Get Parameters- |
-
31 | -- |
- #'- |
-
32 | -- |
- #' Retrieves the parameters and copies the parameter file to the specified- |
-
33 | -- |
- #' output path.- |
-
34 | -- |
- #' @param out_path character. The output path where the parameter file will be- |
-
35 | -- |
- #' copied. Default is current working directory [getwd()]- |
-
36 | -- |
- #' @param iwindow logical. If `TRUE`, prompts the user to select the output- |
-
37 | -- |
- #' directory using a file chooser window. Default is `FALSE`- |
-
38 | -- |
- #' @return character. The path to the copied parameter file.- |
-
39 | -- |
- #' @export- |
-
40 | -- |
- #' @details- |
-
41 | -- |
- #' Using configuration file is an alternative to [sean()]- |
-
42 | -- |
- #'- |
-
43 | -- |
- #' @seealso [set_parameters()]- |
-
44 | -- |
- get_parameters <- function(out_path = getwd(), iwindow = FALSE) {- |
-
45 | -2x | -
- if (interactive() && iwindow) {- |
-
46 | -! | -
- out_path <- .get_directoryfromuser()- |
-
47 | -- |
- }- |
-
48 | -- | - - | -
49 | -2x | -
- params_file_path <- .get_param_file_path()- |
-
50 | -2x | -
- .copy_file(params_file_path, out_path)- |
-
51 | -2x | -
- print("parameters fetched successfully")- |
-
52 | -2x | -
- return(file.path(paste(out_path, "parameters.yaml", sep = "/")))- |
-
53 | -- |
- }- |
-
54 | -- | - - | -
55 | -- |
- #' Set Parameters- |
-
56 | -- |
- #'- |
-
57 | -- |
- #' This function allows you to set the parameters by replacing the existing- |
-
58 | -- |
- #' parameters file with a new one. Use [get_parameters()] to modify the parameter values.- |
-
59 | -- |
- #'- |
-
60 | -- |
- #' @param new_parameters_file The path to the new parameters file.- |
-
61 | -- |
- #' @param iwindow Logical indicating whether to prompt the user to select the- |
-
62 | -- |
- #' new parameters file using a file selection window. Defaults to FALSE.- |
-
63 | -- |
- #' @return None- |
-
64 | -- |
- #' @export- |
-
65 | -- |
- set_parameters <- function(new_parameters_file, iwindow = FALSE) {- |
-
66 | -1x | -
- if (iwindow && interactive()) {- |
-
67 | -! | -
- new_parameters_file <- .open_file_selection_prompt()- |
-
68 | -- |
- }- |
-
69 | -- | - - | -
70 | -1x | -
- current_params_file <- .get_param_file_path()- |
-
71 | -1x | -
- if (.check_yaml_structure(- |
-
72 | -1x | -
- existing_yaml_file = current_params_file,- |
-
73 | -1x | -
- provided_yaml_file = new_parameters_file- |
-
74 | -- |
- )) {- |
-
75 | -1x | -
- .copy_file(new_parameters_file, current_params_file)- |
-
76 | -- |
- }- |
-
77 | -- |
- }- |
-
78 | -- | - - | -
79 | -- |
- #' Load Parameters from YAML File- |
-
80 | -- |
- #'- |
-
81 | -- |
- #' This function loads parameters from a YAML file and stores them in an object.- |
-
82 | -- |
- #'- |
-
83 | -- |
- #' @param filepath Path to the YAML file containing the parameters. By default, it- |
-
84 | -- |
- #' takes the value of ".kparameters_file_type" which is set to "parameters.yaml".- |
-
85 | -- |
- #'- |
-
86 | -- |
- #' @return object with parameters and values- |
-
87 | -- |
- #'- |
-
88 | -- |
- #' @importFrom config get- |
-
89 | -- |
- #'- |
-
90 | -- |
- #' @examples- |
-
91 | -- |
- #' # Load parameters from default file- |
-
92 | -- |
- #' load_parameters()- |
-
93 | -- |
- #'- |
-
94 | -- |
- #' @export- |
-
95 | -- |
- load_parameters <- function(filepath = .get_helper_filepath(.kparameters_file_type)) {- |
-
96 | -3x | -
- return(config::get(file = filepath))- |
-
97 | -- |
- }- |
-
98 | -- | - - | -
99 | -- |
- #' Get resolution value- |
-
100 | -- |
- #'- |
-
101 | -- |
- #' Resolution stored in `parameter.yaml`.- |
-
102 | -- |
- #' If not present it will result default value.- |
-
103 | -- |
- #' @export- |
-
104 | -- |
- #' @seealso [set_reso()]- |
-
105 | -- |
- reso <- function() {- |
-
106 | -! | -
- reso <- the$parameters_config$`CCRI parameters`$Resolution- |
-
107 | -! | -
- reso <- if (is.null(reso) || is.na(reso)) {- |
-
108 | -! | -
- 24- |
-
109 | -- |
- } else {- |
-
110 | -! | -
- reso- |
-
111 | -- |
- }- |
-
112 | -! | -
- return(reso)- |
-
113 | -- |
- }- |
-
114 | -- | - - | -
115 | -- |
- #' Set resolution value- |
-
116 | -- |
- #'- |
-
117 | -- |
- #' Set `resolution` to be used in analysis.- |
-
118 | -- |
- #' It doesn't modify the `parameters.yaml`- |
-
119 | -- |
- #' but instead a currently loaded instance of it.- |
-
120 | -- |
- #' Must be greater than 0 and less than or equal to 48.- |
-
121 | -- |
- #'- |
-
122 | -- |
- #' @param value numeric. Resolution value.- |
-
123 | -- |
- #' @export- |
-
124 | -- |
- #' @examples- |
-
125 | -- |
- #' \dontrun{- |
-
126 | -- |
- #' set_reso(24)- |
-
127 | -- |
- #' }- |
-
128 | -- |
- #'- |
-
129 | -- |
- set_reso <- function(value) {- |
-
130 | -! | -
- stopifnot("Invalid resolution" = is.numeric(value), value <= 0, value >= 48)- |
-
131 | -! | -
- .loadparam_ifnull()- |
-
132 | -! | -
- the$parameters_config$`CCRI parameters`$Resolution <- value- |
-
133 | -! | -
- invisible(TRUE)- |
-
134 | -- |
- }- |
-
1 | -- |
- #' Only meant to global variables- |
-
2 | -- |
- #' @keywords internal- |
-
3 | -- |
- scales <- new.env(parent = emptyenv())- |
-
4 | -- |
- scales$geast <- c(-24, 180, -58, 60)- |
-
5 | -- |
- scales$gwest <- c(-140, -34, -58, 60)- |
-
6 | -- | - - | -
7 | -- |
- #' Global geographical extent- |
-
8 | -- |
- #'- |
-
9 | -- |
- #' See geographical extents used in global analysis.- |
-
10 | -- |
- #' Returns eastern and western hemisphere extents.- |
-
11 | -- |
- #' Each extent is in the form of c(Xmin, Xmax, Ymin, Ymax).- |
-
12 | -- |
- #' @export- |
-
13 | -- |
- #' @seealso [set_global_scales()]- |
-
14 | -- |
- global_scales <- function() {- |
-
15 | -7x | -
- sc <- list(scales$geast, scales$gwest)- |
-
16 | -7x | -
- names(sc) <- c(STR_EAST, STR_WEST)- |
-
17 | -7x | -
- return(sc)- |
-
18 | -- |
- }- |
-
19 | -- | - - | -
20 | -- |
- .global_ext <- function(scales = global_scales()) {- |
-
21 | -4x | -
- return(scales[[STR_EAST]] + scales[[STR_WEST]])- |
-
22 | -- |
- }- |
-
23 | -- | - - | -
24 | -- |
- #' Set global geographical extent- |
-
25 | -- |
- #'- |
-
26 | -- |
- #' Set the geographical extents used in global analysis.- |
-
27 | -- |
- #' Each extent should be in the form of c(Xmin, Xmax, Ymin, Ymax)- |
-
28 | -- |
- #' @param value list. Named list of eastern and western hemisphere extents. See usage.- |
-
29 | -- |
- #' @export- |
-
30 | -- |
- #' @examples- |
-
31 | -- |
- #' \dontrun{- |
-
32 | -- |
- #' set_global_scales(list(east = c(-24, 180, -58, 60), west = c(-140, -34, -58, 60)))- |
-
33 | -- |
- #' }- |
-
34 | -- |
- #' @seealso- |
-
35 | -- |
- #' [global_scales()]- |
-
36 | -- |
- #' [terra::ext()]- |
-
37 | -- |
- set_global_scales <- function(value) {- |
-
38 | -- | - - | -
39 | -! | -
- east <- is.vector(value[[STR_EAST]])- |
-
40 | -! | -
- west <- is.vector(value[[STR_WEST]])- |
-
41 | -! | -
- stopifnot("Not a valid east coordinate" = is.numeric(east), is.vector(east), length(east) == 4)- |
-
42 | -! | -
- stopifnot("Not a valid west coordinate" = is.numeric(west), is.vector(west), length(west) == 4)- |
-
43 | -- | - - | -
44 | -! | -
- scales$geast <- east- |
-
45 | -! | -
- scales$gwest <- west- |
-
46 | -! | -
- value <- c(east, west)- |
-
47 | -! | -
- return(invisible(value))- |
-
48 | -- |
- }- |
-
49 | -- | - - | -
50 | -- |
- #' Get geographical scales from the parameters- |
-
51 | -- |
- #'- |
-
52 | -- |
- #' This function returns a list of geographical scales set in global and custom extent in `parameters.yaml`.- |
-
53 | -- |
- #' If `global` is `TRUE`, the `CustomExt` is ignored.- |
-
54 | -- |
- #' @return A list of geographical scales- |
-
55 | -- |
- #' @export- |
-
56 | -- |
- geoscale_param <- function() {- |
-
57 | -1x | -
- .loadparam_ifnull()- |
-
58 | -1x | -
- xf <- the$parameters_config$`CCRI parameters`$GeoExtent$global- |
-
59 | -1x | -
- if (as.logical(xf) == FALSE) {- |
-
60 | -! | -
- stopifnot("Geographical missing in parameters " =- |
-
61 | -! | -
- length(the$parameters_config$`CCRI parameters`$GeoExtent$customExt) == 4)- |
-
62 | -- |
- }- |
-
63 | -1x | -
- return(as.vector(the$parameters_config$`CCRI parameters`$GeoExtent$customExt))- |
-
64 | -- |
- }- |
-
1 | -- |
- #' Get rasters object from parameters- |
-
2 | -- |
- #'- |
-
3 | -- |
- #' Takes named list of hosts as an input. See host object in [get_parameters()] or [load_parameters()].- |
-
4 | -- |
- #' Function creates 2 raster object - one is a sum of all the crops specified under sources- |
-
5 | -- |
- #' and other using the provided raster file. See [get_crop_raster_fromtif()]- |
-
6 | -- |
- #' @param hosts List of hosts and values. It is synonym to Hosts object in parameters- |
-
7 | -- |
- #' @return List of rasters- |
-
8 | -- |
- #' @examples- |
-
9 | -- |
- #' # Get default rasters- |
-
10 | -- |
- #' \dontrun{- |
-
11 | -- |
- #' get_rasters(list(mapspam = c("wheat"), monfreda = c("avocado"), file = "some_raster.tif"))- |
-
12 | -- |
- #' }- |
-
13 | -- |
- #' @seealso [load_parameters()], [get_parameters()], [get_crop_raster_fromtif()], [get_cropharvest_raster()]- |
-
14 | -- |
- get_rasters <- function(hosts) {- |
-
15 | -3x | -
- t_file <- hosts[["file"]]- |
-
16 | -3x | -
- rasters <- list()- |
-
17 | -3x | -
- if (!is.null(t_file) &&- |
-
18 | -3x | -
- !is.na(t_file)) {- |
-
19 | -2x | -
- rasters <- c(rasters, get_crop_raster_fromtif(t_file))- |
-
20 | -- |
- }- |
-
21 | -2x | -
- rasters <- c(rasters, get_cropharvest_raster_sum(hosts))- |
-
22 | -2x | -
- return(rasters)- |
-
23 | -- |
- }- |
-
24 | -- | - - | -
25 | -- | - - | -
26 | -- |
- #' @title Get sum of rasters for individual crops- |
-
27 | -- |
- #'- |
-
28 | -- |
- #' @description- |
-
29 | -- |
- #' Takes crop names and returns raster object which is sum of raster of individual crops.- |
-
30 | -- |
- #' Currently, only supports crops listed in- |
-
31 | -- |
- #' [geodata::monfredaCrops()], [geodata::spamCrops()]- |
-
32 | -- |
- #' If crop is present in multiple sources, then their mean is calculated.- |
-
33 | -- |
- #' @param crop_names A named list of source along with crop names- |
-
34 | -- |
- #' @return Raster object which is sum of all the individual crop rasters- |
-
35 | -- |
- #' @export- |
-
36 | -- |
- #' @examples- |
-
37 | -- |
- #' \dontrun{- |
-
38 | -- |
- #' get_cropharvest_raster_sum(list(monfreda = c("wheat", "barley"), mapspam = c("wheat", "potato")))- |
-
39 | -- |
- #' }- |
-
40 | -- |
- get_cropharvest_raster_sum <- function(crop_names) {- |
-
41 | -2x | -
- if (!is.list(crop_names) || length(crop_names) == 0) {- |
-
42 | -! | -
- stop("Input 'crop_names' must be a non-empty list of crop names.")- |
-
43 | -- |
- }- |
-
44 | -- | - - | -
45 | -- |
- # output: list("wheat" = c("monfreda", "mapspam"), "barley" = c("monfreda"), "potato" = c("mapspam"))- |
-
46 | -2x | -
- crops <- list()- |
-
47 | -- | - - | -
48 | -2x | -
- for (src in get_supported_sources()) {- |
-
49 | -4x | -
- for (crop_name in crop_names[[src]]) {- |
-
50 | -1x | -
- crops[[crop_name]] <- c(crops[[crop_name]], src)- |
-
51 | -- |
- }- |
-
52 | -- |
- }- |
-
53 | -- | - - | -
54 | -- |
- # calculate sum of rasters- |
-
55 | -- | - - | -
56 | -- |
- # iterate named lists- |
-
57 | -- |
- # crop names- |
-
58 | -2x | -
- nams <- names(crops)- |
-
59 | -2x | -
- cropharvests <- list()- |
-
60 | -2x | -
- for (i in seq_along(crops)) {- |
-
61 | -1x | -
- single_crop_rasters <- list()- |
-
62 | -1x | -
- for (j in crops[[i]]) {- |
-
63 | -1x | -
- cr <- get_cropharvest_raster(nams[i], j)- |
-
64 | -1x | -
- single_crop_rasters <- c(single_crop_rasters, cr)- |
-
65 | -- |
- }- |
-
66 | -1x | -
- len_scr <- length(single_crop_rasters)- |
-
67 | -1x | -
- if (len_scr > 1) {- |
-
68 | -! | -
- cropharvests <- c(cropharvests, terra::app(terra::rast(single_crop_rasters), fun = sum, na.rm = TRUE) / len_scr)- |
-
69 | -- |
- } else {- |
-
70 | -1x | -
- cropharvests <- c(cropharvests, single_crop_rasters)- |
-
71 | -- |
- }- |
-
72 | -- |
- }- |
-
73 | -- | - - | -
74 | -2x | -
- return(Reduce("+", cropharvests))- |
-
75 | -- |
- }- |
-
76 | -- | - - | -
77 | -- |
- #' @title Get raster object for crop- |
-
78 | -- |
- #' @description Get cropland information in a form of raster object from data source for crop- |
-
79 | -- |
- #' @param crop_name Name of the crop- |
-
80 | -- |
- #' @param data_source Data source for cropland information- |
-
81 | -- |
- #' @return Raster object- |
-
82 | -- |
- #' @export- |
-
83 | -- |
- #' @examples- |
-
84 | -- |
- #' get_cropharvest_raster("avocado", "monfreda")- |
-
85 | -- |
- get_cropharvest_raster <- function(crop_name, data_source) {- |
-
86 | -- |
- # supported sources- |
-
87 | -2x | -
- sources <- get_supported_sources()- |
-
88 | -2x | -
- if (!(data_source %in% sources)) {- |
-
89 | -! | -
- stop(paste("data source: ", data_source, " is not supported"))- |
-
90 | -- |
- }- |
-
91 | -2x | -
- cropharvest_r <- .get_cropharvest_raster_helper(crop_name = crop_name, data_source = data_source)- |
-
92 | -- | - - | -
93 | -2x | -
- return(cropharvest_r)- |
-
94 | -- |
- }- |
-
95 | -- | - - | -
96 | -- |
- #' Get raster object from tif file- |
-
97 | -- |
- #'- |
-
98 | -- |
- #' This is a wrapper of [terra::rast()] and generates a raster object if provided with a TIF file.- |
-
99 | -- |
- #'- |
-
100 | -- |
- #' @param path_to_tif TIF file- |
-
101 | -- |
- #' @return Raster object- |
-
102 | -- |
- #' @examples- |
-
103 | -- |
- #' \dontrun{- |
-
104 | -- |
- #' # Generate raster for usage- |
-
105 | -- |
- #' fp <- .get_helper_filepath("avocado")- |
-
106 | -- |
- #' get_crop_raster_fromtif(fp)- |
-
107 | -- |
- #' get_crop_raster_fromtif("avocado_HarvestedAreaFraction.tif")- |
-
108 | -- |
- #'- |
-
109 | -- |
- #' }- |
-
110 | -- |
- get_crop_raster_fromtif <- function(path_to_tif) {- |
-
111 | -2x | -
- .validate_tif(path_to_tif)- |
-
112 | -1x | -
- return(terra::rast(path_to_tif))- |
-
113 | -- |
- }- |
-
114 | -- | - - | -
115 | -- |
- .validate_tif <- function(path_to_tif) {- |
-
116 | -2x | -
- stopifnot(- |
-
117 | -2x | -
- file.exists(path_to_tif),- |
-
118 | -2x | -
- stringr::str_sub(path_to_tif, start = -4) == ".tif"- |
-
119 | -- |
- )- |
-
120 | -- |
- }- |
-
1 | -- |
- # exported functions ------------------------------------------------------- |
-
2 | -- | - - | -
3 | -- |
- #' list stores functions to apply metrics to distance metrics.- |
-
4 | -- |
- #' @keywords internal- |
-
5 | -- |
- metric_funs <- list(- |
-
6 | -- |
- STR_NEAREST_NEIGHBORS_SUM = function(graph, param) nn_sum(graph, param),- |
-
7 | -- |
- STR_NODE_STRENGTH = function(graph, param) node_strength(graph, param),- |
-
8 | -- |
- STR_BETWEENNESS = function(graph, param) betweeness(graph, param),- |
-
9 | -- |
- STR_EIGEN_VECTOR_CENTRALITY = function(graph, param) ev(graph, param),- |
-
10 | -- |
- STR_CLOSENESS_CENTRALITY = function(graph, param) closeness(graph, param),- |
-
11 | -- |
- STR_PAGE_RANK = function(graph, param) pagerank(graph, param),- |
-
12 | -- |
- STR_DEGREE = function(graph, param) degree(graph, param)- |
-
13 | -- |
- )- |
-
14 | -- | - - | -
15 | -- |
- #'- |
-
16 | -- |
- #' Returns metrics currently supported in the analysis.- |
-
17 | -- |
- #'- |
-
18 | -- |
- #'- |
-
19 | -- |
- #' @returns vector of supported metrics.- |
-
20 | -- |
- #' @examples- |
-
21 | -- |
- #' supported_metrics()- |
-
22 | -- |
- #'- |
-
23 | -- |
- #' @family metrics- |
-
24 | -- |
- #' @export- |
-
25 | -- |
- supported_metrics <- function() {- |
-
26 | -60x | -
- return(c(STR_BETWEENNESS,- |
-
27 | -60x | -
- STR_NODE_STRENGTH,- |
-
28 | -60x | -
- STR_NEAREST_NEIGHBORS_SUM,- |
-
29 | -60x | -
- STR_EIGEN_VECTOR_CENTRALITY,- |
-
30 | -60x | -
- STR_CLOSENESS_CENTRALITY,- |
-
31 | -60x | -
- STR_DEGREE,- |
-
32 | -60x | -
- STR_PAGE_RANK))- |
-
33 | -- |
- }- |
-
34 | -- | - - | -
35 | -- |
- #' Get metrics from parameters- |
-
36 | -- |
- #'- |
-
37 | -- |
- #' Get metrics and parameters stored in `parameters.yaml`.- |
-
38 | -- |
- #'- |
-
39 | -- |
- #' @param params R object of [load_parameters()]. Default is `load_parameters()`.- |
-
40 | -- |
- #'- |
-
41 | -- |
- #' @return List of metrics - parameters and values. See usage.- |
-
42 | -- |
- #'- |
-
43 | -- |
- #' @examples- |
-
44 | -- |
- #' # Get metrics from parameters- |
-
45 | -- |
- #' get_param_metrics()- |
-
46 | -- |
- #' get_param_metrics(load_parameters())- |
-
47 | -- |
- #'- |
-
48 | -- |
- #' @export- |
-
49 | -- |
- get_param_metrics <- function(params = load_parameters()) {- |
-
50 | -- | - - | -
51 | -2x | -
- return(list(pl = .refined_mets(params$`CCRI parameters`$NetworkMetrics$InversePowerLaw),- |
-
52 | -2x | -
- ne = .refined_mets(params$`CCRI parameters`$NetworkMetrics$NegativeExponential)))- |
-
53 | -- |
- }- |
-
54 | -- | - - | -
55 | -- |
- #' Calculation on network matrix.- |
-
56 | -- |
- #'- |
-
57 | -- |
- #' @description- |
-
58 | -- |
- #' These are basically an abstraction of functions under the [igraph] package.- |
-
59 | -- |
- #' The functions included in this abstraction are:- |
-
60 | -- |
- #' - `nn_sum()`: Calculates the sum of nearest neighbors [igraph::graph.knn()].- |
-
61 | -- |
- #' - `node_strength()`: Calculates the sum of edge weights of adjacent nodes [igraph::graph.strength()].- |
-
62 | -- |
- #' - `betweeness()`: Calculates the vertex and edge betweenness based on the number of geodesics- |
-
63 | -- |
- #' [igraph::betweenness()].- |
-
64 | -- |
- #' - `ev()`: Calculates the eigenvector centrality of positions within the network [igraph::evcent()].- |
-
65 | -- |
- #' - `closeness()`: measures how many steps is required to access every other vertex from a given vertex- |
-
66 | -- |
- #' [igraph::closeness()].- |
-
67 | -- |
- #' - `degree()`: number of adjacent edges [igraph::degree()].- |
-
68 | -- |
- #' - `page_rank()`: page rank score for vertices [igraph::page_rank()].- |
-
69 | -- |
- #' @param crop_dm Distance matrix.- |
-
70 | -- |
- #' In the internal workflow,- |
-
71 | -- |
- #' the distance matrix comes is a result of operations within [sean()] and risk functions.- |
-
72 | -- |
- #' @param we Weight in percentage.- |
-
73 | -- |
- #' @return Matrix with the mean value based on the assigned weight.- |
-
74 | -- |
- #'- |
-
75 | -- |
- #' @family metrics- |
-
76 | -- |
- #' @export- |
-
77 | -- |
- nn_sum <- function(crop_dm, we) {- |
-
78 | -- | - - | -
79 | -56x | -
- knnpref0 <- igraph::graph.knn(crop_dm, weights = NA)$knn- |
-
80 | -56x | -
- knnpref0[is.na(knnpref0)] <- 0- |
-
81 | -56x | -
- degreematr <- igraph::degree(crop_dm)- |
-
82 | -56x | -
- knnpref <- knnpref0 * degreematr- |
-
83 | -56x | -
- knnprefp <- if (max(knnpref) == 0) {- |
-
84 | -4x | -
- 0- |
-
85 | -56x | -
- } else if (max(knnpref) > 0) {- |
-
86 | -52x | -
- (knnpref / max(knnpref)) * .per_to_real(we)- |
-
87 | -- |
- }- |
-
88 | -56x | -
- return(knnprefp)- |
-
89 | -- |
- }- |
-
90 | -- | - - | -
91 | -- | - - | -
92 | -- |
- #' @rdname nn_sum- |
-
93 | -- |
- node_strength <- function(crop_dm, we) {- |
-
94 | -56x | -
- nodestrength <- igraph::graph.strength(crop_dm)- |
-
95 | -56x | -
- nodestrength[is.na(nodestrength)] <- 0- |
-
96 | -56x | -
- nodestr <- if (max(nodestrength) == 0) {- |
-
97 | -4x | -
- 0- |
-
98 | -56x | -
- } else if (max(nodestrength) > 0) {- |
-
99 | -52x | -
- (nodestrength / max(nodestrength)) * .per_to_real(we)- |
-
100 | -- |
- }- |
-
101 | -- | - - | -
102 | -56x | -
- return(nodestr)- |
-
103 | -- |
- }- |
-
104 | -- | - - | -
105 | -- |
- #' @rdname nn_sum- |
-
106 | -- |
- betweeness <- function(crop_dm, we) {- |
-
107 | -- | - - | -
108 | -56x | -
- between <- igraph::betweenness(crop_dm,- |
-
109 | -56x | -
- weights =- |
-
110 | -56x | -
- (1 - 1 / exp(.get_weight_vector(- |
-
111 | -56x | -
- crop_dm- |
-
112 | -- |
- ))))- |
-
113 | -- | - - | -
114 | -56x | -
- between[is.na(between)] <- 0- |
-
115 | -56x | -
- betweenp <- if (max(between) == 0) {- |
-
116 | -16x | -
- 0- |
-
117 | -56x | -
- } else if (max(between) > 0) {- |
-
118 | -40x | -
- (between / max(between)) * .per_to_real(we)- |
-
119 | -- |
- }- |
-
120 | -56x | -
- return(betweenp)- |
-
121 | -- |
- }- |
-
122 | -- | - - | -
123 | -- |
- #' @rdname nn_sum- |
-
124 | -- |
- ev <- function(crop_dm, we) {- |
-
125 | -56x | -
- eigenvectorvalues <- igraph::evcent(crop_dm)- |
-
126 | -56x | -
- evv <- eigenvectorvalues$vector- |
-
127 | -56x | -
- evv[is.na(evv)] <- 0- |
-
128 | -56x | -
- evp <- if (max(evv) == 0) {- |
-
129 | -! | -
- 0- |
-
130 | -- |
- } else {- |
-
131 | -56x | -
- (evv / max(evv)) * .per_to_real(we)- |
-
132 | -- |
- }- |
-
133 | -- | - - | -
134 | -56x | -
- return(evp)- |
-
135 | -- |
- }- |
-
136 | -- | - - | -
137 | -- |
- #' @rdname nn_sum- |
-
138 | -- |
- degree <- function(crop_dm, we) {- |
-
139 | -! | -
- dmat <- igraph::degree(crop_dm)- |
-
140 | -! | -
- dmat[is.na(dmat)] <- 0- |
-
141 | -! | -
- dmatr <- if (max(dmat) == 0) {- |
-
142 | -! | -
- 0- |
-
143 | -- |
- } else {- |
-
144 | -! | -
- (dmat / max(dmat)) * .per_to_real(we)- |
-
145 | -- |
- }- |
-
146 | -! | -
- return(dmatr)- |
-
147 | -- |
- }- |
-
148 | -- | - - | -
149 | -- |
- #' @rdname nn_sum- |
-
150 | -- |
- closeness <- function(crop_dm, we) {- |
-
151 | -! | -
- cvv <- igraph::closeness(crop_dm)- |
-
152 | -! | -
- cvv[is.na(cvv)] <- 0- |
-
153 | -! | -
- cns <- if (max(cvv) == 0) {- |
-
154 | -! | -
- 0- |
-
155 | -- |
- } else {- |
-
156 | -! | -
- (cvv / max(cvv)) * .per_to_real(we)- |
-
157 | -- |
- }- |
-
158 | -! | -
- return(cns)- |
-
159 | -- |
- }- |
-
160 | -- | - - | -
161 | -- |
- #' @rdname nn_sum- |
-
162 | -- |
- page_rank <- function(crop_dm, we) {- |
-
163 | -! | -
- pr_scores <- igraph::page_rank(crop_dm)- |
-
164 | -! | -
- prv <- pr_scores$vector- |
-
165 | -! | -
- prv[is.na(prv)] <- 0- |
-
166 | -! | -
- prv <- if (max(prv) == 0) {- |
-
167 | -! | -
- 0- |
-
168 | -- |
- } else {- |
-
169 | -! | -
- (prv / max(prv)) * .per_to_real(we)- |
-
170 | -- |
- }- |
-
171 | -! | -
- return(prv)- |
-
172 | -- |
- }- |
-
173 | -- | - - | -
174 | -- |
- # ------------------------------------------Private methods------------------------------------------- |
-
175 | -- | - - | -
176 | -- |
- .validate_weights <- function(me, we) {- |
-
177 | -60x | -
- stopifnot("Sum of metric weights should be 100" = sum(we) == 100)- |
-
178 | -60x | -
- stopifnot("Weights or metrics missing. Each metric should have a weight" = length(me) == length(we))- |
-
179 | -- |
- }- |
-
180 | -- | - - | -
181 | -- |
- .validate_metrics <- function(metrics) {- |
-
182 | -- |
- # check if weights are valid- |
-
183 | -60x | -
- .validate_weights(metrics[[1]], metrics[[2]])- |
-
184 | -- | - - | -
185 | -60x | -
- not_sup <- metrics[[1]][!metrics[[1]] %in% supported_metrics()]- |
-
186 | -60x | -
- if (length(not_sup) > 0) {- |
-
187 | -! | -
- stop(paste("Following metrics are not supported: ", paste(not_sup, collapse = ", ")))- |
-
188 | -- |
- }- |
-
189 | -60x | -
- invisible()- |
-
190 | -- |
- }- |
-
191 | -- | - - | -
192 | -- |
- .per_to_real <- function(we) {- |
-
193 | -200x | -
- return(as.numeric(we) / 100)- |
-
194 | -- |
- }- |
-
195 | -- | - - | -
196 | -- |
- .list_to_vec <- function(metrics) {- |
-
197 | -! | -
- return(do.call(cbind, metrics))- |
-
198 | -- |
- }- |
-
199 | -- | - - | -
200 | -- |
- .get_weight_vector <- function(cropdistancematrix) {- |
-
201 | -56x | -
- weight_vec <- igraph::E(cropdistancematrix)$weight- |
-
202 | -56x | -
- weight_vec[is.na(weight_vec)] <- 0- |
-
203 | -56x | -
- weight_vec <- weight_vec + 1e-10- |
-
204 | -56x | -
- return(weight_vec)- |
-
205 | -- |
- }- |
-
206 | -- | - - | -
207 | -- |
- .refined_mets <- function(mets) {- |
-
208 | -60x | -
- mets <- lapply(mets, tolower)- |
-
209 | -60x | -
- mets[["weights"]] <- as.numeric(mets[["weights"]])- |
-
210 | -60x | -
- if (!is.null(mets)) {- |
-
211 | -60x | -
- .validate_metrics(mets)- |
-
212 | -- |
- }- |
-
213 | -60x | -
- return(mets)- |
-
214 | -- |
- }- |
-
"+y.value+"
";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p