From a54519f71a55dce9fe6b70f86403b3428672e615 Mon Sep 17 00:00:00 2001 From: eduzenli Date: Fri, 25 Oct 2024 16:06:16 +0200 Subject: [PATCH 1/2] leave-one-year-out characteristic for daily data have been added --- R/Intbc.R | 133 ++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 98 insertions(+), 35 deletions(-) diff --git a/R/Intbc.R b/R/Intbc.R index 70043c1..d11b5e1 100644 --- a/R/Intbc.R +++ b/R/Intbc.R @@ -73,7 +73,8 @@ CST_Intbc <- function(exp, obs, exp_cor = NULL, target_grid, bc_method, int_method = NULL, points = NULL, method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", - sdate_dim = "sdate", member_dim = "member", region = NULL, ncores = NULL, ...) + sdate_dim = "sdate", member_dim = "member", time_dim = "time", + region = NULL, ncores = NULL, ...) { if (!inherits(exp,'s2dv_cube')) { stop("Parameter 'exp' must be of the class 's2dv_cube'") @@ -91,7 +92,8 @@ CST_Intbc <- function(exp, obs, exp_cor = NULL, target_grid, bc_method, int_meth source_file_exp = exp$attrs$source_files[1], source_file_obs = obs$attrs$source_files[1], method_point_interp = method_point_interp, lat_dim = lat_dim, lon_dim = lon_dim, - sdate_dim = sdate_dim, member_dim = member_dim, region = region, ncores = ncores, ...) + sdate_dim = sdate_dim, member_dim = member_dim, + time_dim = time_dim, region = region, ncores = ncores, ...) # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data obs$data <- res$obs @@ -363,45 +365,106 @@ Intbc <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo # which(names(dim(obs_ref)) == sdate_dim), 'sdate') #} - if (bc_method == 'qm' | bc_method == 'quantile_mapping') { - - res <- QuantileMapping(exp = exp_interpolated$data, obs = obs_ref, - exp_cor = exp_cor_interpolated$data, - na.rm = TRUE, memb_dim = member_dim, sdate_dim = sdate_dim, - ncores = ncores, ...) + + k_out <- loop <- 1 # the data is NOT hourly/daily/weekly, i.e., time = 1. Apply leave one data out. + daily <- FALSE # for the case at which time > 1 or for the case of forecast dwn (exp_cor) + eval.method <- "leave-one-out" + # if time > 1, apply leave-one-year-out + if ( time_dim %in% names(dim(obs_ref)) & dim(obs_ref)[time_dim] > 1 & is.null(exp_cor)) { + daily <- TRUE + eval.method <- "in-sample" + sdate_num <- as.numeric (dim(obs_ref)[sdate_dim]) + k_out <- as.numeric (dim(obs_ref)[time_dim]) ## in case loocv = TRUE and the data is daily, leave one-year data out + obs_ref <- MergeDims(obs_ref, merge_dims = c(time_dim, sdate_dim), rename_dim = sdate_dim) + exp_interpolated$data <- MergeDims(exp_interpolated$data, + merge_dims = c(time_dim, sdate_dim), rename_dim = sdate_dim) + loop <- dim(obs_ref)[sdate_dim]/k_out + } else if (!is.null(exp_cor)) { + eval.method <- "in-sample" } - else if (bc_method == 'dbc' | bc_method == 'dynamical_bias') { - # Dynamical bias correction is not yet prepared to handle hindcast-forecast data - # It will return only the hindcast downscaled - if (!is.null(exp_cor)) { - warning("For the dynamical bias correction only the hindcast downscaled will be returned.") - } - # the temporal dimension must be only one dimension called "time" - if (all(c(time_dim, sdate_dim) %in% names(dim(exp_interpolated$data)))) { - exp_interpolated$data <- Apply(exp_interpolated$data, target_dims = c(time_dim, sdate_dim), - fun = as.vector, output_dims = "time", ncores = ncores)$output1 - } - if (all(c(time_dim, sdate_dim) %in% names(dim(obs_ref)))) { - obs_ref <- Apply(obs_ref, target_dims = c(time_dim, sdate_dim), fun = as.vector, - output_dims = "time", ncores = ncores)$output1 + + sdate_num <- as.numeric (dim(obs_ref)[sdate_dim]) + res <- list() + + for (yy in 1:loop) { # with cv purposes. leave-one-year-out + + if (loop == 1) { + window <- -(1:sdate_num) + EXP_COR <- NULL + if (!is.null(exp_cor)) { + EXP_COR <- exp_cor_interpolated$data + } + } else { + window <- ((yy-1) * k_out+1):((yy-1) * k_out + k_out) + EXP_COR <- ClimProjDiags::Subset(exp_interpolated$data, along = sdate_dim, + indices = (1:sdate_num)[window], drop = "selected") } - # REMEMBER to add na.rm = T in colMeans in .proxiesattractor - res <- DynBiasCorrection(exp = exp_interpolated$data, obs = obs_ref, ncores = ncores, ...) - } else { - if (dim(exp_interpolated$data)[member_dim] == 1) { - stop('Calibration must not be used with only one ensemble member.') + + OBS_REF <- ClimProjDiags::Subset(obs_ref, along = sdate_dim, + indices = (1:sdate_num)[-window], drop = "selected") + EXP <- ClimProjDiags::Subset(exp_interpolated$data, along = sdate_dim, + indices = (1:sdate_num)[-window], drop = "selected") + + if (bc_method == 'qm' | bc_method == 'quantile_mapping') { + dum <- QuantileMapping(exp = EXP, obs = OBS_REF, + exp_cor = EXP_COR, + na.rm = TRUE, memb_dim = member_dim, sdate_dim = sdate_dim, + ncores = ncores, eval.method = eval.method, ...) + + res[[yy]] <- dum + rm(dum) } - if (dim(obs_ref)[sdate_dim] == 1) { - warning('Simple Bias Correction should not be used with only one observation. Returning NA.') + else if (bc_method == 'dbc' | bc_method == 'dynamical_bias') { + # Dynamical bias correction is not yet prepared to handle hindcast-forecast data + # It will return only the hindcast downscaled + # Dynamical bias correction is not yet prepared to handle leave-one-year-out for daily case + # It will properly work only with the monthly or coarser temporal scale data. + if (!is.null(exp_cor)) { + warning("For the dynamical bias correction only the hindcast downscaled will be returned.") + } + # the temporal dimension must be only one dimension called "time" + if (all(c(time_dim, sdate_dim) %in% names(dim(exp_interpolated$data)))) { + exp_interpolated$data <- Apply(exp_interpolated$data, target_dims = c(time_dim, sdate_dim), + fun = as.vector, output_dims = "time", ncores = ncores)$output1 + } + if (all(c(time_dim, sdate_dim) %in% names(dim(obs_ref)))) { + obs_ref <- Apply(obs_ref, target_dims = c(time_dim, sdate_dim), fun = as.vector, + output_dims = "time", ncores = ncores)$output1 + } + # REMEMBER to add na.rm = T in colMeans in .proxiesattractor + dum <- DynBiasCorrection(exp = EXP, obs = OBS_REF, ncores = ncores, ...) + res[[yy]] <- dum + rm(dum) + } else { + if (dim(exp_interpolated$data)[member_dim] == 1) { + stop('Calibration must not be used with only one ensemble member.') + } + if (dim(obs_ref)[sdate_dim] == 1) { + warning('Simple Bias Correction should not be used with only one observation. Returning NA.') + } + dum <- Calibration(exp = EXP, obs = OBS_REF, + exp_cor = EXP_COR, + memb_dim = member_dim, sdate_dim = sdate_dim, ncores = ncores, + cal.method = bc_method, + eval.method = eval.method) + res[[yy]] <- dum + rm(dum) } - res <- Calibration(exp = exp_interpolated$data, obs = obs_ref, - exp_cor = exp_cor_interpolated$data, - memb_dim = member_dim, sdate_dim = sdate_dim, ncores = ncores, - cal.method = bc_method) } - + + require (abind) + dimnames <- names(dim(res[[1]])) + res <- do.call(abind, c(res, along = which(names(dim(res[[1]])) == "syear"))) + names(dim(res)) <- dimnames + + if (daily) { + res <- SplitDim(res, split_dim = sdate_dim, new_dim_name = time_dim, + indices = rep(1:k_out, sdate_num/k_out)) + obs_ref <- SplitDim(obs_ref, split_dim = sdate_dim, new_dim_name = time_dim, + indices = rep(1:k_out, sdate_num/k_out)) + } + # Return a list of three elements res <- list(data = res, obs = obs_ref, lon = exp_interpolated$lon, lat = exp_interpolated$lat) - return(res) } -- GitLab From 38a962b67ab11c771dd5d1d00d437026b86d88b1 Mon Sep 17 00:00:00 2001 From: eduzenli Date: Fri, 25 Oct 2024 16:06:43 +0200 Subject: [PATCH 2/2] fine tuning --- R/Intlr.R | 2 +- R/LogisticReg.R | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/R/Intlr.R b/R/Intlr.R index 8b5dabe..c0f3b9c 100644 --- a/R/Intlr.R +++ b/R/Intlr.R @@ -861,7 +861,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo "grd" ), fun = function(x, y) { a <- x[, y[1], y[2], , drop = FALSE] - a <- Subset (a, along = c(lat_dim, lon_dim), + a <- ClimProjDiags::Subset (a, along = c(lat_dim, lon_dim), indices = list (1,1), drop = "selected") # drop lat lon dim coming from coar data return(a) diff --git a/R/LogisticReg.R b/R/LogisticReg.R index c0d914a..2a6259b 100644 --- a/R/LogisticReg.R +++ b/R/LogisticReg.R @@ -484,8 +484,6 @@ LogisticReg <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_ref <- obs } - # convert observations to categorical predictands - k_out <- 1 ## in case loocv = TRUE and the data is NOT daily, leave one data out. daily <- FALSE if ( time_dim %in% names(dim(obs_ref)) & dim(obs_ref)[time_dim] > 1) { @@ -500,6 +498,7 @@ LogisticReg <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, } } + # convert observations to categorical predictands obs_cat <- Apply(obs_ref, target_dims = sdate_dim, function(x) { if (!any(!is.na(x))) { rep(NA,length(x)) -- GitLab