Skip to content

Commit

Permalink
Merge pull request #9 from Risk-Team/Dev
Browse files Browse the repository at this point in the history 
v3.2
  • Loading branch information
@RSO9192
RSO9192 committed May 15, 2024
2 parents 5eef86a + 44ec8d2 commit 3bda49d
Show file tree
Hide file tree
Showing 60 changed files with 563 additions and 262 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: CAVAanalytics
Type: Package
Title: R package providing a framework for easy access, processing, and advanced visualization of gridded climate products
Version: 3.1.0
Version: 3.2.0
Authors@R: c(person("Riccardo", "Soldan", , "riccardosoldan@hotmail.it", role = c("aut", "cre"),
comment = c(ORCID = "0000-0003-2583-5840")),
person("Rodrigo", "Manzanas", , "rodrigo.manzanas@unican.es", role = "aut",
Expand Down
11 changes: 10 additions & 1 deletion NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,18 @@ output:
toc_depth: 2
---

## CAVAanalytics 3.2.0

In this new version of CAVAanalytics we improve the bias correction funcionalities of CAVAanalytics. Specifically:

- the scaling method was included among bias-correction options
- Cross-validation was added as an option in model_biases to take care of overfitting
- The window argument allows to bias-correct data monthly or on annual basis
- The memory efficient functions performs interpolation when needed to merge spatial chunks with different resolutions

## CAVAanalytics 3.1.0

In this new version of CAVAanalytics, many improvments have been made:
In this new version of CAVAanalytics, many improvements have been made:

- Users can now select which bias correction method to use (default to Empirical Quantile Mapping)
- Intervals can be specified for customizing breaks in the color palette
Expand Down
70 changes: 57 additions & 13 deletions R/climate_change_signal.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,8 @@
#' @param bias.correction logical
#' @param threshold numerical value with range 0-1. It indicates the threshold for assigning model agreement. For example, 0.6 indicates that model agreement is assigned when 60 percent of the models agree in the sign of the change
#' @param n.sessions numeric, number of sessions to use, default is one. Parallelisation can be useful when multiple scenarios are used (RCPS, SSPs). However, note that parallelising will increase RAM usage
#' @param method character, bias-correction method to use. One of eqm (Empirical Quantile Mapping) or qdm (Quantile Delta Mapping). Default to eqm
#' @param method character, bias-correction method to use. One of eqm (Empirical Quantile Mapping), qdm (Quantile Delta Mapping) or scaling. Default to eqm. When using the scaling method, the multiplicative approach is automatically applied only when the variable is precipitation.
#' @param window character, one of none or monthly. Whether bias correction should be applied on a monthly or annual basis. Monthly is the preferred option when performing bias-correction using daily data
#' @param percentage logical, whether the climate change signal is to be calculated as relative changes (in percentage). Default to FALSE
#' @importFrom magrittr %>%
#' @return list with SpatRaster. To explore the output run attributes(output)
Expand All @@ -30,7 +31,8 @@ climate_change_signal <- function(data,
threshold = 0.6,
n.sessions = 1,
method = "eqm",
percentage = F) {
percentage = F,
window = "monthly") {
# Intermediate functions --------------------------------------------------

# check inputs requirement
Expand All @@ -44,7 +46,8 @@ climate_change_signal <- function(data,
season,
agreement,
method,
percentage) {
percentage,
window) {
stopifnot(is.logical(consecutive))
stopifnot(is.logical(percentage))
stopifnot(is.numeric(threshold), threshold >= 0, threshold <= 1)
Expand All @@ -53,9 +56,24 @@ climate_change_signal <- function(data,
if (!(method == "eqm" || method == "qdm")) {
cli::cli_abort("method must be 'eqm' or qdm")
}
if (!(window == "none" || window == "monthly")) {
cli::cli_abort("window must be one of 'none', or 'monthly'")
}
if (!(duration == "max" || is.numeric(duration))) {
cli::cli_abort("duration must be 'max' or a number")
}
if (bias.correction & window == "monthly") {
dates <- data[[1]]$models_mbrs[[1]]$Dates$start
dates <- as.Date(dates)
# calculate the differences between consecutive dates
diffs <- diff(dates)
# check if the differences are equal to 1
if (any(diffs == 1)) {

} else {
cli::cli_abort(c("x" = "Data is monthly or greater, set window to none"))
}
}
if (!any(stringr::str_detect(colnames(data[[1]]), "obs")) &
isTRUE(bias.correction)) {
cli::cli_abort(
Expand All @@ -66,9 +84,9 @@ climate_change_signal <- function(data,
!is.null(uppert))
cli::cli_abort(c("x" = "Specify only one threshold argument"))
if ((is.null(lowert) &
is.null(uppert)) & bias.correction)
is.null(uppert)) & bias.correction & method == "scaling")
cli::cli_abort(
c("x" = "Bias correction can change the results of the climate change signal only for the calculation of indicators. Specify lowert or uppert aguments to use this option")
c("x" = "Bias correction with the scaling method can change the results of the climate change signal only for the calculation of indicators. Specify lowert or uppert aguments to use this option")
)
if (!is.null(lowert) |
!is.null(uppert)) {
Expand All @@ -83,6 +101,20 @@ climate_change_signal <- function(data,
cli::cli_abort(c("x" = "Data is monthly or greater, thresholds cannot be calculated. Set as NULL"))
}
}

if (bias.correction & window == "monthly") {
dates <- data[[1]]$models_mbrs[[1]]$Dates$start
dates <- as.Date(dates)
# calculate the differences between consecutive dates
diffs <- diff(dates)
# check if the differences are equal to 1
if (any(diffs == 1)) {

} else {
cli::cli_abort(c("x" = "Data is monthly or greater, set window to none"))
}
}

if (consecutive &
is.null(uppert) &
is.null(lowert))
Expand Down Expand Up @@ -203,7 +235,8 @@ climate_change_signal <- function(data,
frequency,
threshold,
method,
percentage) {
percentage,
window) {
season_name <-
convert_vector_to_month_initials(season)
data_list <- datasets %>%
Expand All @@ -212,9 +245,11 @@ climate_change_signal <- function(data,
cli::cli_text(
paste(
"{cli::symbol$arrow_right}",
" Performing bias correction with the ",
" Performing ",
ifelse(window == "monthly", "monthly", ""),
" bias correction with the ",
method,
" method, for each model and month separately. This can take a while. Season",
" method, for each model separately. This can take a while. Season",
glue::glue_collapse(season, "-")
)
)
Expand All @@ -226,12 +261,19 @@ climate_change_signal <- function(data,
downscaleR::biasCorrection(
y = obs[[1]],
x = mod,
precipitation = ifelse(var == "pr", TRUE, FALSE),
precipitation = if (var == "pr")
TRUE
else
FALSE,
method = method,
window = if (any(diffs == 1))
scaling.type = if (var == "pr")
"multiplicative"
else
"additive",
window = if (window == "monthly")
c(30, 30)
else
c(1, 1),
NULL,
extrapolation = "constant"
)
)
Expand Down Expand Up @@ -426,7 +468,8 @@ climate_change_signal <- function(data,
season,
threshold,
method,
percentage
percentage,
window
)

# retrieve information
Expand Down Expand Up @@ -486,7 +529,8 @@ climate_change_signal <- function(data,
frequency,
threshold,
method,
percentage
percentage,
window
)
cli::cli_progress_done()
# return results
Expand Down
55 changes: 32 additions & 23 deletions R/load_data_and_climate_change_signal.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,8 @@
#' @param threshold numerical value with range 0-1. It indicates the threshold for assigning model agreement. For example, 0.6 indicates that model agreement is assigned when 60 percent of the models agree in the sign of the change
#' @param overlap numeric, amount of overlap needed to create the composite. Default 0.25
#' @param percentage logical, whether the climate change signal is to be calculated as relative changes (in percentage). Default to FALSE
#' @param method character, bias-correction method to use. One of eqm (Empirical Quantile Mapping) or qdm (Quantile Delta Mapping). Default to eqm
#' @param method character, bias-correction method to use. One of eqm (Empirical Quantile Mapping), qdm (Quantile Delta Mapping) or scaling. Default to eqm. When using the scaling method, the multiplicative approach is automatically applied only when the variable is precipitation.
#' @param window character, one of none or monthly. Whether bias correction should be applied on a monthly or annual basis. Monthly is the preferred option when performing bias-correction using daily data
#' @importFrom magrittr %>%
#' @return list with SpatRaster. To explore the output run attributes(output)
#' @export
Expand Down Expand Up @@ -53,7 +54,8 @@ load_data_and_climate_change_signal <-
threshold = 0.6,
n.sessions = 6,
method = "eqm",
percentage = F) {
percentage = F,
window="monthly") {
# calculate number of chunks based on xlim and ylim
if (missing(chunk.size) | missing(season)) {
cli::cli_abort("chunk.size and season must be specified")
Expand All @@ -64,9 +66,7 @@ load_data_and_climate_change_signal <-
)
}

cli::cli_alert_warning(
"Interpolation may be required to merge the rasters. Be aware that interpolation can introduce slight discrepancies in the data, potentially affecting the consistency of results across different spatial segments."
)


country_shp = if (!is.null(country) &
!inherits(country, "sf")) {
Expand Down Expand Up @@ -170,7 +170,8 @@ load_data_and_climate_change_signal <-
threshold = threshold,
n.sessions = 1,
method = method,
percentage = percentage
percentage = percentage,
window=window
)
)

Expand All @@ -193,28 +194,36 @@ load_data_and_climate_change_signal <-
dplyr::summarise(value = median(value, na.rm = T))

merge_rasters <- function(rst_list) {
# Determine the smallest (finest) resolution among all rasters
# Determine the resolution of each raster in the list
resolutions <- sapply(rst_list, function(r)
terra::res(r))
common_res <- max(resolutions)

# Resample all rasters to the common resolution
resampled_rasters <- lapply(rst_list, function(r) {
terra::resample(r,
terra::rast(
terra::ext(r),
resolution = common_res,
crs = terra::crs(r)
),
method = "mode")
})

# Merge the resampled rasters

# Check if all rasters have the same resolution
if (length(unique(resolutions)) > 1) {
cli::cli_alert_warning("Interpolation was required to merge the rasters.")
# If resolutions differ, determine the smallest (finest) resolution among all rasters
common_res <- max(resolutions)

# Resample all rasters to the common resolution
rst_list <- lapply(rst_list, function(r) {
terra::resample(
r,
terra::rast(
terra::ext(r),
resolution = common_res,
crs = terra::crs(r)
),
method = "mode"
)
})
}

# Merge rasters
merged_raster <-
Reduce(function(x, y)
terra::merge(x, y), resampled_rasters)
terra::merge(x, y), rst_list)

#Set names from the first raster in the list
# Set names from the first raster in the list
names <- names(rst_list[[1]])
setNames(merged_raster, names)
}
Expand Down
57 changes: 34 additions & 23 deletions R/load_data_and_model_biases.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,9 @@
#' @param n.sessions numeric, number of sessions to use in parallel processing for loading the data. Default to 6. Increasing the number of sessions will not necessarily results in better performances. Leave as default unless necessary
#' @param chunk.size numeric, indicating the number of chunks. The smaller the better when working with limited RAM
#' @param overlap numeric, amount of overlap needed to create the composite. Default 0.25
#' @param method character, bias-correction method to use. One of eqm (Empirical Quantile Mapping) or qdm (Quantile Delta Mapping). Default to eqm
#' @param method character, bias-correction method to use. One of eqm (Empirical Quantile Mapping), qdm (Quantile Delta Mapping) or scaling. Default to eqm. When using the scaling method, the multiplicative approach is automatically applied only when the variable is precipitation.
#' @param cross_validation character, one of none or 3fold. Whether 3-fold cross validation should be used to avoid overfitting during bias-correction. Default to "none"
#' @param window character, one of none or monthly. Whether bias correction should be applied on a monthly or annual basis. Monthly is the preferred option when performing bias-correction using daily data
#' @importFrom magrittr %>%
#' @return list with SpatRaster. To explore the output run attributes(output)
#' @export
Expand All @@ -49,7 +51,9 @@ load_data_and_model_biases <-
bias.correction = F,
domain = NULL,
n.sessions = 6,
method = "eqm") {
method = "eqm",
cross_validation = "none",
window = "monthly") {
# calculate number of chunks based on xlim and ylim
if (missing(chunk.size) | missing(season)) {
cli::cli_abort("chunk.size and season must be specified")
Expand All @@ -61,9 +65,6 @@ load_data_and_model_biases <-
)
}

cli::cli_alert_warning(
"Interpolation may be required to merge the rasters. Be aware that interpolation can introduce slight discrepancies in the data, potentially affecting the consistency of results across different spatial segments."
)
country_shp = if (!is.null(country) &
!inherits(country, "sf")) {
suppressMessages(
Expand Down Expand Up @@ -163,7 +164,9 @@ load_data_and_model_biases <-
duration = duration,
frequency = frequency,
n.sessions = 1,
method = method
method = method,
cross_validation = cross_validation,
window = window
)
)

Expand All @@ -180,28 +183,36 @@ load_data_and_model_biases <-
rst_mbrs <- lapply(out_list, `[[`, 2)
# Merge the extracted rasters using `Reduce` and set their names
merge_rasters <- function(rst_list) {
# Determine the smallest (finest) resolution among all rasters
# Determine the resolution of each raster in the list
resolutions <- sapply(rst_list, function(r)
terra::res(r))
common_res <- max(resolutions)

# Resample all rasters to the common resolution
resampled_rasters <- lapply(rst_list, function(r) {
terra::resample(r,
terra::rast(
terra::ext(r),
resolution = common_res,
crs = terra::crs(r)
),
method = "mode")
})

# Merge the resampled rasters

# Check if all rasters have the same resolution
if (length(unique(resolutions)) > 1) {
cli::cli_alert_warning("Interpolation was required to merge the rasters.")
# If resolutions differ, determine the smallest (finest) resolution among all rasters
common_res <- max(resolutions)

# Resample all rasters to the common resolution
rst_list <- lapply(rst_list, function(r) {
terra::resample(
r,
terra::rast(
terra::ext(r),
resolution = common_res,
crs = terra::crs(r)
),
method = "mode"
)
})
}

# Merge rasters
merged_raster <-
Reduce(function(x, y)
terra::merge(x, y), resampled_rasters)
terra::merge(x, y), rst_list)

#Set names from the first raster in the list
# Set names from the first raster in the list
names <- names(rst_list[[1]])
setNames(merged_raster, names)
}
Expand Down
Loading

0 comments on commit 3bda49d

Please sign in to comment.