From 62479ca38e4c5426681972661da0e1052821938e Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Fri, 17 Nov 2023 17:03:49 +0100 Subject: [PATCH 01/17] adapted to hcst-fcst --- R/Intbc.R | 44 +++++++++++++++++++++++++++++++++++++++----- 1 file changed, 39 insertions(+), 5 deletions(-) diff --git a/R/Intbc.R b/R/Intbc.R index dc5d050..d6b8dce 100644 --- a/R/Intbc.R +++ b/R/Intbc.R @@ -187,7 +187,7 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point #'res <- Intbc(exp = exp, obs = obs, exp_lats = exp_lats, exp_lons = exp_lons, obs_lats = obs_lats, #'obs_lons = obs_lons, target_grid = 'r1280x640', bc_method = 'simple_bias', int_method = 'conservative') #'@export -Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, bc_method, int_method = NULL, +Intbc <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, bc_method, int_method = NULL, points = NULL, method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", member_dim = "member", source_file_exp = NULL, source_file_obs = NULL, region = NULL, ncores = NULL, ...) { @@ -238,6 +238,27 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, "'crps_min', 'rpc-based'. The abbreviations 'dbc','qm' can also be used.") } + if (!is.null(exp_cor)) { + if (is.na(match(lon_dim, names(dim(exp_cor))))) { + stop("Missing longitude dimension in 'exp_cor', or does not match the parameter ", + "'lon_dim'") + } + + if (is.na(match(lat_dim, names(dim(exp_cor))))) { + stop("Missing latitude dimension in 'exp_cor', or does not match the parameter ", + "'lat_dim'") + } + + if (is.na(match(sdate_dim, names(dim(exp_cor))))) { + stop("Missing start date dimension in 'exp_cor', or does not match the parameter ", + "'sdate_dim'") + } + + if (is.na(match(member_dim, names(dim(exp_cor))))) { + stop("Missing member dimension in 'exp_cor', or does not match the parameter 'member_dim'") + } + } + # When observations are pointwise if (!is.null(points) & !is.na(match("location", names(dim(obs))))) { point_obs <- T @@ -276,6 +297,14 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, method_remap = int_method, points = points, source_file = source_file_exp, lat_dim = lat_dim, lon_dim = lon_dim, method_point_interp = method_point_interp, region = region, ncores = ncores) + + if (!is.null(exp_cor)) { + exp_cor_interpolated <- Interpolation(exp = exp_cor, lats = exp_lats, lons = exp_lons, target_grid = target_grid, + method_remap = int_method, points = points, source_file = source_file_exp, + lat_dim = lat_dim, lon_dim = lon_dim, method_point_interp = method_point_interp, + region = region, ncores = ncores) + + } # If after interpolating 'exp' data the coordinates do not match, the obs data is interpolated to # the same grid to force the matching @@ -306,10 +335,15 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, if (bc_method == 'qm' | bc_method == 'quantile_mapping') { - res <- QuantileMapping(exp = exp_interpolated$data, obs = obs_ref, na.rm = TRUE, - memb_dim = member_dim, sdate_dim = sdate_dim, ncores = ncores, ...) + 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, ...) } 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), @@ -328,8 +362,8 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, if (dim(obs_ref)[sdate_dim] == 1) { warning('Simple Bias Correction should not be used with only one observation. Returning NA.') } - res <- Calibration(exp = exp_interpolated$data, obs = obs_ref, memb_dim = member_dim, - sdate_dim = sdate_dim, ncores = ncores, cal.method = bc_method) + 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) } # Return a list of three elements -- GitLab From 89a7953cf122afa243ea47683638989d66ea7590 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Wed, 29 Nov 2023 17:23:35 +0100 Subject: [PATCH 02/17] adapted to hindcast-forecast --- R/Intbc.R | 43 ++++++--- R/Intlr.R | 275 +++++++++++++++++++++++++++++++++++++++++++----------- 2 files changed, 252 insertions(+), 66 deletions(-) diff --git a/R/Intbc.R b/R/Intbc.R index d6b8dce..58ae18f 100644 --- a/R/Intbc.R +++ b/R/Intbc.R @@ -20,6 +20,10 @@ #'@param obs an 's2dv object' containing the observational field. The object #'must have, at least, the dimensions latitude, longitude and start date. The object is #'expected to be already subset for the desired region. +#'@param exp_cor an optional 's2dv_cube' object with named dimensions containing the seasonal +#'forecast experiment data. If the forecast is provided, it will be downscaled using the +#'hindcast and observations; if not provided, the hindcast will be downscaled instead. The +#'default value is NULL. #'@param target_grid a character vector indicating the target grid to be passed to CDO. #'It must be a grid recognised by CDO or a NetCDF file. #'@param bc_method a character vector indicating the bias adjustment method to be applied after @@ -67,7 +71,7 @@ #'res <- CST_Intbc(exp = exp, obs = obs, target_grid = 'r1280x640', bc_method = 'simple_bias', int_method = 'conservative') #'@export -CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, points = NULL, +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, ...) { @@ -79,25 +83,36 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point stop("Parameter 'obs' must be of the class 's2dv_cube'") } - res <- Intbc(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], exp_lons = exp$coords[[lon_dim]], - obs_lats = obs$coords[[lat_dim]], obs_lons = obs$coords[[lon_dim]], target_grid = target_grid, + res <- Intbc(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, exp_lats = exp$attrs[[lat_dim]], + exp_lons = exp$attrs[[lon_dim]], obs_lats = obs$attrs[[lat_dim]], + obs_lons = obs$attrs[[lon_dim]], target_grid = target_grid, int_method = int_method, bc_method = bc_method, points = points, 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, ...) # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data - exp$data <- res$data - exp$dims <- dim(exp$data) - exp$coords[[lon_dim]] <- res$lon - exp$coords[[lat_dim]] <- res$lat - obs$data <- res$obs obs$dims <- dim(obs$data) - obs$coords[[lon_dim]] <- res$lon - obs$coords[[lat_dim]] <- res$lat + obs$attrs[[lon_dim]] <- res$lon + obs$attrs[[lat_dim]] <- res$lat + + if (is.null(exp_cor)) { + exp$data <- res$data + exp$dims <- dim(exp$data) + exp$attrs[[lon_dim]] <- res$lon + exp$attrs[[lat_dim]] <- res$lat + + res_s2dv <- list(exp = exp, obs = obs) + } else { + exp_cor$data <- res$data + exp_cor$dims <- dim(exp_cor$data) + exp_cor$attrs[[lon_dim]] <- res$lon + exp_cor$attrs[[lat_dim]] <- res$lat + + res_s2dv <- list(exp = exp_cor, obs = obs) + } - res_s2dv <- list(exp = exp, obs = obs) return(res_s2dv) } @@ -123,7 +138,11 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point #''region'. #'@param obs an array with named dimensions containing the observational field. The object #'must have, at least, the dimensions latitude, longitude and start date. The object is -#'expected to be already subset for the desired region. +#'expected to be already subset for the desired region. +#'@param exp_cor an optional array with named dimensions containing the seasonal forecast +#'experiment data. If the forecast is provided, it will be downscaled using the hindcast and +#'observations; if not provided, the hindcast will be downscaled instead. The default value +#'is NULL. #'@param exp_lats a numeric vector containing the latitude values in 'exp'. Latitudes must #'range from -90 to 90. #'@param exp_lons a numeric vector containing the longitude values in 'exp'. Longitudes diff --git a/R/Intlr.R b/R/Intlr.R index 36a7f11..9ce7bea 100644 --- a/R/Intlr.R +++ b/R/Intlr.R @@ -22,6 +22,10 @@ #'@param obs an 's2dv object' containing the observational field. The object #'must have, at least, the dimensions latitude, longitude and start date. The object is #'expected to be already subset for the desired region. +#'@param exp_cor an optional 's2dv_cube' object with named dimensions containing the seasonal +#'forecast experiment data. If the forecast is provided, it will be downscaled using the +#'hindcast and observations; if not provided, the hindcast will be downscaled instead. The +#'default value is NULL. #'@param lr_method a character vector indicating the linear regression method to be applied #'after the interpolation. Accepted methods are 'basic', 'large-scale' and '4nn'. The 'basic' #'method fits a linear regression using high resolution observations as predictands and the @@ -86,7 +90,7 @@ #'obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons) #'res <- CST_Intlr(exp = exp, obs = obs, target_grid = 'r1280x640', lr_method = 'basic', int_method = 'conservative') #'@export -CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, int_method = NULL, +CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, points = NULL, int_method = NULL, method_point_interp = NULL, predictors = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", member_dim = "member", large_scale_predictor_dimname = 'vars', loocv = TRUE, region = NULL, ncores = NULL) { @@ -98,28 +102,54 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in if (!inherits(obs,'s2dv_cube')) { stop("Parameter 'obs' must be of the class 's2dv_cube'") } - - res <- Intlr(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], exp_lons = exp$coords[[lon_dim]], - obs_lats = obs$coords[[lat_dim]], obs_lons = obs$coords[[lon_dim]], points = points, - source_file_exp = exp$attrs$source_files[1], source_file_obs = obs$attrs$source_files[1], + + if (identical(lr_method, 'basic') | identical(lr_method, '4nn')) { + if (!inherits(exp_cor,'s2dv_cube')) { + stop("Parameter 'exp_cor' must be of the class 's2dv_cube'") + } + + exp_cor_aux <- exp_cor$data + # when large-scale is selected, the forecast object does not have to be an s2dv_cube object + } else if (identical(lr_method, 'large-scale')) { + exp_cor_aux <- exp_cor + } + + res <- Intlr(exp = exp$data, obs = obs$data, exp_cor = exp_cor_aux, exp_lats = exp$attrs[[lat_dim]], + exp_lons = exp$attrs[[lon_dim]], obs_lats = obs$attrs[[lat_dim]], obs_lons = obs$attrs[[lon_dim]], + points = points, source_file_exp = exp$attrs$source_files[1], source_file_obs = obs$attrs$source_files[1], target_grid = target_grid, lr_method = lr_method, int_method = int_method, method_point_interp = method_point_interp, predictors = predictors, lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, time_dim = time_dim, member_dim = member_dim, large_scale_predictor_dimname = large_scale_predictor_dimname, loocv = loocv, region = region, ncores = ncores) - + # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data - exp$data <- res$data - exp$dims <- dim(exp$data) - exp$coords[[lon_dim]] <- res$lon - exp$coords[[lat_dim]] <- res$lat - obs$data <- res$obs obs$dims <- dim(obs$data) - obs$coords[[lon_dim]] <- res$lon - obs$coords[[lat_dim]] <- res$lat - - res_s2dv <- list(exp = exp, obs = obs) + obs$attrs[[lon_dim]] <- res$lon + obs$attrs[[lat_dim]] <- res$lat + + if (is.null(exp_cor)) { + exp$data <- res$data + exp$dims <- dim(exp$data) + exp$attrs[[lon_dim]] <- res$lon + exp$attrs[[lat_dim]] <- res$lat + + res_s2dv <- list(exp = exp, obs = obs) + } else { + if (identical(lr_method, 'basic') | identical(lr_method, '4nn')) { + exp_cor$data <- res$data + exp_cor$dims <- dim(exp_cor$data) + exp_cor$attrs[[lon_dim]] <- res$lon + exp_cor$attrs[[lat_dim]] <- res$lat + # when large-scale is selected, the forecast object does not have to be an s2dv_cube object + } else if (identical(lr_method, 'large-scale')) { + exp_cor <- suppressWarnings(s2dv_cube(res$data, lat = res$lat, lon = res$lon)) + } + + res_s2dv <- list(exp = exp_cor, obs = obs) + } + return(res_s2dv) } @@ -147,6 +177,10 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in #'@param obs an array with named dimensions containing the observational field. The object #'must have, at least, the dimensions latitude, longitude and start date. The object is #'expected to be already subset for the desired region. +#'@param exp_cor an optional array with named dimensions containing the seasonal forecast +#'experiment data. If the forecast is provided, it will be downscaled using the hindcast and +#'observations; if not provided, the hindcast will be downscaled instead. The default value +#'is NULL. #'@param exp_lats a numeric vector containing the latitude values in 'exp'. Latitudes must #'range from -90 to 90. #'@param exp_lons a numeric vector containing the longitude values in 'exp'. Longitudes @@ -222,7 +256,7 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in #'res <- Intlr(exp = exp, obs = obs, exp_lats = exp_lats, exp_lons = exp_lons, obs_lats = obs_lats, #'obs_lons = obs_lons, target_grid = 'r1280x640', lr_method = 'basic', int_method = 'conservative') #'@export -Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, target_grid = NULL, points = NULL, +Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, target_grid = NULL, points = NULL, int_method = NULL, method_point_interp = NULL, source_file_exp = NULL, source_file_obs = NULL, predictors = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", member_dim = "member", region = NULL, large_scale_predictor_dimname = 'vars', loocv = TRUE, @@ -274,6 +308,35 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t "'sdate_dim'") } + if (!is.null(exp_cor)) { + + if (is.na(match(sdate_dim, names(dim(exp_cor))))) { + stop("Missing start date dimension in 'exp_cor', or does not match the parameter ", + "'sdate_dim'") + } + + if (is.na(match(member_dim, names(dim(exp_cor))))) { + stop("Missing member dimension in 'exp_cor', or does not match the parameter 'member_dim'") + } + + if (loocv) { # loocv equal to false to train with the whole hindcast and predict with the forecast + loocv <- FALSE + } + + if (is.null(predictors)) { + + if (is.na(match(lon_dim, names(dim(exp_cor))))) { + stop("Missing longitude dimension in 'exp_cor', or does not match the parameter ", + "'lon_dim'") + } + + if (is.na(match(lat_dim, names(dim(exp_cor))))) { + stop("Missing latitude dimension in 'exp_cor', or does not match the parameter ", + "'lat_dim'") + } + } + } + # When observations are pointwise if (!is.null(points) & !is.na(match("location", names(dim(obs))))) { point_obs <- T @@ -295,10 +358,6 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t "through the parameter 'method_point_interp'.") } - # sdate must be the time dimension in the input data - stopifnot(sdate_dim %in% names(dim(exp))) - stopifnot(sdate_dim %in% names(dim(obs))) - # the code is not yet prepared to handle members in the observations restore_ens <- FALSE if (member_dim %in% names(dim(obs))) { @@ -321,6 +380,17 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t stopifnot(sdate_dim %in% names(dim(predictors))) stopifnot(dim(predictors)[sdate_dim] == dim(exp)[sdate_dim]) } + # forecasts for the large scale predictors + if (!is.null(exp_cor)) { + if (is.na(match(large_scale_predictor_dimname, names(dim(exp_cor))))) { + stop("Missing large scale predictor dimension in 'exp_cor', or does not match the parameter ", + "'large_scale_predictor_dimname'") + } + if (!identical(dim(exp_cor)[names(dim(exp_cor)) == large_scale_predictor_dimname], + dim(predictors)[names(dim(predictors)) == large_scale_predictor_dimname])) { + stop("Large scale predictor dimension in exp_cor and predictors must be identical.") + } + } } ## ncores if (!is.null(ncores)) { @@ -351,6 +421,13 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t points = points, method_point_interp = method_point_interp, source_file = source_file_exp, lat_dim = lat_dim, lon_dim = lon_dim, method_remap = int_method, region = region, ncores = ncores) + + if (!is.null(exp_cor) & is.null(predictors)) { + exp_cor_interpolated <- Interpolation(exp = exp_cor, lats = exp_lats, lons = exp_lons, target_grid = target_grid, + points = points, method_point_interp = method_point_interp, + source_file = source_file_exp, lat_dim = lat_dim, lon_dim = lon_dim, + method_remap = int_method, region = region, ncores = ncores) + } # If after interpolating 'exp' data the coordinates do not match, the obs data is interpolated to # the same grid to force the matching @@ -383,6 +460,17 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t target_dims_predictor <- sdate_dim target_dims_predictand <- sdate_dim + + if (!is.null(exp_cor)) { + aux_dim <- NULL + forecast <- exp_cor_interpolated$data + target_dims_predictor <- c(sdate_dim, member_dim) + target_dims_forecast <- c(sdate_dim, member_dim) + } else { + forecast <- NULL + target_dims_forecast <- NULL + } + } # (Multi) linear regression with large-scale predictors @@ -398,6 +486,13 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t target_dims_predictor <- c(sdate_dim, large_scale_predictor_dimname) target_dims_predictand <- sdate_dim + + if (!is.null(exp_cor)) { + aux_dim <- large_scale_predictor_dimname + forecast <- exp_cor + target_dims_predictor <- c(sdate_dim, large_scale_predictor_dimname, member_dim) + target_dims_forecast <- c(sdate_dim, large_scale_predictor_dimname, member_dim) + } } # Multi-linear regression with the four nearest neighbours @@ -406,8 +501,14 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t else if (lr_method == '4nn') { predictor <- .find_nn(coar = exp, lats_hres = obs_lats, lons_hres = obs_lons, lats_coar = exp_lats, - lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, nn = 4, ncores = ncores) - + lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, nn = 4, ncores = ncores) + + if (!is.null(exp_cor)) { + aux_dim <- 'nn' + forecast <- .find_nn(coar = exp_cor, lats_hres = obs_lats, lons_hres = obs_lons, lats_coar = exp_lats, + lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, nn = 4, ncores = ncores) + } + if (is.null(points) | ("location" %in% names(dim(obs)))) { if (!is.null(target_grid)) { warning("Interpolating to the 'obs' grid") @@ -433,27 +534,34 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t predictor <- predictor$data predictand <- predictand$data } - - target_dims_predictor <- c(sdate_dim,'nn') - target_dims_predictand <- sdate_dim + + target_dims_predictand <- sdate_dim + + if (!is.null(exp_cor)) { + target_dims_predictor <- c(sdate_dim,'nn', member_dim) + target_dims_forecast <- c(sdate_dim,'nn', member_dim) + } else { + target_dims_predictor <- c(sdate_dim,'nn') + } } else { stop(paste0(lr_method, " method is not implemented yet")) } - - print(paste0('dim predictor',dim(predictor))) - print(paste0('dim predictand',dim(predictand))) - print(dim(list(predictor[1]))) - # Apply the linear regressions - - - res <- Apply(list(predictor, predictand), target_dims = list(target_dims_predictor, target_dims_predictand), - fun = .intlr, loocv = loocv, ncores = ncores)$output1 - - names(dim(res))[1] <- sdate_dim - # names(dim(res))[which(names(dim(res)) == '')] + # Apply the linear regressions + # case hindcast - forecast + if (!is.null(exp_cor)) { + res <- Apply(list(predictor, predictand, forecast), + target_dims = list(target_dims_predictor, target_dims_predictand, target_dims_forecast), + fun = .intlr, loocv = loocv, aux_dim = aux_dim, ncores = ncores)$output1 + } + # case hindcast only + else { + res <- Apply(list(predictor, predictand), target_dims = list(target_dims_predictor, target_dims_predictand), + fun = .intlr, loocv = loocv, ncores = ncores)$output1 + names(dim(res))[1] <- sdate_dim + } # restore ensemble dimension in observations if it existed originally if (restore_ens) { @@ -463,30 +571,56 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t # Return a list of three elements res <- list(data = res, obs = predictand, lon = lons, lat = lats) + # for testing 4nn + #x <- predictor[10,10,1,1,1,,,] + #x <- Reorder(x,c(2,3,1)) + #y <- predictand[1,1,1,10,,10] + #f <- forecast[10,10,1,1,1,,,] + #f <- Reorder(f,c(2,1)) + #f <- InsertDim(f,1,1,'sdate') + + # large-scale + #x <- Reorder(predictor,c(1,3,2)) + #y <- predictand[1,1,1,10,,10] + #f <- Reorder(forecast,c(1,3,2)) + return(res) } #----------------------------------- # Atomic function to generate and apply the linear regressions #----------------------------------- -.intlr <- function(x, y, loocv) { - - tmp_df <- data.frame(x = x, y = y) +.intlr <- function(x, y, f = NULL, loocv, aux_dim = NULL) { + + if (!is.null(f)) { + if (!is.null(aux_dim)) { + tmp_df <- data.frame(x = adply(x,.margins = 3, .id = NULL), y = y) + } else { + tmp_df <- data.frame(x = as.vector(x), y = y) + } + } else { + tmp_df <- data.frame(x = x, y = y) + } # if the data is all NA, force return return NA if (all(is.na(tmp_df)) | (sum(apply(tmp_df, 2, function(x) !all(is.na(x)))) == 1)) { - - n <- nrow(tmp_df) - res <- rep(NA, n) - + if (is.null(f)) { + res <- rep(NA, nrow(tmp_df)) + } else { + if (!is.null(aux_dim)) { + res <- array(NA, dim(f)[names(dim(f)) != aux_dim]) + } else { + res <- array(NA, dim(f)) + } + } } else { # training lm1 <- .train_lm(df = tmp_df, loocv = loocv) - + # prediction - res <- .pred_lm(lm1 = lm1, df = tmp_df, loocv = loocv) + res <- .pred_lm(lm1 = lm1, df = tmp_df, f = f, loocv = loocv, aux_dim = aux_dim) } - - return(res) + + return(res) } #----------------------------------- @@ -517,7 +651,7 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t #----------------------------------- # Function to apply the linear regressions. #----------------------------------- -.pred_lm <- function(df, lm1, loocv) { +.pred_lm <- function(df, lm1, f, loocv, aux_dim) { if (loocv) { pred_vals <- sapply(1:nrow(df), function(j) { @@ -527,13 +661,46 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t return(predict(lm1[[j]], df[j,])) }}) } else { - if (!is.na(lm1)) { - pred_vals_ls <- lapply(lm1, predict, data = df) - pred_vals <- unlist(pred_vals_ls) + if (!is.na(lm1)) { + # case to downscale hindcasts + if (is.null(f)) { + pred_vals_ls <- lapply(lm1, predict, data = df) + pred_vals <- unlist(pred_vals_ls) + # case to downscale forecasts + } else { + if (!is.null(aux_dim)) { + fcst_df <- adply(f,.margins = 3, .id = NULL) + # 4nn + if (identical(aux_dim, 'nn')) { + colnames(fcst_df) <- paste0('x.',1:4) + pred_vals_ls <- lapply(lm1, predict, newdata = fcst_df) + pred_vals <- unlist(pred_vals_ls) + dim(pred_vals) <- dim(f)[names(dim(f)) != aux_dim] + } else { + pred_vals_ls <- lapply(lm1, predict, newdata = data.frame(x = fcst_df)) + pred_vals <- unlist(pred_vals_ls) + dim(pred_vals) <- dim(f)[names(dim(f)) != aux_dim] + } + # basic + } else { + pred_vals_ls <- lapply(lm1, predict, newdata = data.frame(x = as.vector(f))) + pred_vals <- unlist(pred_vals_ls) + dim(pred_vals) <- dim(f) + } + } } else { - pred_vals <- rep(NA, nrow(df)) + if (is.null(f)) { + pred_vals <- rep(NA, nrow(df)) + } else { + if (!is.null(aux_dim)) { + pred_vals <- array(NA, dim(f)[names(dim(f)) != aux_dim]) + } else { + pred_vals <- array(NA, dim(f)) + } + } } - } + } + return(pred_vals) } -- GitLab From 72405abd8002fae8a8061b25bf98f43cdb320853 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Wed, 29 Nov 2023 17:24:11 +0100 Subject: [PATCH 03/17] corrected coords by attrs for lat lon values --- R/Analogs.R | 14 +++++++------- R/Interpolation.R | 6 +++--- R/LogisticReg.R | 14 +++++++------- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/R/Analogs.R b/R/Analogs.R index 1cbea15..9c444b7 100644 --- a/R/Analogs.R +++ b/R/Analogs.R @@ -99,9 +99,9 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo stop("Parameter 'obs' must be of the class 's2dv_cube'") } - res <- Analogs(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], - exp_lons = exp$coords[[lon_dim]], obs_lats = obs$coords[[lat_dim]], - obs_lons = obs$coords[[lon_dim]], grid_exp = grid_exp, nanalogs = nanalogs, + res <- Analogs(exp = exp$data, obs = obs$data, exp_lats = exp$attrs[[lat_dim]], + exp_lons = exp$attrs[[lon_dim]], obs_lats = obs$attrs[[lat_dim]], + obs_lons = obs$attrs[[lon_dim]], grid_exp = grid_exp, nanalogs = nanalogs, fun_analog = fun_analog, lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, time_dim = time_dim, member_dim = member_dim, region = region, return_indices = return_indices, loocv_window = loocv_window, @@ -110,13 +110,13 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data exp$data <- res$data exp$dims <- dim(exp$data) - exp$coords[[lon_dim]] <- res$lon - exp$coords[[lat_dim]] <- res$lat + exp$attrs[[lon_dim]] <- res$lon + exp$attrs[[lat_dim]] <- res$lat obs$data <- res$obs obs$dims <- dim(obs$data) - obs$coords[[lon_dim]] <- res$lon - obs$coords[[lat_dim]] <- res$lat + obs$attrs[[lon_dim]] <- res$lon + obs$attrs[[lat_dim]] <- res$lat res_s2dv <- list(exp = exp, obs = obs) return(res_s2dv) diff --git a/R/Interpolation.R b/R/Interpolation.R index ed79f4f..906ba26 100644 --- a/R/Interpolation.R +++ b/R/Interpolation.R @@ -68,7 +68,7 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr stop("The name of the latitude/longitude dimensions in 'exp$data' must match the parametres 'lat_dim' and 'lon_dim'") } - res <- Interpolation(exp = exp$data, lats = exp$coords[[lat_dim]], lons = exp$coords[[lon_dim]], + res <- Interpolation(exp = exp$data, lats = exp$attrs[[lat_dim]], lons = exp$attrs[[lon_dim]], source_file = exp$attrs$source_files[1], points = points, method_remap = method_remap, target_grid = target_grid, lat_dim = lat_dim, lon_dim = lon_dim, region = region, method_point_interp = method_point_interp, ncores = ncores) @@ -76,8 +76,8 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data exp$data <- res$data exp$dims <- dim(exp$data) - exp$coords[[lon_dim]] <- res$lon - exp$coords[[lat_dim]] <- res$lat + exp$attrs[[lon_dim]] <- res$lon + exp$attrs[[lat_dim]] <- res$lat res_s2dv <- list(exp = exp, obs = NULL) return(res_s2dv) diff --git a/R/LogisticReg.R b/R/LogisticReg.R index a85a1b3..60aa573 100644 --- a/R/LogisticReg.R +++ b/R/LogisticReg.R @@ -99,9 +99,9 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me stop("Parameter 'obs' must be of the class 's2dv_cube'") } - res <- LogisticReg(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], - exp_lons = exp$coords[[lon_dim]], obs_lats = obs$coords[[lat_dim]], - obs_lons = obs$coords[[lon_dim]], target_grid = target_grid, + res <- LogisticReg(exp = exp$data, obs = obs$data, exp_lats = exp$attrs[[lat_dim]], + exp_lons = exp$attrs[[lon_dim]], obs_lats = obs$attrs[[lat_dim]], + obs_lons = obs$attrs[[lon_dim]], target_grid = target_grid, probs_cat = probs_cat, return_most_likely_cat = return_most_likely_cat, int_method = int_method, log_reg_method = log_reg_method, points = points, method_point_interp = method_point_interp, lat_dim = lat_dim, @@ -112,13 +112,13 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data exp$data <- res$data exp$dims <- dim(exp$data) - exp$coords[[lon_dim]] <- res$lon - exp$coords[[lat_dim]] <- res$lat + exp$attrs[[lon_dim]] <- res$lon + exp$attrs[[lat_dim]] <- res$lat obs$data <- res$obs obs$dims <- dim(obs$data) - obs$coords[[lon_dim]] <- res$lon - obs$coords[[lat_dim]] <- res$lat + obs$attrs[[lon_dim]] <- res$lon + obs$attrs[[lat_dim]] <- res$lat res_s2dv <- list(exp = exp, obs = obs) return(res_s2dv) -- GitLab From 644756c78cd1a33d7408f7b329e29fd792fc6ec6 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Wed, 13 Dec 2023 15:49:00 +0100 Subject: [PATCH 04/17] added hcst-fcst for Interpolation, Intbc, Intlr and LogisticReg and changed coords by attrs --- R/Analogs.R | 14 +-- R/Intbc.R | 20 ++-- R/Interpolation.R | 8 +- R/Intlr.R | 18 +-- R/LogisticReg.R | 190 +++++++++++++++++++++++-------- inst/examples/analogs.R | 6 +- inst/examples/interpolation-bc.R | 45 +++++--- inst/examples/interpolation-lr.R | 70 +++++++++--- inst/examples/interpolation.R | 7 +- inst/examples/logistic-reg.R | 42 +++++-- 10 files changed, 294 insertions(+), 126 deletions(-) diff --git a/R/Analogs.R b/R/Analogs.R index 9c444b7..1cbea15 100644 --- a/R/Analogs.R +++ b/R/Analogs.R @@ -99,9 +99,9 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo stop("Parameter 'obs' must be of the class 's2dv_cube'") } - res <- Analogs(exp = exp$data, obs = obs$data, exp_lats = exp$attrs[[lat_dim]], - exp_lons = exp$attrs[[lon_dim]], obs_lats = obs$attrs[[lat_dim]], - obs_lons = obs$attrs[[lon_dim]], grid_exp = grid_exp, nanalogs = nanalogs, + res <- Analogs(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], + exp_lons = exp$coords[[lon_dim]], obs_lats = obs$coords[[lat_dim]], + obs_lons = obs$coords[[lon_dim]], grid_exp = grid_exp, nanalogs = nanalogs, fun_analog = fun_analog, lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, time_dim = time_dim, member_dim = member_dim, region = region, return_indices = return_indices, loocv_window = loocv_window, @@ -110,13 +110,13 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data exp$data <- res$data exp$dims <- dim(exp$data) - exp$attrs[[lon_dim]] <- res$lon - exp$attrs[[lat_dim]] <- res$lat + exp$coords[[lon_dim]] <- res$lon + exp$coords[[lat_dim]] <- res$lat obs$data <- res$obs obs$dims <- dim(obs$data) - obs$attrs[[lon_dim]] <- res$lon - obs$attrs[[lat_dim]] <- res$lat + obs$coords[[lon_dim]] <- res$lon + obs$coords[[lat_dim]] <- res$lat res_s2dv <- list(exp = exp, obs = obs) return(res_s2dv) diff --git a/R/Intbc.R b/R/Intbc.R index 58ae18f..e78308b 100644 --- a/R/Intbc.R +++ b/R/Intbc.R @@ -83,32 +83,32 @@ CST_Intbc <- function(exp, obs, exp_cor = NULL, target_grid, bc_method, int_meth stop("Parameter 'obs' must be of the class 's2dv_cube'") } - res <- Intbc(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, exp_lats = exp$attrs[[lat_dim]], - exp_lons = exp$attrs[[lon_dim]], obs_lats = obs$attrs[[lat_dim]], - obs_lons = obs$attrs[[lon_dim]], target_grid = target_grid, + res <- Intbc(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, exp_lats = exp$coords[[lat_dim]], + exp_lons = exp$coords[[lon_dim]], obs_lats = obs$coords[[lat_dim]], + obs_lons = obs$coords[[lon_dim]], target_grid = target_grid, int_method = int_method, bc_method = bc_method, points = points, - source_file_exp = exp$attrs$source_files[1], source_file_obs = obs$attrs$source_files[1], + source_file_exp = exp$coords$source_files[1], source_file_obs = obs$coords$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, ...) # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data obs$data <- res$obs obs$dims <- dim(obs$data) - obs$attrs[[lon_dim]] <- res$lon - obs$attrs[[lat_dim]] <- res$lat + obs$coords[[lon_dim]] <- res$lon + obs$coords[[lat_dim]] <- res$lat if (is.null(exp_cor)) { exp$data <- res$data exp$dims <- dim(exp$data) - exp$attrs[[lon_dim]] <- res$lon - exp$attrs[[lat_dim]] <- res$lat + exp$coords[[lon_dim]] <- res$lon + exp$coords[[lat_dim]] <- res$lat res_s2dv <- list(exp = exp, obs = obs) } else { exp_cor$data <- res$data exp_cor$dims <- dim(exp_cor$data) - exp_cor$attrs[[lon_dim]] <- res$lon - exp_cor$attrs[[lat_dim]] <- res$lat + exp_cor$coords[[lon_dim]] <- res$lon + exp_cor$coords[[lat_dim]] <- res$lat res_s2dv <- list(exp = exp_cor, obs = obs) } diff --git a/R/Interpolation.R b/R/Interpolation.R index 906ba26..5d5f70b 100644 --- a/R/Interpolation.R +++ b/R/Interpolation.R @@ -68,16 +68,16 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr stop("The name of the latitude/longitude dimensions in 'exp$data' must match the parametres 'lat_dim' and 'lon_dim'") } - res <- Interpolation(exp = exp$data, lats = exp$attrs[[lat_dim]], lons = exp$attrs[[lon_dim]], - source_file = exp$attrs$source_files[1], points = points, + res <- Interpolation(exp = exp$data, lats = exp$coords[[lat_dim]], lons = exp$coords[[lon_dim]], + source_file = exp$coords$source_files[1], points = points, method_remap = method_remap, target_grid = target_grid, lat_dim = lat_dim, lon_dim = lon_dim, region = region, method_point_interp = method_point_interp, ncores = ncores) # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data exp$data <- res$data exp$dims <- dim(exp$data) - exp$attrs[[lon_dim]] <- res$lon - exp$attrs[[lat_dim]] <- res$lat + exp$coords[[lon_dim]] <- res$lon + exp$coords[[lat_dim]] <- res$lat res_s2dv <- list(exp = exp, obs = NULL) return(res_s2dv) diff --git a/R/Intlr.R b/R/Intlr.R index 9ce7bea..0a05060 100644 --- a/R/Intlr.R +++ b/R/Intlr.R @@ -114,9 +114,9 @@ CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, p exp_cor_aux <- exp_cor } - res <- Intlr(exp = exp$data, obs = obs$data, exp_cor = exp_cor_aux, exp_lats = exp$attrs[[lat_dim]], - exp_lons = exp$attrs[[lon_dim]], obs_lats = obs$attrs[[lat_dim]], obs_lons = obs$attrs[[lon_dim]], - points = points, source_file_exp = exp$attrs$source_files[1], source_file_obs = obs$attrs$source_files[1], + res <- Intlr(exp = exp$data, obs = obs$data, exp_cor = exp_cor_aux, exp_lats = exp$coords[[lat_dim]], + exp_lons = exp$coords[[lon_dim]], obs_lats = obs$coords[[lat_dim]], obs_lons = obs$coords[[lon_dim]], + points = points, source_file_exp = exp$coords$source_files[1], source_file_obs = obs$coords$source_files[1], target_grid = target_grid, lr_method = lr_method, int_method = int_method, method_point_interp = method_point_interp, predictors = predictors, lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, time_dim = time_dim, @@ -126,22 +126,22 @@ CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, p # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data obs$data <- res$obs obs$dims <- dim(obs$data) - obs$attrs[[lon_dim]] <- res$lon - obs$attrs[[lat_dim]] <- res$lat + obs$coords[[lon_dim]] <- res$lon + obs$coords[[lat_dim]] <- res$lat if (is.null(exp_cor)) { exp$data <- res$data exp$dims <- dim(exp$data) - exp$attrs[[lon_dim]] <- res$lon - exp$attrs[[lat_dim]] <- res$lat + exp$coords[[lon_dim]] <- res$lon + exp$coords[[lat_dim]] <- res$lat res_s2dv <- list(exp = exp, obs = obs) } else { if (identical(lr_method, 'basic') | identical(lr_method, '4nn')) { exp_cor$data <- res$data exp_cor$dims <- dim(exp_cor$data) - exp_cor$attrs[[lon_dim]] <- res$lon - exp_cor$attrs[[lat_dim]] <- res$lat + exp_cor$coords[[lon_dim]] <- res$lon + exp_cor$coords[[lat_dim]] <- res$lat # when large-scale is selected, the forecast object does not have to be an s2dv_cube object } else if (identical(lr_method, 'large-scale')) { exp_cor <- suppressWarnings(s2dv_cube(res$data, lat = res$lat, lon = res$lon)) diff --git a/R/LogisticReg.R b/R/LogisticReg.R index 60aa573..875fcbe 100644 --- a/R/LogisticReg.R +++ b/R/LogisticReg.R @@ -21,6 +21,10 @@ #'@param obs an 's2dv object' with named dimensions containing the observational field. #'The object must have, at least, the dimensions latitude, longitude and start date. The #'object is expected to be already subset for the desired region. +#'@param exp_cor an optional array with named dimensions containing the seasonal forecast +#'experiment data. If the forecast is provided, it will be downscaled using the hindcast and +#'observations; if not provided, the hindcast will be downscaled instead. The default value +#'is NULL. #'@param target_grid a character vector indicating the target grid to be passed to CDO. #'It must be a grid recognised by CDO or a NetCDF file. #'@param int_method a character vector indicating the regridding method to be passed to CDORemap. @@ -86,7 +90,7 @@ #'res <- CST_LogisticReg(exp = exp, obs = obs, int_method = 'bil', target_grid = 'r1280x640', #'probs_cat = c(1/3, 2/3)) #'@export -CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_method = "ens_mean", +CST_LogisticReg <- function(exp, obs, exp_cor = NULL, target_grid, int_method = NULL, log_reg_method = "ens_mean", probs_cat = c(1/3,2/3), return_most_likely_cat = FALSE, points = NULL, method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", member_dim = "member", region = NULL, loocv = TRUE, ncores = NULL) { @@ -99,28 +103,38 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me stop("Parameter 'obs' must be of the class 's2dv_cube'") } - res <- LogisticReg(exp = exp$data, obs = obs$data, exp_lats = exp$attrs[[lat_dim]], - exp_lons = exp$attrs[[lon_dim]], obs_lats = obs$attrs[[lat_dim]], - obs_lons = obs$attrs[[lon_dim]], target_grid = target_grid, + res <- LogisticReg(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, exp_lats = exp$coords[[lat_dim]], + exp_lons = exp$coords[[lon_dim]], obs_lats = obs$coords[[lat_dim]], + obs_lons = obs$coords[[lon_dim]], target_grid = target_grid, probs_cat = probs_cat, return_most_likely_cat = return_most_likely_cat, int_method = int_method, log_reg_method = log_reg_method, points = points, method_point_interp = method_point_interp, lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, - source_file_exp = exp$attrs$source_files[1], source_file_obs = obs$attrs$source_files[1], + source_file_exp = exp$coords$source_files[1], source_file_obs = obs$coords$source_files[1], region = region, loocv = loocv, ncores = ncores) # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data - exp$data <- res$data - exp$dims <- dim(exp$data) - exp$attrs[[lon_dim]] <- res$lon - exp$attrs[[lat_dim]] <- res$lat - obs$data <- res$obs obs$dims <- dim(obs$data) - obs$attrs[[lon_dim]] <- res$lon - obs$attrs[[lat_dim]] <- res$lat + obs$coords[[lon_dim]] <- res$lon + obs$coords[[lat_dim]] <- res$lat + + if (is.null(exp_cor)) { + exp$data <- res$data + exp$dims <- dim(exp$data) + exp$coords[[lon_dim]] <- res$lon + exp$coords[[lat_dim]] <- res$lat + + res_s2dv <- list(exp = exp, obs = obs) + } else { + exp_cor$data <- res$data + exp_cor$dims <- dim(exp_cor$data) + exp_cor$coords[[lon_dim]] <- res$lon + exp_cor$coords[[lat_dim]] <- res$lat + + res_s2dv <- list(exp = exp_cor, obs = obs) + } - res_s2dv <- list(exp = exp, obs = obs) return(res_s2dv) } @@ -147,6 +161,10 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me #'@param obs an array with named dimensions containing the observational field. The object #'must have, at least, the dimensions latitude, longitude and start date. The object is #'expected to be already subset for the desired region. +#'@param exp_cor an optional array with named dimensions containing the seasonal forecast +#'experiment data. If the forecast is provided, it will be downscaled using the hindcast and +#'observations; if not provided, the hindcast will be downscaled instead. The default value +#'is NULL. #'@param exp_lats a numeric vector containing the latitude values in 'exp'. Latitudes must #'range from -90 to 90. #'@param exp_lons a numeric vector containing the longitude values in 'exp'. Longitudes @@ -223,7 +241,7 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me #'obs_lats = obs_lats, obs_lons = obs_lons, int_method = 'bil', target_grid = 'r1280x640', #'probs_cat = c(1/3, 2/3)) #'@export -LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, +LogisticReg <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, int_method = NULL, log_reg_method = "ens_mean", probs_cat = c(1/3,2/3), return_most_likely_cat = FALSE, points = NULL, method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", member_dim = "member", @@ -292,6 +310,32 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target "'member_dim'") } + if (!is.null(exp_cor)) { + + if (is.na(match(sdate_dim, names(dim(exp_cor))))) { + stop("Missing start date dimension in 'exp_cor', or does not match the parameter ", + "'sdate_dim'") + } + + if (is.na(match(member_dim, names(dim(exp_cor))))) { + stop("Missing member dimension in 'exp_cor', or does not match the parameter 'member_dim'") + } + + if (is.na(match(lon_dim, names(dim(exp_cor))))) { + stop("Missing longitude dimension in 'exp_cor', or does not match the parameter ", + "'lon_dim'") + } + + if (is.na(match(lat_dim, names(dim(exp_cor))))) { + stop("Missing latitude dimension in 'exp_cor', or does not match the parameter ", + "'lat_dim'") + } + + if (loocv) { # loocv equal to false to train with the whole hindcast and predict with the forecast + loocv <- FALSE + } + } + # When observations are pointwise if (!is.null(points) & !is.na(match("location", names(dim(obs))))) { point_obs <- T @@ -343,12 +387,24 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target lat_dim = lat_dim, lon_dim = lon_dim, method_point_interp = method_point_interp, region = region, ncores = ncores) + if (!is.null(exp_cor)) { + exp_cor_interpolated <- Interpolation(exp = exp_cor, lats = exp_lats, lons = exp_lons, target_grid = target_grid, + points = points, method_point_interp = method_point_interp, + source_file = source_file_exp, lat_dim = lat_dim, lon_dim = lon_dim, + method_remap = int_method, region = region, ncores = ncores) + } + # compute ensemble mean anomalies if (log_reg_method == "ens_mean") { predictor <- .get_ens_mean_anom(obj_ens = exp_interpolated$data, member_dim = member_dim, sdate_dim = sdate_dim, ncores = ncores) target_dims_predictor <- sdate_dim + if (!is.null(exp_cor)) { + forecast <- .get_ens_mean_anom(obj_ens = exp_cor_interpolated$data, member_dim = member_dim, sdate_dim = sdate_dim, + ncores = ncores) + target_dims_forecast <- sdate_dim + } } else if (log_reg_method == "ens_mean_sd") { @@ -363,6 +419,17 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target names(dim(predictor)) <- c("pred", names(dim(ens_mean_anom))) target_dims_predictor <- c(sdate_dim, "pred") + + if (!is.null(exp_cor)) { + forecast_mean_anom <- .get_ens_mean_anom(obj_ens = exp_cor_interpolated$data, member_dim = member_dim, + sdate_dim = sdate_dim, ncores = ncores) + forecast_sd <- .get_ens_sd(obj_ens = exp_cor_interpolated$data, member_dim = member_dim, ncores = ncores) + forecast <- abind(forecast_mean_anom, forecast_sd, along = 1/2) + names(dim(forecast)) <- c("pred", names(dim(forecast_mean_anom))) + + target_dims_forecast <- c(sdate_dim, "pred") + } + } else if (log_reg_method == "sorted_members") { predictor <- .sort_members(obj_ens = exp_interpolated$data, member_dim = member_dim, ncores = ncores) @@ -385,22 +452,28 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target } # 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)) - } else { - terc <- convert2prob(as.vector(x), prob = probs_cat) - apply(terc, 1, function(r) which (r == 1))}}, - output_dims = sdate_dim, ncores = ncores)$output1 - - - res <- Apply(list(predictor, obs_cat), target_dims = list(target_dims_predictor, sdate_dim), fun = function(x, y) - .log_reg(x = x, y = y, loocv = loocv,probs_cat=probs_cat), output_dims = c(sdate_dim, "category"), ncores = ncores)$output1 - - if (return_most_likely_cat) { - res <- Apply(res, target_dims = c(sdate_dim, "category"), .most_likely_category, - output_dims = sdate_dim, ncores = ncores)$output1 + obs_cat <- Apply(obs_ref, target_dims = sdate_dim, function(x) { + if (!any(!is.na(x))) { + rep(NA,length(x)) + } else { + terc <- convert2prob(as.vector(x), prob = probs_cat) + apply(terc, 1, function(r) which (r == 1)) + } + }, output_dims = sdate_dim, ncores = ncores)$output1 + + if (is.null(exp_cor)) { + res <- Apply(list(predictor, obs_cat), target_dims = list(target_dims_predictor, sdate_dim), fun = function(x, y) + .log_reg(x = x, y = y, loocv = loocv,probs_cat=probs_cat), + output_dims = c(sdate_dim, "category"), ncores = ncores)$output1 + + if (return_most_likely_cat) { + res <- Apply(res, target_dims = c(sdate_dim, "category"), .most_likely_category, + output_dims = sdate_dim, ncores = ncores)$output1 + } + } else { + res <- Apply(list(predictor, obs_cat, forecast), target_dims = list(target_dims_predictor, sdate_dim, target_dims_forecast), + fun = function(x, y, f) .log_reg(x = x, y = y, f = f, loocv = loocv, probs_cat = probs_cat), + output_dims = c(sdate_dim, "category"), ncores = ncores)$output1 } # restore ensemble dimension in observations if it existed originally @@ -456,15 +529,23 @@ obs_cat <- Apply(obs_ref, target_dims = sdate_dim, function(x) { } # atomic functions for logistic regressions -.log_reg <- function(x, y, loocv,probs_cat) { +.log_reg <- function(x, y, f = NULL, loocv, probs_cat) { tmp_df <- data.frame(x = x, y = y) # if the data is all NA, force return return NA if (all(is.na(tmp_df)) | (sum(apply(tmp_df, 2, function(x) !all(is.na(x)))) == 1) | all(is.na(tmp_df$y))) { + if (is.null(f)) { + n1 <- nrow(tmp_df) + } else { + if (is.null(dim(f))) { + n1 <- length(f) + } else { + n1 <- dim(f)[1] + } + } - n1 <- nrow(tmp_df) - n2<- length(probs_cat)+1 + n2 <- length(probs_cat) + 1 res <- matrix(NA, nrow = n1, ncol = n2) } else { @@ -472,7 +553,7 @@ obs_cat <- Apply(obs_ref, target_dims = sdate_dim, function(x) { lm1 <- .train_lr(df = tmp_df, loocv = loocv) # prediction - res <- pred_lr(lm1 = lm1, df = tmp_df, loocv = loocv,probs_cat=probs_cat) + res <- .pred_lr(lm1 = lm1, df = tmp_df, f = f, loocv = loocv, probs_cat = probs_cat) } return(res) } @@ -503,7 +584,7 @@ obs_cat <- Apply(obs_ref, target_dims = sdate_dim, function(x) { #----------------------------------- # Function to apply the logistic regressions. #----------------------------------- -pred_lr <- function(df, lm1, loocv,probs_cat) { +.pred_lr <- function(df, lm1, f, loocv, probs_cat) { require(plyr) @@ -528,22 +609,39 @@ pred_lr <- function(df, lm1, loocv,probs_cat) { colnames(pred_vals)<-c(1,2) } - } else { + } else { # type = class, probs #pred_vals_ls <- lapply(lm1, predict, data = df, type = "probs") #pred_vals <- unlist(pred_vals_ls) - pred_vals <- predict(lm1[[1]], df, type = "probs") + if (is.null(f)) { + pred_vals <- predict(lm1[[1]], df, type = "probs") - if( length(probs_cat)+1==2) - { - pred_vals_dum<-array(NA,dim=c(nrow(df),2)) - pred_vals_dum[,2]<-pred_vals - pred_vals_dum[,1]<-1-pred_vals - pred_vals<-pred_vals_dum - colnames(pred_vals)<-c(1,2) - } - + if( length(probs_cat) + 1 == 2) { + pred_vals_dum<-array(NA,dim=c(nrow(df),2)) + pred_vals_dum[,2]<-pred_vals + pred_vals_dum[,1]<-1-pred_vals + pred_vals<-pred_vals_dum + colnames(pred_vals)<-c(1,2) + } + } else { + if (is.null(dim(f))) { + pred_vals <- predict(lm1[[1]], newdata = data.frame(x = as.vector(f)), type = "probs") + } else { + pred_vals <- predict(lm1[[1]], newdata = data.frame(x = f), type = "probs") + } + if (length(probs_cat) + 1 == 2) { + if (is.null(dim(f))) { + pred_vals_dum <- matrix(NA, nrow = length(f), ncol = 2) + } else { + pred_vals_dum <- matrix(NA, nrow = dim(f)[1], ncol = 2) + } + pred_vals_dum[,2] <- pred_vals + pred_vals_dum[,1] <- 1 - pred_vals + pred_vals <- pred_vals_dum + colnames(pred_vals) <- c(1,2) + } + } } return(pred_vals) diff --git a/inst/examples/analogs.R b/inst/examples/analogs.R index 3446647..a4dcbc2 100644 --- a/inst/examples/analogs.R +++ b/inst/examples/analogs.R @@ -44,8 +44,7 @@ obs_era5 <- startR::Start(dat = '/esarchive/recon/ecmwf/era5/daily_mean/$var$_f1 num_procs = 1, retrieve = TRUE) # Create 's2dv_cube' object -obs_s2dv <- s2dv_cube(obs_era5, lat = attr(obs_era5, "Variables")$dat1$lat, - lon = attr(obs_era5, "Variables")$dat1$lon, source_files = attr(obs_era5, "Files")[1,1,]) +obs_s2dv <- as.s2dv_cube(obs_era5) # ERA5-Land observations: high-resolution observations # change to time = indices(1:28) @@ -90,8 +89,7 @@ exp <- startR::Start(dat = '/esarchive/exp/ecmwf/system5c3s/daily_mean/$var$_f6h num_procs = 1, retrieve = TRUE) # Create 's2dv_cube' object -exp_s2dv <- s2dv_cube(exp, lat = attr(exp, "Variables")$dat1$lat, lon = attr(exp, "Variables")$dat1$lon, - source_files = attr(exp, "Files")[1,1,]) +exp_s2dv <- as.s2dv_cube(exp) #--------------------------- # Analogs with no 'window' dimension in the initial object and with cv diff --git a/inst/examples/interpolation-bc.R b/inst/examples/interpolation-bc.R index 2d8ef6b..0a68e8d 100644 --- a/inst/examples/interpolation-bc.R +++ b/inst/examples/interpolation-bc.R @@ -22,7 +22,8 @@ latmax <- 44.1 #sdates <- c('20000201','20010201','20020201','20030201','20040201','20050201','20060201','20070201') sdates <- format(ymd("20000201") + rep(years(0:20), each=1),"%Y%m%d") -sdates <- format(seq(ymd("20000201"), ymd("20200201"), "1 month"), "%Y%m%d") +#sdates <- format(seq(ymd("20000201"), ymd("20200201"), "1 month"), "%Y%m%d") +sdates_fcst <- format(seq(ymd("20230201"), ymd("20230401"), "1 month"), "%Y%m%d") obs <- startR::Start(dat = '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$var$_$sdate$.nc', var = 'tas', time = indices(1), lat = values(list(latmin, latmax)), @@ -33,7 +34,7 @@ obs <- startR::Start(dat = '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r num_procs = 1, retrieve = TRUE) exp <- startR::Start(dat = '/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc', - var = 'tas', time = indices(1), member = indices(1:3), sdate = sdates, + var = 'tas', time = indices(1), member = 'all', sdate = sdates, lat = values(list(latmin, latmax)), lat_reorder = Sort(decreasing = FALSE), lon = values(list(lonmin, lonmax)), lon_reorder = CircularSort(-180, 180), synonims = list(var = c('var','variable'), lon = c('lon', 'longitude'), @@ -41,13 +42,22 @@ exp <- startR::Start(dat = '/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_f return_vars = list(lat = 'dat', lon = 'dat'), num_procs = 1, retrieve = TRUE) +fcst <- startR::Start(dat = '/esarchive/exp/ecmwf/system51c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', time = indices(1), member = 'all', sdate = sdates_fcst, + lat = values(list(latmin, latmax)), lat_reorder = Sort(decreasing = FALSE), + lon = values(list(lonmin, lonmax)), lon_reorder = CircularSort(-180, 180), + synonims = list(var = c('var','variable'), lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), member = c('member','ensemble')), + return_vars = list(lat = 'dat', lon = 'dat'), + num_procs = 1, retrieve = TRUE) + #------------------------------------ # Downscaling with Intbc #------------------------------------ down_1_sbc <- Intbc(exp = exp, obs = obs, exp_lats = attr(exp, "Variables")$dat1$lat, exp_lons = attr(exp, "Variables")$dat1$lon, obs_lats = attr(obs, "Variables")$dat1$lat, obs_lons = attr(obs, "Variables")$dat1$lon, target_grid = target_grid, int_method = 'con', - bc_method = 'simple_bias', ncores = 4) + bc_method = 'evmos', ncores = 4) down_1_cal <- Intbc(exp = exp, obs = obs, exp_lats = attr(exp, "Variables")$dat1$lat, exp_lons = attr(exp, "Variables")$dat1$lon, obs_lats = attr(obs, "Variables")$dat1$lat, @@ -69,19 +79,9 @@ down_1_dbc <- Intbc(exp = exp, obs = obs, exp_lats = attr(exp, "Variables")$dat1 #------------------------------------ # Transform exp and obs into s2dv_objects #------------------------------------ -obs <- s2dv_cube(obs, lat = attr(obs, "Variables")$dat1$lat, lon = attr(obs, "Variables")$dat1$lon, - source_files = attr(obs, "Files")[1,1,]) -attr(obs$lon,"first_lon") <- obs$lon[1] -attr(obs$lon,"last_lon") <- obs$lon[length(obs$lon)] -attr(obs$lat,"first_lat") <- obs$lat[1] -attr(obs$lat,"last_lat") <- obs$lat[length(obs$lat)] - -exp <- s2dv_cube(exp, lat = attr(exp, "Variables")$dat1$lat, lon = attr(exp, "Variables")$dat1$lon, - source_files = attr(exp, "Files")[1,1,]) -attr(exp$lon,"first_lon") <- exp$lon[1] -attr(exp$lon,"last_lon") <- exp$lon[length(exp$lon)] -attr(exp$lat,"first_lat") <- exp$lat[1] -attr(exp$lat,"last_lat") <- exp$lat[length(exp$lat)] +obs <- as.s2dv_cube(obs) +exp <- as.s2dv_cube(exp) +fcst <- as.s2dv_cube(fcst) #------------------------------------ # Downscaling with CST_Intbc @@ -128,6 +128,19 @@ down_points2_qm <- Intbc(exp = exp, obs = obs_point, exp_lats = attr(exp, "Vari source_file_exp = attr(exp, "Files")[1,1,1], source_file_obs = NULL, target_grid = target_grid, bc_method = 'quantile_mapping', ncores = 4) +#------------------------------------ +# Downscaling forecasts using hindcasts to train +#------------------------------------ +down_1_sbc <- Intbc(exp = exp, obs = obs, exp_cor = fcst, exp_lats = attr(exp, "Variables")$dat1$lat, + exp_lons = attr(exp, "Variables")$dat1$lon, obs_lats = attr(obs, "Variables")$dat1$lat, + obs_lons = attr(obs, "Variables")$dat1$lon, target_grid = target_grid, int_method = 'con', + bc_method = 'evmos', ncores = 4) + +#------------------------------------ +# Downscaling forecasts using hindcasts to train with CST_Intbc +#------------------------------------ +down_2_sbc <- CST_Intbc(exp = exp, obs = obs, exp_cor = fcst, target_grid = target_grid, int_method = 'con', + bc_method = 'evmos', ncores = 4) diff --git a/inst/examples/interpolation-lr.R b/inst/examples/interpolation-lr.R index 012dc3a..25446c8 100644 --- a/inst/examples/interpolation-lr.R +++ b/inst/examples/interpolation-lr.R @@ -4,6 +4,7 @@ library(startR) library(s2dv) library(lubridate) library(multiApply) +library(plyr) source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Interpolation.R') source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Intlr.R') source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Utils.R') @@ -23,6 +24,7 @@ latmax <- 44.1 #sdates <- format(ymd("19811201") + rep(years(0:36),each=1),"%Y%m%d") sdates <- c('20000201','20010201','20020201','20030201','20040201','20050201','20060201','20070201') +sdates_fcst <- c('20230201', '20230401') obs1 <- startR::Start(dat = '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$var$_$sdate$.nc', var = 'tas', time = indices(1), lat = values(list(latmin, latmax)), @@ -33,7 +35,16 @@ obs1 <- startR::Start(dat = '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h- num_procs = 1, retrieve = TRUE) exp1 <- startR::Start(dat = '/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc', - var = 'tas', time = indices(1), member = indices(1:3), sdate = sdates, + var = 'tas', time = indices(1), member = 'all', sdate = sdates, + lat = values(list(latmin, latmax)), lat_reorder = Sort(decreasing = FALSE), + lon = values(list(lonmin, lonmax)), lon_reorder = CircularSort(-180, 180), + synonims = list(var = c('var','variable'), lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), member = c('member','ensemble')), + return_vars = list(lat = 'dat', lon = 'dat'), + num_procs = 1, retrieve = TRUE) + +fcst <- startR::Start(dat = '/esarchive/exp/ecmwf/system51c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', time = indices(1), member = 'all', sdate = sdates_fcst, lat = values(list(latmin, latmax)), lat_reorder = Sort(decreasing = FALSE), lon = values(list(lonmin, lonmax)), lon_reorder = CircularSort(-180, 180), synonims = list(var = c('var','variable'), lon = c('lon', 'longitude'), @@ -44,9 +55,9 @@ exp1 <- startR::Start(dat = '/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_ #---------------------------------- # Downscaling with Intlr #---------------------------------- -ind_rdm <- array(NA, dim = c('sdate' = length(sdates), member = 3, 'vars' = 2)) -ind_rdm[,,1] <- rnorm(n = prod(length(sdates), 3), mean = 0, sd = 1) -ind_rdm[,,2] <- rnorm(n = prod(length(sdates), 3), mean = 0, sd = 1) +ind_rdm <- array(NA, dim = c('sdate' = length(sdates), member = 25, 'vars' = 2)) +ind_rdm[,,1] <- rnorm(n = prod(length(sdates),25), mean = 0, sd = 1) +ind_rdm[,,2] <- rnorm(n = prod(length(sdates),25), mean = 0, sd = 1) down_1_bas <- Intlr(exp = exp1, obs = obs1, exp_lats = attr(exp1, "Variables")$dat1$lat, exp_lons = attr(exp1, "Variables")$dat1$lon, obs_lats = attr(obs1, "Variables")$dat1$lat, @@ -92,19 +103,9 @@ down_2_nn <- Intlr(exp = exp1, obs = obs1, exp_lats = attr(exp1, "Variables")$da #---------------------------------- # Create s2dv objects #---------------------------------- -obs1 <- s2dv_cube(obs1, lat = attr(obs1, "Variables")$dat1$lat, lon = attr(obs1, "Variables")$dat1$lon, - source_files = attr(obs, "Files")[1,1,]) -attr(obs1$lon,"first_lon") <- obs1$lon[1] -attr(obs1$lon,"last_lon") <- obs1$lon[length(obs1$lon)] -attr(obs1$lat,"first_lat") <- obs1$lat[1] -attr(obs1$lat,"last_lat") <- obs1$lat[length(obs1$lat)] - -exp1 <- s2dv_cube(exp1, lat = attr(exp1, "Variables")$dat1$lat, lon = attr(exp1, "Variables")$dat1$lon, - source_files = attr(exp1, "Files")[1,1,]) -attr(exp1$lon,"first_lon") <- exp1$lon[1] -attr(exp1$lon,"last_lon") <- exp1$lon[length(exp1$lon)] -attr(exp1$lat,"first_lat") <- exp1$lat[1] -attr(exp1$lat,"last_lat") <- exp1$lat[length(exp1$lat)] +obs1 <- as.s2dv_cube(obs1) +exp1 <- as.s2dv_cube(exp1) +fcst <- as.s2dv_cube(fcst) #---------------------------------- # Downscaling with CST_Intlr @@ -130,5 +131,40 @@ down_4_bas <- Intlr(exp = exp1, obs = obs_point, exp_lats = attr(exp1, "Variable int_method = 'bilinear', predictors = NULL, source_file_exp = attr(exp1, "Files")[1,1,1], source_file_obs = attr(obs1, "Files")[1,1,1], loocv = TRUE, ncores = 1) +#------------------------------------ +# Downscaling forecasts using hindcasts to train +#------------------------------------ +ind_rdm <- array(NA, dim = c('sdate' = length(sdates), member = 25, 'vars' = 2)) +ind_rdm[,,1] <- rnorm(n = prod(length(sdates),25), mean = 0, sd = 1) +ind_rdm[,,2] <- rnorm(n = prod(length(sdates),25), mean = 0, sd = 1) +ind_rdm_fcst <- array(NA, dim = c('sdate' = length(sdates_fcst), member = 51, 'vars' = 2)) +ind_rdm_fcst[,,1] <- rnorm(n = prod(length(sdates_fcst),51), mean = 0, sd = 1) +ind_rdm_fcst[,,2] <- rnorm(n = prod(length(sdates_fcst),51), mean = 0, sd = 1) +down_5_bas <- Intlr(exp = exp1, obs = obs1, exp_cor = fcst, exp_lats = attr(exp1, "Variables")$dat1$lat, + exp_lons = attr(exp1, "Variables")$dat1$lon, obs_lats = attr(obs1, "Variables")$dat1$lat, + obs_lons = attr(obs1, "Variables")$dat1$lon, target_grid = target_grid, lr_method = "basic", + int_method = "bilinear", predictors = NULL, loocv = TRUE, ncores = 4) +down_5_lsc <- Intlr(exp = exp1, obs = obs1, exp_cor = ind_rdm_fcst, exp_lats = attr(exp1, "Variables")$dat1$lat, + exp_lons = attr(exp1, "Variables")$dat1$lon, obs_lats = attr(obs1, "Variables")$dat1$lat, + obs_lons = attr(obs1, "Variables")$dat1$lon, target_grid = target_grid, lr_method = "large-scale", + int_method = "bilinear", predictors = ind_rdm, loocv = TRUE, ncores = 4) +down_5_nn <- Intlr(exp = exp1, obs = obs1, exp_cor = fcst, exp_lats = attr(exp1, "Variables")$dat1$lat, + exp_lons = attr(exp1, "Variables")$dat1$lon, obs_lats = attr(obs1, "Variables")$dat1$lat, + obs_lons = attr(obs1, "Variables")$dat1$lon, target_grid = target_grid, lr_method = "4nn", + int_method = "bilinear", predictors = NULL, loocv = TRUE, ncores = 4) + +#------------------------------------ +# Downscaling forecasts using hindcasts to train with CST_Intlr +#------------------------------------ +down_6_bas <- CST_Intlr(exp = exp1, obs = obs1, exp_cor = fcst, target_grid = target_grid, lr_method = "basic", + int_method = "bilinear", predictors = NULL, loocv = TRUE, ncores = 4) + +down_6_lsc <- CST_Intlr(exp = exp1, obs = obs1, exp_cor = ind_rdm_fcst, target_grid = target_grid, lr_method = "large-scale", + int_method = "bilinear", predictors = ind_rdm, loocv = TRUE, ncores = 4) + +down_6_nn <- CST_Intlr(exp = exp1, obs = obs1, exp_cor = fcst, target_grid = target_grid, lr_method = "4nn", + int_method = "bilinear", predictors = NULL, loocv = TRUE, ncores = 4) + + diff --git a/inst/examples/interpolation.R b/inst/examples/interpolation.R index 0784398..1761929 100644 --- a/inst/examples/interpolation.R +++ b/inst/examples/interpolation.R @@ -34,13 +34,8 @@ exp1 <- startR::Start(dat = '/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_ member = c('member','ensemble')), return_vars = list(lat = 'dat', lon = 'dat'), num_procs = 1, retrieve = TRUE) -exp2 <- s2dv_cube(exp1, lat = attr(exp1, "Variables")$dat1$lat, lon = attr(exp1, "Variables")$dat1$lon, - source_files = attr(exp1, "Files")[1,1,]) -#attr(exp2$lon,"first_lon") <- exp2$attrs$lon[1] -#attr(exp2$lon,"last_lon") <- exp2$attrs$lon[length(exp2$attrs$lon)] -#attr(exp2$lat,"first_lat") <- exp2$attrs$lat[1] -#attr(exp2$lat,"last_lat") <- exp2$attrs$lat[length(exp2$attrs$lat)] +exp2 <- as.s2dv_cube(exp1) #-------------------------------- # Downscaling with Interpolation diff --git a/inst/examples/logistic-reg.R b/inst/examples/logistic-reg.R index f2e4ba0..6af0041 100644 --- a/inst/examples/logistic-reg.R +++ b/inst/examples/logistic-reg.R @@ -23,6 +23,7 @@ latmax <- 44.1 #sdates <- c('20000201','20010201','20020201','20030201','20040201','20050201','20060201','20070201') sdates <- format(ymd("20000701") + rep(years(0:20), each=1),"%Y%m%d") +sdates_fcst <- c('20230201', '20230401') obs <- startR::Start(dat = '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$var$_$sdate$.nc', var = 'tas', time = indices(1), lat = values(list(latmin, latmax)), @@ -33,7 +34,7 @@ obs <- startR::Start(dat = '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r num_procs = 1, retrieve = TRUE) exp <- startR::Start(dat = '/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc', - var = 'tas', time = indices(1), member = indices(1:3), sdate = sdates, + var = 'tas', time = indices(1), member = 'all', sdate = sdates, lat = values(list(latmin, latmax)), lat_reorder = Sort(decreasing = FALSE), lon = values(list(lonmin, lonmax)), lon_reorder = CircularSort(-180, 180), synonims = list(var = c('var','variable'), lon = c('lon', 'longitude'), @@ -41,6 +42,15 @@ exp <- startR::Start(dat = '/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_f return_vars = list(lat= 'dat', lon = 'dat'), num_procs = 1, retrieve = TRUE) +fcst <- startR::Start(dat = '/esarchive/exp/ecmwf/system51c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', time = indices(1), member = 'all', sdate = sdates_fcst, + lat = values(list(latmin, latmax)), lat_reorder = Sort(decreasing = FALSE), + lon = values(list(lonmin, lonmax)), lon_reorder = CircularSort(-180, 180), + synonims = list(var = c('var','variable'), lon = c('lon', 'longitude'), + lat = c('lat', 'latitude'), member = c('member','ensemble')), + return_vars = list(lat = 'dat', lon = 'dat'), + num_procs = 1, retrieve = TRUE) + # downscaling to fields w LogisticReg() res <- LogisticReg(exp = exp, obs = obs, @@ -55,10 +65,8 @@ res <- LogisticReg(exp = exp, return_most_likely_cat = TRUE) # downscaling to fields w CST_LogisticReg() -obs <- s2dv_cube(obs, lat = attr(obs, "Variables")$dat1$lat, lon = attr(obs, "Variables")$dat1$lon, - source_files = attr(obs, "Files")[1,1,]) -exp <- s2dv_cube(exp, lat = attr(exp, "Variables")$dat1$lat, lon = attr(exp, "Variables")$dat1$lon, - source_files = attr(exp, "Files")[1,1,]) +obs <- as.s2dv_cube(obs) +exp <- as.s2dv_cube(exp) res_cst1 <- CST_LogisticReg(exp = exp, obs = obs, @@ -87,8 +95,9 @@ res_cst3 <- CST_LogisticReg(exp = exp, return_most_likely_cat = TRUE, loocv = TRUE) - +#---------------------------------- # Plot Most Likely Tercile +#---------------------------------- colors <- c("#33BFD1", "#b5b5b5", "#FF764D") PlotEquiMap(res_cst1$data[21,1,1,1,1:29,2:58], lat = res_cst1$lat[1:29], lon = res_cst1$lon[2:58], filled.continents = F, brks = c(0,1,2,3), cols = colors, @@ -108,8 +117,9 @@ PlotEquiMap(exp_anom[1,1,1,,,21], lat = attr(exp, "Variables")$dat1$lat, lon = a filled.continents = F, fileout = "/esarchive/scratch/jramon/downscaling/plots/exp_tas_202007.png", width = 7, height = 5) - +#---------------------------------- # downscaling to point locations +#---------------------------------- points <- list(lat = c(36, 36.1), lon = c(0.6, 1)) source_file <- '/esarchive/exp/ecmwf/system5c3s/monthly_mean/tas_f6h/tas_20000201.nc' method_point_interp <- 'bilinear' @@ -118,3 +128,21 @@ res <- LogisticReg(exp = exp, obs = obs, exp_lats = attr(exp, "Variables")$dat1$ obs_lons = attr(obs, "Variables")$dat1$lon, target_grid = target_grid, int_method = "bil", points = points, source_file = source_file, method_point_interp = method_point_interp) + +#---------------------------------- +# downscaling training with hindcasts and providing forecasts +#---------------------------------- +res_2_em <- LogisticReg(exp = exp, obs = obs, exp_cor = fcst, exp_lats = attr(exp, "Variables")$dat1$lat, + exp_lons = attr(exp, "Variables")$dat1$lon, obs_lats = attr(obs, "Variables")$dat1$lat, + obs_lons = attr(obs, "Variables")$dat1$lon, target_grid = target_grid, int_method = "bilinear", + log_reg_method = "ens_mean", return_most_likely_cat = TRUE, loocv = TRUE, probs_cat = 1:2/3, + ncores = 1) + +res_2_esd <- LogisticReg(exp = exp, obs = obs, exp_cor = fcst, exp_lats = attr(exp, "Variables")$dat1$lat, + exp_lons = attr(exp, "Variables")$dat1$lon, obs_lats = attr(obs, "Variables")$dat1$lat, + obs_lons = attr(obs, "Variables")$dat1$lon, target_grid = target_grid, int_method = "bilinear", + log_reg_method = "ens_mean_sd", return_most_likely_cat = TRUE, loocv = TRUE, probs_cat = 1:2/3, + ncores = 1) + + + -- GitLab From 51ce706e23b9fe0d11cfd6b8d8342dd6b00d3ba7 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Thu, 14 Dec 2023 19:21:18 +0100 Subject: [PATCH 05/17] added error for sorted members, which cannot be used because hindcasts and forecasts have generally different ensemble members --- R/LogisticReg.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/R/LogisticReg.R b/R/LogisticReg.R index 875fcbe..deb587c 100644 --- a/R/LogisticReg.R +++ b/R/LogisticReg.R @@ -430,7 +430,11 @@ LogisticReg <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, target_dims_forecast <- c(sdate_dim, "pred") } - } else if (log_reg_method == "sorted_members") { + } else if (log_reg_method == "sorted_members") { + + if (!is.null(exp_cor)) { + stop('sorted_members method cannot be used to downscale forecasts since the ensemble members are generally exchangeable') + } predictor <- .sort_members(obj_ens = exp_interpolated$data, member_dim = member_dim, ncores = ncores) target_dims_predictor <- c(sdate_dim, member_dim) -- GitLab From 6612627ab72dd47c2f7fb8f9b20a033374ebedf3 Mon Sep 17 00:00:00 2001 From: eduzenli Date: Thu, 4 Jan 2024 16:03:38 +0100 Subject: [PATCH 06/17] small mistakes are corrected --- R/Intbc.R | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/R/Intbc.R b/R/Intbc.R index e78308b..9ba9840 100644 --- a/R/Intbc.R +++ b/R/Intbc.R @@ -66,9 +66,9 @@ #'dim(obs) <- c(lat = 12, lon = 15, sdate = 5) #'obs_lons <- seq(1,5, 4/14) #'obs_lats <- seq(1,4, 3/11) -#'exp <- s2dv_cube(data = exp, lat = exp_lats, lon = exp_lons) -#'obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons) -#'res <- CST_Intbc(exp = exp, obs = obs, target_grid = 'r1280x640', bc_method = 'simple_bias', int_method = 'conservative') +#'exp <- s2dv_cube(data = exp, coords = list(lat = exp_lats, lon = exp_lons)) +#'obs <- s2dv_cube(data = obs, coords = list(lat = obs_lats, lon = obs_lons)) +#'res <- CST_Intbc(exp = exp, obs = obs, target_grid = 'r1280x640', bc_method = 'bias', int_method = 'conservative') #'@export CST_Intbc <- function(exp, obs, exp_cor = NULL, target_grid, bc_method, int_method = NULL, points = NULL, @@ -317,13 +317,16 @@ Intbc <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo lat_dim = lat_dim, lon_dim = lon_dim, method_point_interp = method_point_interp, region = region, ncores = ncores) + + exp_cor_interpolated <- NULL + if (!is.null(exp_cor)) { exp_cor_interpolated <- Interpolation(exp = exp_cor, lats = exp_lats, lons = exp_lons, target_grid = target_grid, method_remap = int_method, points = points, source_file = source_file_exp, lat_dim = lat_dim, lon_dim = lon_dim, method_point_interp = method_point_interp, region = region, ncores = ncores) - } + } # If after interpolating 'exp' data the coordinates do not match, the obs data is interpolated to # the same grid to force the matching -- GitLab From 93ecfc0d25a5035e411192495749c3f3cac0f439 Mon Sep 17 00:00:00 2001 From: eduzenli Date: Thu, 4 Jan 2024 16:09:04 +0100 Subject: [PATCH 07/17] example with s2dv_cube object is updated to prevent coords & attrs confusion --- R/Interpolation.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/Interpolation.R b/R/Interpolation.R index 5d5f70b..4ff9ad7 100644 --- a/R/Interpolation.R +++ b/R/Interpolation.R @@ -48,7 +48,7 @@ #'dim(exp) <- c(member = 5, lat = 4, lon = 5, sdate = 5, time = 1) #'lons <- 1:5 #'lats <- 1:4 -#'exp <- s2dv_cube(data = exp, lat = lats, lon = lons) +#'exp <- s2dv_cube(data = exp, coords = list(lat = lats, lon = lons)) #'res <- CST_Interpolation(exp = exp, method_remap = 'conservative', target_grid = 'r1280x640') #'@export CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_grid = NULL, -- GitLab From 5c82579a9c3f6b2038bea22e28316ac8e42d91ad Mon Sep 17 00:00:00 2001 From: eduzenli Date: Fri, 5 Jan 2024 15:17:49 +0100 Subject: [PATCH 08/17] example with s2dv_cube object is updated to prevent coords & attrs confusion --- R/Analogs.R | 4 ++-- R/Intlr.R | 4 ++-- R/LogisticReg.R | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/R/Analogs.R b/R/Analogs.R index 1cbea15..9bde35b 100644 --- a/R/Analogs.R +++ b/R/Analogs.R @@ -81,8 +81,8 @@ #'dim(obs) <- c(lat = 12, lon = 15, sdate = 5, time = 30) #'obs_lons <- seq(0,6, 6/14) #'obs_lats <- seq(0,6, 6/11) -#'exp <- s2dv_cube(data = exp, lat = exp_lats, lon = exp_lons) -#'obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons) +#'exp <- s2dv_cube(data = exp, coords = list(lat = exp_lats, lon = exp_lons)) +#'obs <- s2dv_cube(data = obs, coords = list(lat = obs_lats, lon = obs_lons)) #'downscaled_field <- CST_Analogs(exp = exp, obs = obs, grid_exp = 'r360x180') #'@export CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analog = NULL, lat_dim = "lat", diff --git a/R/Intlr.R b/R/Intlr.R index 0a05060..b2e6b09 100644 --- a/R/Intlr.R +++ b/R/Intlr.R @@ -86,8 +86,8 @@ #'dim(obs) <- c(lat = 12, lon = 15, sdate = 5) #'obs_lons <- seq(1,5, 4/14) #'obs_lats <- seq(1,4, 3/11) -#'exp <- s2dv_cube(data = exp, lat = exp_lats, lon = exp_lons) -#'obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons) +#'exp <- s2dv_cube(data = exp, coords = list(lat = exp_lats, lon = exp_lons)) +#'obs <- s2dv_cube(data = obs, coords = list(lat = obs_lats, lon = obs_lons)) #'res <- CST_Intlr(exp = exp, obs = obs, target_grid = 'r1280x640', lr_method = 'basic', int_method = 'conservative') #'@export CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, points = NULL, int_method = NULL, diff --git a/R/LogisticReg.R b/R/LogisticReg.R index deb587c..6071190 100644 --- a/R/LogisticReg.R +++ b/R/LogisticReg.R @@ -85,8 +85,8 @@ #'dim(obs) <- c(lat = 12, lon = 15, sdate = 15) #'obs_lons <- seq(1,5, 4/14) #'obs_lats <- seq(1,4, 3/11) -#'exp <- s2dv_cube(data = exp, lat = exp_lats, lon = exp_lons) -#'obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons) +#'exp <- s2dv_cube(data = exp, coords = list(lat = exp_lats, lon = exp_lons)) +#'obs <- s2dv_cube(data = obs, coords = list(lat = obs_lats, lon = obs_lons)) #'res <- CST_LogisticReg(exp = exp, obs = obs, int_method = 'bil', target_grid = 'r1280x640', #'probs_cat = c(1/3, 2/3)) #'@export -- GitLab From b146058467f9e3a828bafa0ec862b6437ecf5a80 Mon Sep 17 00:00:00 2001 From: eduzenli Date: Wed, 17 Jan 2024 10:51:13 +0100 Subject: [PATCH 09/17] Analogs in the branch updated with the Analogs in the master --- R/Analogs.R | 440 +++++++++++++++++++++++++++++++++------------------- 1 file changed, 284 insertions(+), 156 deletions(-) diff --git a/R/Analogs.R b/R/Analogs.R index 9bde35b..062425f 100644 --- a/R/Analogs.R +++ b/R/Analogs.R @@ -1,16 +1,24 @@ #'@rdname CST_Analogs -#'@title Downscaling using Analogs based on large scale fields. +#'@title Downscaling using Analogs based on coarse scale fields. #' #'@author J. Ramon, \email{jaume.ramon@bsc.es} #' #'@description This function performs a downscaling using Analogs. To compute -#'the analogs given a coarse-scale field, the function looks for days with -#'similar conditions in the historical observations. The coarse scale and -#'observation data can be either global or regional. In the latter case, the -#'region is defined by the user. In principle, the coarse and observation data -#'should be of the same variable, although different variables can also be admitted. -#'The analogs function will find the N best analogs based in Minimum Euclidean -#'distance. +#'the analogs given a coarse-scale field, the function looks for days with similar conditions +#'in the historical observations. The analogs function determines the N best analogs based +#'on Euclidian distance, distance correlation, or Spearman's correlation metrics. To downscale +#'a local-scale variable, either the variable itself or another large-scale variable +#'can be utilized as the predictor. In the first scenario, analogs are examined between +#'the observation and model data of the same local-scale variable. In the latter scenario, +#'the function identifies the day in the observation data that closely resembles +#'the large-scale pattern of interest in the model. When it identifies the date of +#'the best analog, the function extracts the corresponding local-scale variable for that day +#'from the observation of the local scale variable. The used local-scale and large-scale +#'variables can be retrieved from independent regions. The input data for the first case must +#'include 'exp' and 'obs,' while in the second case, 'obs,' 'obsL,' and 'exp' are the +#'required input fields. Users can perform the downscaling process over the subregions +#'that can be identified through the 'region' argument, instead of focusing +#'on the entire area of the loaded data. #' #'The search of analogs must be done in the longest dataset posible, but might #'require high-memory computational resources. This is important since it is @@ -21,25 +29,24 @@ #'is intended to downscale climate prediction data (i.e., sub-seasonal, seasonal #'and decadal predictions) but can admit climate projections or reanalyses. It does #'not have constrains of specific region or variables to downscale. -#'@param exp an 's2dv' object with named dimensions containing the experimental field on -#'the coarse scale for which the downscaling is aimed. The object must have, at least, -#'the dimensions latitude, longitude, start date and time. The object is expected to be -#'already subset for the desired region. Data can be in one or two integrated regions, e.g., -#'crossing the Greenwich meridian. To get the correct results in the latter case, -#'the borders of the region should be specified in the parameter 'region'. See parameter -#''region'. -#'@param obs an 's2dv' object with named dimensions containing the observational field. -#'The object must have, at least, the dimensions latitude, longitude and start date. -#'The object is expected to be already subset for the desired region. -#'@param obs2 an 's2dv' object with named dimensions containing a different observational -#'field to that in 'obs'. If provided, these observations will be used in the training, -#'i.e. the searching for analogs, so that they should be in a coarser grid to those in -#''obs'. Training with observations on a grid with a spatial resolution closer to that -#'in 'exp', will in principle ensure better results. The object must have, at least, the -#'dimensions latitude, longitude and start date. The object is expected to be already -#'subset for the desired region. -#'@param grid_exp a character vector with a path to an example file of the exp data. -#'It can be either a path to another NetCDF file which to read the target grid from +#'@param exp an 's2dv_cube' object with named dimensions containing the experimental field +#'on the coarse scale for the variable targeted for downscaling (in case obsL is not provided) +#'or for the large-scale variable used as the predictor (if obsL is provided). +#'The object must have, at least, the dimensions latitude, longitude, start date and time. +#'The object is expected to be already subset for the desired region. Data can be in one +#'or two integrated regions, e.g., crossing the Greenwich meridian. To get the correct +#'results in the latter case, the borders of the region should be specified in the parameter +#''region'. See parameter 'region'. +#'@param obs an 's2dv_cube' object with named dimensions containing the observational field +#'for the variable targeted for downscaling. The object must have, at least, the dimensions +#'latitude, longitude and start date. The object is expected to be already subset for the +#'desired region. +#'@param obsL an 's2dv_cube' object with named dimensions containing the observational +#'field of the large-scale variable. The object must have, at least, the dimensions latitude, +#'longitude and start date. The object is expected to be already subset for the desired region. +#'@param grid_exp a character vector with a path to an example file of the exp (if the +#'predictor is the local scale variable) or expL (if the predictor is a large scale variable) +#'data. It can be either a path to another NetCDF file which to read the target grid from #'(a single grid must be defined in such file) or a character vector indicating the #'coarse grid to be passed to CDO, and it must be a grid recognised by CDO. #'@param nanalogs an integer indicating the number of analogs to be searched @@ -56,6 +63,9 @@ #''data' in exp and obs. Default set to "time". #'@param member_dim a character vector indicating the member dimension name in the element #''data' in exp and obs. Default set to "member". +#'@param metric a character vector to select the analog specification method. Only these +#'options are valid: "dist" (i.e., Euclidian distance), "dcor" (i.e., distance correlation) +#'or "cor" (i.e., Spearman's .correlation). The default metric is "dist". #'@param region a numeric vector indicating the borders of the region defined in exp. #'It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers #'to the left border, while lonmax refers to the right border. latmin indicates the lower @@ -68,7 +78,7 @@ #'of the window. It is recommended to be set to TRUE. Default to TRUE. #'@param ncores an integer indicating the number of cores to use in parallel computation. #'The default value is NULL. -#'@return An 's2dv' object. The element 'data' contains the dowscaled field, 'lat' the +#'@return An 's2dv_cube' object. The element 'data' contains the dowscaled field, 'lat' the #'downscaled latitudes, and 'lon' the downscaled longitudes. If fun_analog is set to NULL #'(default), the output array in 'data' also contains the dimension 'analog' with the best #'analog days. @@ -85,27 +95,36 @@ #'obs <- s2dv_cube(data = obs, coords = list(lat = obs_lats, lon = obs_lons)) #'downscaled_field <- CST_Analogs(exp = exp, obs = obs, grid_exp = 'r360x180') #'@export -CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analog = NULL, lat_dim = "lat", - lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", member_dim = "member", - region = NULL, return_indices = FALSE, loocv_window = TRUE, ncores = NULL) { - - # input exp and obs must be s2dv_cube objects - if (!inherits(exp,'s2dv_cube')) { - stop("Parameter 'exp' must be of the class 's2dv_cube'") - } +CST_Analogs <- function(exp, obs, obsL = NULL, grid_exp, nanalogs = 3, + fun_analog = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", + time_dim = "time", member_dim = "member", metric = "dist", region = NULL, + return_indices = FALSE, loocv_window = TRUE, ncores = NULL) { # input exp and obs must be s2dv_cube objects if (!inherits(obs,'s2dv_cube')) { stop("Parameter 'obs' must be of the class 's2dv_cube'") } - res <- Analogs(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], - exp_lons = exp$coords[[lon_dim]], obs_lats = obs$coords[[lat_dim]], - obs_lons = obs$coords[[lon_dim]], grid_exp = grid_exp, nanalogs = nanalogs, - fun_analog = fun_analog, lat_dim = lat_dim, lon_dim = lon_dim, - sdate_dim = sdate_dim, time_dim = time_dim, member_dim = member_dim, - region = region, return_indices = return_indices, loocv_window = loocv_window, - ncores = ncores) + if (!is.null(obsL)) { + # input obs must be s2dv_cube objects + if (!inherits(obsL,'s2dv_cube')) { + stop("Parameter 'obsL' must be of the class 's2dv_cube'") + } + } + # input exp and obs must be s2dv_cube objects + if (!inherits(exp,'s2dv_cube')) { + stop("Parameter 'exp' must be of the class 's2dv_cube'") + } + + res <- Analogs(exp = exp$data, obs = obs$data, obsL = obsL$data, + exp_lats = exp$coords[[lat_dim]], exp_lons = exp$coords[[lon_dim]], + obs_lats = obs$coords[[lat_dim]], obs_lons = obs$coords[[lon_dim]], + obsL_lats = obsL$coords[[lat_dim]], obsL_lons = obsL$coords[[lon_dim]], + grid_exp = grid_exp, nanalogs = nanalogs, fun_analog = fun_analog, + lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, + time_dim = time_dim, member_dim = member_dim, metric = metric, + region = region, return_indices = return_indices, + loocv_window = loocv_window, ncores = ncores) # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data exp$data <- res$data @@ -129,13 +148,21 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo #'@author Ll. Lledó, \email{llorenc.lledo@ecmwf.int} #' #'@description This function performs a downscaling using Analogs. To compute -#'the analogs given a coarse-scale field, the function looks for days with -#'similar conditions in the historical observations. The coarse scale and -#'observation data can be either global or regional. In the latter case, the -#'region is defined by the user. In principle, the coarse and observation data -#'should be of the same variable, although different variables can also be admitted. -#'The analogs function will find the N best analogs based in Minimum Euclidean -#'distance. +#'the analogs given a coarse-scale field, the function looks for days with similar conditions +#'in the historical observations. The analogs function determines the N best analogs based +#'on RMSE, distance correlation, or Spearman's correlation metrics. To downscale +#'a local-scale variable, either the variable itself or another large-scale variable +#'can be utilized as the predictor. In the first scenario, analogs are examined between +#'the observation and model data of the same local-scale variable. In the latter scenario, +#'the function identifies the day in the observation data that closely resembles +#'the large-scale pattern of interest in the model. When it identifies the date of +#'the best analog, the function extracts the corresponding local-scale variable for that day +#'from the observation of the local scale variable. The used local-scale and large-scale +#'variables can be retrieved from independent regions. The input data for the first case must +#'include 'exp' and 'obs,' while in the second case, 'obs,' 'obsL,' and 'expL' are the +#'required input fields. Users can perform the downscaling process over the subregions +#'that can be identified through the 'region' and 'regionL' arguments, instead of focusing +#'on the entire area of the loaded data. #' #'The search of analogs must be done in the longest dataset posible, but might #'require high-memory computational resources. This is important since it is @@ -146,19 +173,20 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo #'is intended to downscale climate prediction data (i.e., sub-seasonal, seasonal #'and decadal predictions) but can admit climate projections or reanalyses. It does #'not have constrains of specific region or variables to downscale. -#'@param exp an array with named dimensions containing the experimental field on the -#'coarse scale for which the downscaling is aimed. The object must have, at least, -#'the dimensions latitude, longitude, start date and time. The object is expected to be -#'already subset for the desired region. Data can be in one or two integrated regions, e.g., -#'crossing the Greenwich meridian. To get the correct results in the latter case, -#'the borders of the region should be specified in the parameter 'region'. See parameter -#''region'. -#'@param obs an array with named dimensions containing the observational field. The object -#'must have, at least, the dimensions latitude, longitude, start date and time. The object -#'is expected to be already subset for the desired region. Optionally, 'obs' can have the -#'dimension 'window', containing the sampled fields into which the function will look for -#'the analogs. See function 'generate_window()'. Otherwise, the function will look for -#'analogs using all the possible fields contained in obs. +#'@param exp an array with named dimensions containing the experimental field +#'on the coarse scale for the variable targeted for downscaling (in case obsL is not provided) +#'or for the large-scale variable used as the predictor (if obsL is provided). +#'The object must have, at least, the dimensions latitude, longitude, start date and time. +#'The object is expected to be already subset for the desired region. Data can be in one +#'or two integrated regions, e.g., crossing the Greenwich meridian. To get the correct +#'results in the latter case, the borders of the region should be specified in the parameter +#''region'. See parameter 'region'. +#'@param obs an array with named dimensions containing the observational field for the variable +#'targeted for downscaling. The object must have, at least, the dimensions latitude, longitude, +#'start date and time. The object is expected to be already subset for the desired region. +#'Optionally, 'obs' can have the dimension 'window', containing the sampled fields into which +#'the function will look for the analogs. See function 'generate_window()'. Otherwise, +#'the function will look for analogs using all the possible fields contained in obs. #'@param exp_lats a numeric vector containing the latitude values in 'exp'. Latitudes must #'range from -90 to 90. #'@param exp_lons a numeric vector containing the longitude values in 'exp'. Longitudes @@ -167,17 +195,21 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo #'range from -90 to 90. #'@param obs_lons a numeric vector containing the longitude values in 'obs'. Longitudes #'can range from -180 to 180 or from 0 to 360. -#'@param grid_exp a character vector with a path to an example file of the exp data. +#'@param obsL an 's2dv_cube' object with named dimensions containing the observational +#'field of the large-scale variable.The object must have, at least, the dimensions latitude, +#'longitude and start date. The object is expected to be already subset for the desired region. +#'Optionally, 'obsL' can have the dimension 'window', containing the sampled fields into which +#'the function will look for the analogs. See function 'generate_window()'. Otherwise, +#'the function will look for analogs using all the possible fields contained in obs. +#'@param obsL_lats a numeric vector containing the latitude values in 'obsL'. Latitudes must +#'range from -90 to 90. +#'@param obsL_lons a numeric vector containing the longitude values in 'obsL'. Longitudes +#'can range from -180 to 180 or from 0 to 360. +#'@param grid_exp a character vector with a path to an example file of the exp (if the +#'predictor is the local scale variable) or expL (if the predictor is a large scale variable) data. #'It can be either a path to another NetCDF file which to read the target grid from #'(a single grid must be defined in such file) or a character vector indicating the #'coarse grid to be passed to CDO, and it must be a grid recognised by CDO. -#'@param obs2 an 's2dv' object with named dimensions containing a different observational -#'field to that in 'obs'. If provided, these observations will be used in the training, -#'i.e. the searching for analogs, so that they should be in a coarser grid to those in -#''obs'. Training with observations on a grid with a spatial resolution closer to that -#'in 'exp', will in principle ensure better results. The object must have, at least, the -#'dimensions latitude, longitude and start date. The object is expected to be already -#'subset for the desired region. #'@param nanalogs an integer indicating the number of analogs to be searched. #'@param fun_analog a function to be applied over the found analogs. Only these options #'are valid: "mean", "wmean", "max", "min", "median" or NULL. If set to NULL (default), @@ -192,6 +224,9 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo #''data' in exp and obs. Default set to "time". #'@param member_dim a character vector indicating the member dimension name in the element #''data' in exp and obs. Default set to "member". +#'@param metric a character vector to select the analog specification method. Only these +#'options are valid: "dist" (i.e., Euclidian distance), "dcor" (i.e., distance correlation) +#'or "cor" (i.e., Spearman's .correlation). The default metric is "dist". #'@param region a numeric vector indicating the borders of the region defined in exp. #'It consists of four elements in this order: lonmin, lonmax, latmin, latmax. lonmin refers #'to the left border, while lonmax refers to the right border. latmin indicates the lower @@ -209,7 +244,7 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo #'@import multiApply #'@import CSTools #'@importFrom s2dv InsertDim CDORemap -#'@importFrom FNN get.knnx +#'@importFrom energy dcor #' #'@seealso \code{\link[s2dverification]{CDORemap}} #' @@ -229,11 +264,11 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo #'downscaled_field <- Analogs(exp = exp, obs = obs, exp_lats = exp_lats, exp_lons = exp_lons, #'obs_lats = obs_lats, obs_lons = obs_lons, grid_exp = 'r360x180') #'@export -Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp, obs2 = NULL, - obs2_lats = NULL, obs2_lons = NULL, nanalogs = 3, fun_analog = NULL, - lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", - member_dim = "member", region = NULL, return_indices = FALSE, - loocv_window = TRUE, ncores = NULL) { +Analogs <- function(exp, obs, exp_lats = NULL, exp_lons = NULL, obs_lats, obs_lons, + grid_exp, obsL = NULL, obsL_lats = NULL, obsL_lons = NULL, nanalogs = 3, + fun_analog = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", + time_dim = "time", member_dim = "member", metric = "dist", region = NULL, + return_indices = FALSE, loocv_window = TRUE, ncores = NULL) { #----------------------------------- # Checkings #----------------------------------- @@ -264,7 +299,7 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp, stop("Parameter 'time_dim' must be of the class 'character'") } - # Do not allow synonims for lat (latitude), lon (longitude) and time (sdate) dimension names + # Do not allow synonims for lat (latitude), lon (longitude) and time (sdate) dimension names if (is.na(match(lon_dim, names(dim(exp)))) | is.na(match(lon_dim, names(dim(obs))))) { stop("Missing longitude dimension in 'exp' and/or 'obs', or does not match the parameter ", "'lon_dim'") @@ -305,38 +340,41 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp, } } - if (!is.null(obs2)) { + if (!is.null(obsL) ) { + # the code is not yet prepared to handle members in the observations - if (member_dim %in% names(dim(obs2))) { - if (identical(as.numeric(dim(obs2)[member_dim]), 1)) { - obs2 <- ClimProjDiags::Subset(x = obs2, along = member_dim, indices = 1, drop = 'selected') + if (member_dim %in% names(dim(obsL))) { + if (identical(as.numeric(dim(obsL)[member_dim]), 1)) { + obsL <- ClimProjDiags::Subset(x = obsL, along = member_dim, indices = 1, drop = 'selected') } else { - stop("Not implemented for observations with members ('obs2' can have 'member_dim', ", + stop("Not implemented for observations with members ('obsL' can have 'member_dim', ", "but it should be of length = 1).") } } - if (is.null(obs2_lats) | is.null(obs2_lons)) { + if (is.null(obsL_lats) | is.null(obsL_lons)) { stop("Missing latitudes and/or longitudes for the provided training observations. Please ", - "provide them with the parametres 'obs2_lats' and 'obs2_lons'") + "provide them with the parametres 'obsL_lats' and 'obsL_lons'") } - if (is.na(match(lon_dim, names(dim(obs2))))) { - stop("Missing longitude dimension in 'obs2', or does not match the parameter 'lon_dim'") + if (is.na(match(lon_dim, names(dim(obsL))))) { + stop("Missing longitude dimension in 'obsL', or does not match the parameter 'lon_dim'") } - if (is.na(match(lat_dim, names(dim(obs2))))) { - stop("Missing latitude dimension in 'obs2', or does not match the parameter 'lat_dim'") + if (is.na(match(lat_dim, names(dim(obsL))))) { + stop("Missing latitude dimension in 'obsL', or does not match the parameter 'lat_dim'") } - if (is.na(match(sdate_dim, names(dim(obs2))))) { - stop("Missing start date dimension in 'obs2', or does not match the parameter 'sdate_dim'") + if (is.na(match(sdate_dim, names(dim(obsL))))) { + stop("Missing start date dimension in 'obsL', or does not match the parameter 'sdate_dim'") } - if (is.na(match(time_dim, names(dim(obs2))))) { - stop("Missing time dimension in 'obs2', or does not match the parameter 'time_dim'") + if (is.na(match(time_dim, names(dim(obsL))))) { + stop("Missing time dimension in 'obsL', or does not match the parameter 'time_dim'") } + } + ## ncores if (!is.null(ncores)) { if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) | @@ -350,15 +388,20 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp, stopifnot(fun_analog %in% c("mean", "wmean", "max", "min", "median")) } - if (!is.null(obs2)) { - obs_train <- obs2 - obs_train_lats <- obs2_lats - obs_train_lons <- obs2_lons + # metric method to be used to specify the analogs + stopifnot(metric %in% c("cor", "dcor", "dist")) + + + + if (!is.null(obsL)) { + obs_train <- obsL + obs_train_lats <- obsL_lats + obs_train_lons <- obsL_lons } else { obs_train <- obs obs_train_lats <- obs_lats obs_train_lons <- obs_lons - } + } # Correct indices later if cross-validation loocv_correction <- FALSE @@ -366,24 +409,40 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp, loocv_correction <- TRUE } - #----------------------------------- - # Interpolate high-res observations to the coarse grid - #----------------------------------- + # crop downscaling region, if the argument region is provided. + if (!is.null(region) & is.null(obsL)) { + # if a border is equally distant from two different grids, the map will be cropped from the grid having smaller coordinate + + a <- which.min(abs((region[1]-obs_lons))) + b <- which.min(abs((region[2]-obs_lons))) + c <- which.min(abs((region[3]-obs_lats))) + d <- which.min(abs((region[4]-obs_lats))) + obs <- ClimProjDiags::Subset(x = obs, along = list(lon_dim,lat_dim), + indices = list(a:b,c:d), drop = 'selected') + } + if (is.null(region)) { - warning("The borders of the downscaling region have not been provided. Assuming the four borders of the ", - "downscaling region are defined by the first and last elements of the parametres 'exp_lats' and ", - "'exp_lons'.") - region <- c(exp_lons[1], exp_lons[length(exp_lons)], exp_lats[1], exp_lats[length(exp_lats)]) + warning("The borders of the downscaling region have not been provided. ", + "Assuming the four borders of the downscaling region are defined by the ", + "first and last elements of the parameters 'exp_lats' and 'exp_lons'.") + region <- c(exp_lons[1], exp_lons[length(exp_lons)], exp_lats[1], + exp_lats[length(exp_lats)]) } - obs_interpolated <- Interpolation(exp = obs_train, lats = obs_train_lats, lons = obs_train_lons, - target_grid = grid_exp, lat_dim = lat_dim, lon_dim = lon_dim, - method_remap = "conservative", region = region, ncores = ncores) + obs_interpolated <- Interpolation(exp = obs_train, lats = obs_train_lats, lons = obs_train_lons, + target_grid = grid_exp, lat_dim = lat_dim, lon_dim = lon_dim, + method_remap = "conservative", region = region, + ncores = ncores) + # If after interpolating 'obs' data the coordinates do not match, the exp data is interpolated to # the same grid to force the matching - if (!.check_coords(lat1 = as.numeric(obs_interpolated$lat), lat2 = exp_lats, lon1 = as.numeric(obs_interpolated$lon), lon2 = exp_lons)) { - exp_interpolated <- Interpolation(exp = exp, lats = exp_lats, lons = exp_lons, target_grid = grid_exp, - lat_dim = lat_dim, lon_dim = lon_dim, method_remap = "conservative", + if (!.check_coords(lat1 = as.numeric(obs_interpolated$lat), + lat2 = exp_lats, + lon1 = as.numeric(obs_interpolated$lon), + lon2 = exp_lons)) { + exp_interpolated <- Interpolation(exp = exp, lats = exp_lats, lons = exp_lons, + target_grid = grid_exp, lat_dim = lat_dim, + lon_dim = lon_dim, method_remap = "conservative", region = region, ncores = ncores)$data } else { exp_interpolated <- exp @@ -391,40 +450,57 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp, # Create window if user does not have it in the training observations if ( !("window" %in% names(dim(obs_interpolated$data))) ) { - obs_train_interpolated <- .generate_window(obj = obs_interpolated$data, sdate_dim = sdate_dim, - time_dim = time_dim, loocv = loocv_window, ncores = ncores) - obs_hres <- .generate_window(obj = obs, sdate_dim = sdate_dim, time_dim = time_dim, loocv = loocv_window, ncores = ncores) + obs_train_interpolated <- .generate_window(obj = obs_interpolated$data, sdate_dim = sdate_dim, + time_dim = time_dim, loocv = loocv_window, + ncores = ncores) + if (!is.null(obsL)) { + if ( ("window" %in% names(dim(obs))) ) { + stop("Either both obs and obsL should include 'window' dimension or none.") + } + } + obs_hres <- .generate_window(obj = obs, sdate_dim = sdate_dim, time_dim = time_dim, + loocv = loocv_window, ncores = ncores) + } else { obs_train_interpolated <- obs_interpolated$data - dim(obs_train_interpolated) <- dim(obs_train_interpolated)[-which(names(dim(obs_train_interpolated))=="time")] + dim(obs_train_interpolated) <- dim(ClimProjDiags::Subset(x = obs_train_interpolated, + along = time_dim, indices = 1, drop = 'selected')) + + if (!is.null(obsL)) { + if ( !("window" %in% names(dim(obs))) ) { + stop("Either both obs and obsL should include 'window' dimension or none.") + } + } obs_hres <- obs - dim(obs_hres) <- dim(obs_hres)[-which(names(dim(obs_hres))=="time")] + dim(obs_hres) <- dim(ClimProjDiags::Subset(x = obs_hres, + along = time_dim, indices = 1, drop = 'selected')) + } #----------------------------------- # Reshape train and test #----------------------------------- - res.data <- Apply(list(obs_train_interpolated, exp_interpolated, obs_hres), - target_dims = list(c("window", lat_dim, lon_dim), c(lat_dim, lon_dim), - c("window", lat_dim, lon_dim)), - fun = function(tr, te, ob) .analogs(train = tr, test = te, obs_hres = ob, k = nanalogs, - fun_analog = fun_analog), ncores = ncores)$output1 - # Return the indices of the best analogs - if (return_indices) { - res.ind <- Apply(list(obs_train_interpolated, exp_interpolated, obs_hres), + RES <- Apply(list(obs_train_interpolated, exp_interpolated, obs_hres), target_dims = list(c("window", lat_dim, lon_dim), c(lat_dim, lon_dim), c("window", lat_dim, lon_dim)), - fun = function(tr, te, ob) .analogs(train = tr, test = te, obs_hres = ob, k = nanalogs, - fun_analog = fun_analog, return_indices = TRUE), ncores = ncores, output_dims = 'ind')$output1 + fun = function(tr, te, ob) .analogs(train = tr, test = te, obs_hres = ob, + k = nanalogs, metric = metric, fun_analog = fun_analog), + ncores = ncores) ## output1 -> data, output2 -> index, output3 -> metric + + res.data <- RES$output1 + + # Return the indices of the best analogs + if (return_indices) { + res.ind <- RES$output2 # If cross-validation has been applied, correct the indices if (loocv_correction) { nsdates <- dim(res.ind)[names(dim(res.ind)) == sdate_dim] ntimes <- dim(res.ind)[names(dim(res.ind)) == time_dim] - res.ind <- Apply(res.ind, target_dims = c("ind", sdate_dim), function(x) + res.ind <- Apply(res.ind, target_dims = c("index", sdate_dim), function(x) sapply(1:nsdates, function(s) seq(ntimes * nsdates)[ - (ntimes * (s - 1) + 1:ntimes)][x[, s]]), - output_dims = c("ind", sdate_dim), ncores = ncores)$output1 + output_dims = c("index", sdate_dim), ncores = ncores)$output1 } # restore ensemble dimension in observations if it existed originally @@ -433,69 +509,86 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp, } res <- list(data = res.data, ind = res.ind, obs = obs, lon = obs_lons, lat = obs_lats) - } - else { + } else { # restore ensemble dimension in observations if it existed originally - if (restore_ens) { - obs <- s2dv::InsertDim(obs, posdim = 1, lendim = 1, name = member_dim) - } - - res <- list(data = res.data, obs = obs, lon = obs_lons, lat = obs_lats) + if (restore_ens) { + obs <- s2dv::InsertDim(obs, posdim = 1, lendim = 1, name = member_dim) + } + res <- list(data = res.data, obs = obs, lon = obs_lons, lat = obs_lats) } return(res) } # For each element in test, find the indices of the k nearest neigbhors in train -.analogs <- function(train, test, obs_hres, k, fun_analog, return_indices = FALSE) { - - require(FNN) +.analogs <- function(train, test, obs_hres, k, fun_analog, metric = NULL, return_indices = FALSE) { - # train and obs_hres dim: 3 dimensions window, lat and lon (in this order) + # train and obs_rs dim: 3 dimensions window, lat and lon (in this order) # test dim: 2 dimensions lat and lon (in this order) # Number of lats/lons of the high-resolution data space_dims_hres <- dim(obs_hres)[c(2,3)] # Reformat train and test as an array with (time, points) - train <- apply(train, 1, as.vector); names(dim(train))[1] <- "space" + train <- apply(train, 1, as.vector); names(dim(train))[1] <- "space" test <- as.vector(test) obs_hres <- apply(obs_hres, 1, as.vector); names(dim(obs_hres))[1] <- "space" - + # Identify and remove NA's - idx_na_tr <- is.na(train[ , 1]) + dum<-which(!apply(train,2,function (x) all(is.na(x))))[1] ## the column in which NA in space will be investigated. it shouldn't be all-NA time-step + idx_na_tr <- is.na(train[ , dum]) # NA in space + idy_na_tr <- is.na(train[1, ]) # NA in time idx_na_te <- is.na(test) idx_na <- idx_na_tr | idx_na_te - tr_wo_na <- t(train[!idx_na , ]) + tr_wo_na <- t(train[!idx_na , !idy_na_tr ]) te_wo_na <- test[!idx_na] te_wo_na <- InsertDim(data = te_wo_na, posdim = 1, lendim = 1, name = "time") - - knn.ind <- get.knnx(tr_wo_na, te_wo_na, k) - dist <- knn.ind$nn.dist - idx <- knn.ind$nn.index - - # Either return only the indices or the analogs - if (return_indices) { - res <- as.numeric(idx) + if (all(is.na(test))) { + res <- array(NA, space_dims_hres) } else { + if (metric == "dist") { + dist_all <- sqrt(rowSums((sweep(tr_wo_na, 2, te_wo_na[1,]))^2)) # euc dist + best_vals <- head(sort(dist_all), k) + idx <- match(best_vals, dist_all) + } else if (metric == "cor") { + cor_all <- apply(tr_wo_na, 1, function (x) cor(x,te_wo_na[1, ], method = "spearman")) # cor + best_vals <- head(sort(cor_all, decreasing = TRUE), k) + idx <- match(best_vals, cor_all) + } else if (metric == "dcor") { +# require(energy) + dcor_all <- apply(tr_wo_na, 1, function (x) .dcor(x,te_wo_na[1, ])) # dcor + best_vals <- head(sort(dcor_all, decreasing = TRUE), k,) + idx <- match(best_vals, dcor_all) + } + if (isTRUE(any(idy_na_tr))) { + dum <-(1:length(idy_na_tr))[!idy_na_tr] + idx <- dum[idx] + } + res_ind <- array(idx, k) + names(dim(res_ind)) <- c("index") + res_metric <- array(best_vals, c(k)) + names(dim(res_metric)) <- c("metric") res <- obs_hres[ , idx] dim(res) <- c(space_dims_hres, analogs = k) if (!is.null(fun_analog)) { if (fun_analog == "wmean") { - weight <- 1 / dist + if (metric == "dist") { + weight <- 1 / best_vals + } else { + weight <- best_vals + } res <- apply(res, c(1,2), function(x) weighted.mean(x, weight)) } else if (fun_analog == "min") { - res<-res[,,which.min(dist)] + res <- res[, , which.min(best_vals)] } else if (fun_analog == "max") { - res<-res[,,which.max(dist)] + res <- res[, , which.max(best_vals)] } else { res <- apply(res, c(1,2), fun_analog) } } + return(list(res, res_ind, res_metric)) } - - return(res) } # Add the dimension window to an array that contains, at least, the start date and time @@ -547,3 +640,38 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp, return(obj_window) } + +# Distance correlation function +.dcor <- function(x, y) { + n <- length(x) + + # Calculate Euclidean distances for x + dist_x <- as.matrix(dist(matrix(x))) + # Calculate Euclidean distances for y + dist_y <- as.matrix(dist(matrix(y))) + + # Centering matrices + H <- diag(n) - 1/n + + # Centered distance matrices + dist_centered_x <- H %*% dist_x %*% H + dist_centered_y <- H %*% dist_y %*% H + + # Calculate the product of mean-centered distance matrices + B <- dist_centered_x %*% t(dist_centered_y) + C <- dist_centered_x %*% t(dist_centered_x) + D <- dist_centered_y %*% t(dist_centered_y) + + # Calculate the distance covariance + dcov_xy <- sqrt(sum(diag(B))) + + # Calculate the variances + cov_xx <- sqrt(sum(diag(C))) + cov_yy <- sqrt(sum(diag(D))) + + # Calculate the distance correlation + dcor_xy <- dcov_xy / sqrt(cov_xx * cov_yy) + + return(dcor_xy) +} + -- GitLab From 2e0317305a20f09184be13dd3067496c0d0622ca Mon Sep 17 00:00:00 2001 From: eduzenli Date: Wed, 17 Jan 2024 12:37:48 +0100 Subject: [PATCH 10/17] forecast data can be downscaled via the analogs method --- R/Analogs.R | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/R/Analogs.R b/R/Analogs.R index 062425f..86f9a3d 100644 --- a/R/Analogs.R +++ b/R/Analogs.R @@ -36,7 +36,8 @@ #'The object is expected to be already subset for the desired region. Data can be in one #'or two integrated regions, e.g., crossing the Greenwich meridian. To get the correct #'results in the latter case, the borders of the region should be specified in the parameter -#''region'. See parameter 'region'. +#''region'. See parameter 'region'. Also, the object can be either hindcast or forecast data. +#'However, if forecast data is provided, the loocv_window parameter should be selected as FALSE. #'@param obs an 's2dv_cube' object with named dimensions containing the observational field #'for the variable targeted for downscaling. The object must have, at least, the dimensions #'latitude, longitude and start date. The object is expected to be already subset for the @@ -180,7 +181,8 @@ CST_Analogs <- function(exp, obs, obsL = NULL, grid_exp, nanalogs = 3, #'The object is expected to be already subset for the desired region. Data can be in one #'or two integrated regions, e.g., crossing the Greenwich meridian. To get the correct #'results in the latter case, the borders of the region should be specified in the parameter -#''region'. See parameter 'region'. +#''region'. See parameter 'region'. Also, the object can be either hindcast or forecast data. +#'However, if forecast data is provided, the loocv_window parameter should be selected as FALSE. #'@param obs an array with named dimensions containing the observational field for the variable #'targeted for downscaling. The object must have, at least, the dimensions latitude, longitude, #'start date and time. The object is expected to be already subset for the desired region. @@ -390,8 +392,13 @@ Analogs <- function(exp, obs, exp_lats = NULL, exp_lons = NULL, obs_lats, obs_lo # metric method to be used to specify the analogs stopifnot(metric %in% c("cor", "dcor", "dist")) - + if ( (dim(exp)[sdate_dim] != dim(obs)[sdate_dim]) & loocv_window ) { + loocv_window <- FALSE + warning("The lengths of sdate_dim in the data provided through exp and obs are ", + "not the same. Thus, exp is considered as forecast data, and loocv_window ", + "is set to FALSE.") + } if (!is.null(obsL)) { obs_train <- obsL -- GitLab From 741fd89952b972a90cdad87c641b1b7fbcbe9c99 Mon Sep 17 00:00:00 2001 From: eduzenli Date: Tue, 30 Jan 2024 16:13:57 +0100 Subject: [PATCH 11/17] return value problem for entire-space-NA was solved --- R/Analogs.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/R/Analogs.R b/R/Analogs.R index 86f9a3d..7b484b0 100644 --- a/R/Analogs.R +++ b/R/Analogs.R @@ -552,6 +552,10 @@ Analogs <- function(exp, obs, exp_lats = NULL, exp_lons = NULL, obs_lats, obs_lo if (all(is.na(test))) { res <- array(NA, space_dims_hres) + res_ind <- array(NA, k) + names(dim(res_ind)) <- c("index") + res_metric <- array(NA, k) + names(dim(res_metric)) <- c("metric") } else { if (metric == "dist") { dist_all <- sqrt(rowSums((sweep(tr_wo_na, 2, te_wo_na[1,]))^2)) # euc dist @@ -594,8 +598,8 @@ Analogs <- function(exp, obs, exp_lats = NULL, exp_lons = NULL, obs_lats, obs_lo res <- apply(res, c(1,2), fun_analog) } } - return(list(res, res_ind, res_metric)) } + return(list(res, res_ind, res_metric)) } # Add the dimension window to an array that contains, at least, the start date and time -- GitLab From da95fd7d0532cf6a534ba1e73dad22c0807d084f Mon Sep 17 00:00:00 2001 From: eduzenli Date: Thu, 29 Feb 2024 13:19:27 +0100 Subject: [PATCH 12/17] arrangements about warnings etc --- R/Intlr.R | 41 ++++++++++++++++++++++++++++++++--------- R/LogisticReg.R | 2 +- 2 files changed, 33 insertions(+), 10 deletions(-) diff --git a/R/Intlr.R b/R/Intlr.R index b2e6b09..b2de7bb 100644 --- a/R/Intlr.R +++ b/R/Intlr.R @@ -90,10 +90,11 @@ #'obs <- s2dv_cube(data = obs, coords = list(lat = obs_lats, lon = obs_lons)) #'res <- CST_Intlr(exp = exp, obs = obs, target_grid = 'r1280x640', lr_method = 'basic', int_method = 'conservative') #'@export -CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, points = NULL, int_method = NULL, - method_point_interp = NULL, predictors = NULL, lat_dim = "lat", lon_dim = "lon", - sdate_dim = "sdate", time_dim = "time", member_dim = "member", - large_scale_predictor_dimname = 'vars', loocv = TRUE, region = NULL, ncores = NULL) { +CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, points = NULL, + int_method = NULL, method_point_interp = NULL, predictors = NULL, + lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", + member_dim = "member", large_scale_predictor_dimname = 'vars', + loocv = TRUE, region = NULL, ncores = NULL) { if (!inherits(exp,'s2dv_cube')) { stop("Parameter 'exp' must be of the class 's2dv_cube'") @@ -107,11 +108,14 @@ CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, p if (!inherits(exp_cor,'s2dv_cube')) { stop("Parameter 'exp_cor' must be of the class 's2dv_cube'") } - - exp_cor_aux <- exp_cor$data + exp_cor_aux <- exp_cor$data # when large-scale is selected, the forecast object does not have to be an s2dv_cube object } else if (identical(lr_method, 'large-scale')) { - exp_cor_aux <- exp_cor + if (!inherits(exp_cor,'s2dv_cube')) { + exp_cor_aux <- exp_cor + } else { + exp_cor_aux <- exp_cor$data + } } res <- Intlr(exp = exp$data, obs = obs$data, exp_cor = exp_cor_aux, exp_lats = exp$coords[[lat_dim]], @@ -144,7 +148,14 @@ CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, p exp_cor$coords[[lat_dim]] <- res$lat # when large-scale is selected, the forecast object does not have to be an s2dv_cube object } else if (identical(lr_method, 'large-scale')) { - exp_cor <- suppressWarnings(s2dv_cube(res$data, lat = res$lat, lon = res$lon)) + if (!inherits(exp_cor,'s2dv_cube')) { + exp_cor <- suppressWarnings(s2dv_cube(res$data, lat = res$lat, lon = res$lon)) + } else { + exp_cor$data <- res$data + exp_cor$dims <- dim(exp_cor$data) + exp_cor$coords[[lon_dim]] <- res$lon + exp_cor$coords[[lat_dim]] <- res$lat + } } res_s2dv <- list(exp = exp_cor, obs = obs) @@ -321,6 +332,8 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo if (loocv) { # loocv equal to false to train with the whole hindcast and predict with the forecast loocv <- FALSE + warning("Forecast data will be downscaled. Thus, loocv_window is set to FALSE ", + "to train with the whole hindcast and predict with the forecast.") } if (is.null(predictors)) { @@ -337,6 +350,12 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo } } + if (lr_method == '4nn') { + warning("4nn method skips the interpolation step since the method itself inherently downscale ", + "the coarse resolution data to the reference data resolution without the need ", + "for interpolation.") + } + # When observations are pointwise if (!is.null(points) & !is.na(match("location", names(dim(obs))))) { point_obs <- T @@ -489,7 +508,11 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo if (!is.null(exp_cor)) { aux_dim <- large_scale_predictor_dimname - forecast <- exp_cor + if (!inherits(exp_cor,'s2dv_cube')) { + forecast <- exp_cor + } else { + forecast <- exp_cor$data + } target_dims_predictor <- c(sdate_dim, large_scale_predictor_dimname, member_dim) target_dims_forecast <- c(sdate_dim, large_scale_predictor_dimname, member_dim) } diff --git a/R/LogisticReg.R b/R/LogisticReg.R index 6071190..59b1bef 100644 --- a/R/LogisticReg.R +++ b/R/LogisticReg.R @@ -110,7 +110,7 @@ CST_LogisticReg <- function(exp, obs, exp_cor = NULL, target_grid, int_method = int_method = int_method, log_reg_method = log_reg_method, points = points, method_point_interp = method_point_interp, lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, - source_file_exp = exp$coords$source_files[1], source_file_obs = obs$coords$source_files[1], + source_file_exp = exp$attrs$source_files[1], source_file_obs = obs$attrs$source_files[1], region = region, loocv = loocv, ncores = ncores) # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data -- GitLab From 1fbce1d95286a7d8567ba6d34401f8ef83e4509d Mon Sep 17 00:00:00 2001 From: eduzenli Date: Wed, 6 Mar 2024 14:15:17 +0100 Subject: [PATCH 13/17] .find_nn function was corrected --- R/Intlr.R | 109 +++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 76 insertions(+), 33 deletions(-) diff --git a/R/Intlr.R b/R/Intlr.R index 1c6ef08..b4fe394 100644 --- a/R/Intlr.R +++ b/R/Intlr.R @@ -538,13 +538,16 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo predictor <- .find_nn(coar = exp, lats_hres = obs_lats, lons_hres = obs_lons, lats_coar = exp_lats, lons_coar = exp_lons, lat_dim = lat_dim, - lon_dim = lon_dim, nn = 4, ncores = ncores) + lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, + nn = 4, ncores = ncores) if (!is.null(exp_cor)) { aux_dim <- 'nn' forecast <- .find_nn(coar = exp_cor, lats_hres = obs_lats, lons_hres = obs_lons, lats_coar = exp_lats, lons_coar = exp_lons, lat_dim = lat_dim, - lon_dim = lon_dim, nn = 4, ncores = ncores) + lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, + nn = 4, ncores = ncores) + } if (is.null(points) | ("location" %in% names(dim(obs)))) { @@ -555,21 +558,31 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo lats <- obs_lats lons <- obs_lons - } + } else { # If the downscaling is to point locations: Once the 4 nearest neighbours have been found, # interpolate to point locations - else { + predictor <- Interpolation(exp = predictor, lats = obs_lats, lons = obs_lons, - target_grid = NULL, points = points, - method_point_interp = method_point_interp, + lon_dim = lon_dim, lat_dim = lat_dim, target_grid = NULL, + points = points, method_point_interp = method_point_interp, source_file = source_file_obs, method_remap = NULL, region = region, ncores = ncores) - predictand <- Interpolation(exp = obs, lats = obs_lats, lons = obs_lons, target_grid = NULL, + predictand <- Interpolation(exp = obs, lats = obs_lats, lons = obs_lons, + lon_dim = lon_dim, lat_dim = lat_dim, target_grid = NULL, points = points, method_point_interp = method_point_interp, source_file = source_file_obs, method_remap = NULL, region = region, ncores = ncores) - + + if (!is.null(exp_cor)) { + forecast <- Interpolation(exp = forecast, lats = obs_lats, lons = obs_lons, + lon_dim = lon_dim, lat_dim = lat_dim, target_grid = NULL, + points = points, method_point_interp = method_point_interp, + source_file = source_file_obs, method_remap = NULL, + region = region, ncores = ncores) + forecast <- forecast$data + } + lats <- predictor$lat lons <- predictor$lon predictor <- predictor$data @@ -584,9 +597,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo } else { target_dims_predictor <- c(sdate_dim,'nn') } - } - - else { + } else { stop(paste0(lr_method, " method is not implemented yet")) } @@ -725,21 +736,22 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo if (!is.null(aux_dim)) { fcst_df <- adply(f,.margins = 3, .id = NULL) # 4nn - if (identical(aux_dim, 'nn')) { - colnames(fcst_df) <- paste0('x.',1:4) +# if (identical(aux_dim, 'nn')) { + colnames(fcst_df) <- paste0('x.',1:ncol(fcst_df)) pred_vals_ls <- lapply(lm1, predict, newdata = fcst_df) pred_vals <- unlist(pred_vals_ls) dim(pred_vals) <- dim(f)[names(dim(f)) != aux_dim] - } else { - pred_vals_ls <- lapply(lm1, predict, newdata = data.frame(x = fcst_df)) - pred_vals <- unlist(pred_vals_ls) - dim(pred_vals) <- dim(f)[names(dim(f)) != aux_dim] - } +# } else { +# colnames(fcst_df) <- paste0('x.',1:ncol(fcst_df)) +# pred_vals_ls <- lapply(lm1, predict, newdata = data.frame(x = fcst_df)) +# pred_vals <- unlist(pred_vals_ls) +# dim(pred_vals) <- dim(f)[names(dim(f)) != aux_dim] +# } # basic } else { pred_vals_ls <- lapply(lm1, predict, newdata = data.frame(x = as.vector(f))) pred_vals <- unlist(pred_vals_ls) - dim(pred_vals) <- dim(f) + dim(pred_vals) <- dim(f ) } } } else { @@ -763,19 +775,50 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo # 'coar' is an array with named dimensions #----------------------------------- .find_nn <- function(coar, lats_hres, lons_hres, lats_coar, lons_coar, - lat_dim, lon_dim, nn = 4, ncores = NULL) { - - # Sort the distances from closest to furthest - idx_lat <- as.array(sapply(lats_hres, function(x) order(abs(lats_coar - x))[1:nn])) - idx_lon <- as.array(sapply(lons_hres, function(x) order(abs(lons_coar - x))[1:nn])) - - names(dim(idx_lat)) <- c('nn', lat_dim) - names(dim(idx_lon)) <- c('nn', lon_dim) - - # obtain the values of the nearest neighbours - nearest <- Apply(list(coar, idx_lat, idx_lon), - target_dims = list(c(lat_dim, lon_dim), lat_dim, lon_dim), - fun = function(x, y, z) x[y, z], ncores = ncores)$output1 - + lat_dim, lon_dim, sdate_dim, member_dim, nn = 4, ncores = NULL) { + + # order lats and lons of the high res and coarse res data in 2-dimension in order to calculate the distances to each other. + hres_mat <- expand.grid(as.numeric(lats_hres), as.numeric(lons_hres)) + coar_mat <- expand.grid(as.numeric(lats_coar), as.numeric(lons_coar)) + + # calculate distances and obtain the closest 4 coarse res grid for each high res grid + dist_id <- apply (proxy::dist(hres_mat, coar_mat), 1, order)[1:nn,] + + names(dim(dist_id)) <- c("nn", lon_dim) + + idh <- SplitDim(dist_id, split_dim = lon_dim, + indices = rep(1:(length(lons_hres)),length(lats_hres)), + new_dim_name = c(lat_dim)) + + idh <- aperm (idh, c(1,3,2)) # dimension: nn, lat_dim, lon_dim + + #coar + idc_lat_grid <- match( coar_mat[,1], as.numeric(lats_coar)) + idc_lon_grid <- match( coar_mat[,2], as.numeric(lons_coar)) + cgrid <- cbind(idc_lat_grid,idc_lon_grid) # For the coarse res data, include the lat and lon order each grids ordered in 2-dimension. + + idh_dum <- array(NA, dim = c(4, 2, length(lats_hres), length(lons_hres))) + + for (i in 1:length(lats_hres)) { + for (j in 1:length(lons_hres)) { + idh_dum[,,i,j] <- cgrid[idh[,i,j],] + } + } + + idh <- idh_dum + rm(idh_dum) + names(dim(idh)) <- c("nn", "grd", lat_dim, lon_dim) + + # idh includes the lat and lon order of the nearest neighbours in the coarse res data specified for each lat lon of the high res data. + + nearest <- Apply(list(coar, idh), + target_dims = list(c(sdate_dim, lat_dim, lon_dim, member_dim), + "grd" ), + fun = function(x, y) { + a <- x[,y[1],y[2],] + return(a) + }, + ncores = ncores)$output1 + return(nearest) } -- GitLab From bc96bc322864a424174b029d7bdef49ec55ce3e6 Mon Sep 17 00:00:00 2001 From: eduzenli Date: Wed, 6 Mar 2024 16:00:45 +0100 Subject: [PATCH 14/17] dwn of forecast feature has been added. --- R/Interpolation.R | 4 +++- R/Intlr.R | 28 ++++++++++++---------------- R/LogisticReg.R | 2 +- 3 files changed, 16 insertions(+), 18 deletions(-) diff --git a/R/Interpolation.R b/R/Interpolation.R index b372c1c..51de9ec 100644 --- a/R/Interpolation.R +++ b/R/Interpolation.R @@ -679,7 +679,9 @@ Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, me # sdate, time, ensemble, num_procs, etc) #====================== .get_model_data <- function(weights.df, mdata, ncores = NULL) { - + + require(plyr) + #----------------- # Get data for all combinations of i and j. # (inefficient, getting many unneded pairs). diff --git a/R/Intlr.R b/R/Intlr.R index b4fe394..554131f 100644 --- a/R/Intlr.R +++ b/R/Intlr.R @@ -547,7 +547,6 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo lats_coar = exp_lats, lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, nn = 4, ncores = ncores) - } if (is.null(points) | ("location" %in% names(dim(obs)))) { @@ -598,6 +597,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo target_dims_predictor <- c(sdate_dim,'nn') } } else { + stop(paste0(lr_method, " method is not implemented yet")) } @@ -734,21 +734,14 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo # case to downscale forecasts } else { if (!is.null(aux_dim)) { + # 4nn & large-scale fcst_df <- adply(f,.margins = 3, .id = NULL) - # 4nn -# if (identical(aux_dim, 'nn')) { - colnames(fcst_df) <- paste0('x.',1:ncol(fcst_df)) - pred_vals_ls <- lapply(lm1, predict, newdata = fcst_df) - pred_vals <- unlist(pred_vals_ls) - dim(pred_vals) <- dim(f)[names(dim(f)) != aux_dim] -# } else { -# colnames(fcst_df) <- paste0('x.',1:ncol(fcst_df)) -# pred_vals_ls <- lapply(lm1, predict, newdata = data.frame(x = fcst_df)) -# pred_vals <- unlist(pred_vals_ls) -# dim(pred_vals) <- dim(f)[names(dim(f)) != aux_dim] -# } + colnames(fcst_df) <- paste0('x.',1:ncol(fcst_df)) + pred_vals_ls <- lapply(lm1, predict, newdata = fcst_df) + pred_vals <- unlist(pred_vals_ls) + dim(pred_vals) <- dim(f)[names(dim(f)) != aux_dim] + } else { # basic - } else { pred_vals_ls <- lapply(lm1, predict, newdata = data.frame(x = as.vector(f))) pred_vals <- unlist(pred_vals_ls) dim(pred_vals) <- dim(f ) @@ -771,7 +764,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo } #----------------------------------- -# Function to find N nearest neighbours. +# Function to ind N nearest neighbours. # 'coar' is an array with named dimensions #----------------------------------- .find_nn <- function(coar, lats_hres, lons_hres, lats_coar, lons_coar, @@ -815,7 +808,10 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo target_dims = list(c(sdate_dim, lat_dim, lon_dim, member_dim), "grd" ), fun = function(x, y) { - a <- x[,y[1],y[2],] + a <- x[, y[1], y[2], , drop = FALSE] + a <- Subset (a, along = c(lat_dim, lon_dim), + indices = list (1,1), + drop = "selected") # drop lat lon dim coming from coar data return(a) }, ncores = ncores)$output1 diff --git a/R/LogisticReg.R b/R/LogisticReg.R index f0477bc..5484faa 100644 --- a/R/LogisticReg.R +++ b/R/LogisticReg.R @@ -586,7 +586,7 @@ LogisticReg <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, # prediction res <- .pred_lr(lm1 = lm1, df = tmp_df, f = f, loocv = loocv, probs_cat = probs_cat) - if ( !(is.null(f)) ) { ## if forecast is provided, but syear of forecast is 1. + if ( !(is.null(f)) ) { ## if forecast is provided, and syear of forecast is 1. if ( nrow(f) == 1 ) { res <- array(res, dim = c(1, length(probs_cat) + 1)) } -- GitLab From 470f195564071d9e4a02cc16bddca80f302c646f Mon Sep 17 00:00:00 2001 From: eduzenli Date: Mon, 11 Mar 2024 11:04:25 +0100 Subject: [PATCH 15/17] documentation revision --- R/Intlr.R | 4 ++-- R/LogisticReg.R | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/R/Intlr.R b/R/Intlr.R index 554131f..0d4f10d 100644 --- a/R/Intlr.R +++ b/R/Intlr.R @@ -339,8 +339,8 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo if (loocv) { # loocv equal to false to train with the whole hindcast and predict with the forecast loocv <- FALSE - warning("Forecast data will be downscaled. Thus, loocv_window is set to FALSE ", - "to train with the whole hindcast and predict with the forecast.") + warning("Forecast data will be downscaled. 'loocv' is set to FALSE ", + "to train the mode with the whole hindcast and predict with the forecast.") } if (is.null(predictors)) { diff --git a/R/LogisticReg.R b/R/LogisticReg.R index 5484faa..bf47a7d 100644 --- a/R/LogisticReg.R +++ b/R/LogisticReg.R @@ -337,8 +337,8 @@ LogisticReg <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, if (loocv) { # loocv equal to false to train with the whole hindcast and predict with the forecast loocv <- FALSE - warning("Forecast data will be downscaled. Thus, loocv_window is set to FALSE ", - "to train with the whole hindcast and predict with the forecast.") + warning("Forecast data will be downscaled. 'loocv' is set to FALSE ", + "to train the model with the whole hindcast and predict with the forecast.") } } -- GitLab From be860cc2660c6e8c0e5392ecb9fe6e8c223dbb03 Mon Sep 17 00:00:00 2001 From: eduzenli Date: Tue, 12 Mar 2024 11:59:22 +0100 Subject: [PATCH 16/17] fine tuning --- R/LogisticReg.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/LogisticReg.R b/R/LogisticReg.R index bf47a7d..f74f2ee 100644 --- a/R/LogisticReg.R +++ b/R/LogisticReg.R @@ -479,14 +479,14 @@ 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)) } else { terc <- convert2prob(as.vector(x), prob = probs_cat) - apply(terc, 1, function(r) which (r == 1)) - } - }, output_dims = sdate_dim, ncores = ncores)$output1 + as.integer(apply(terc, 1, function(r) which (r == 1)))}}, + output_dims = sdate_dim, ncores = ncores)$output1 if (is.null(exp_cor)) { res <- Apply(list(predictor, obs_cat), target_dims = list(target_dims_predictor, sdate_dim), -- GitLab From 8914229ca8e16fa1facb38e4e2dac28afdafcd03 Mon Sep 17 00:00:00 2001 From: eduzenli Date: Tue, 12 Mar 2024 12:00:55 +0100 Subject: [PATCH 17/17] 4nn method is converted into 9nn with principal component analysis --- R/Intlr.R | 101 ++++++++++++++++++++++++++++++++---------------------- 1 file changed, 60 insertions(+), 41 deletions(-) diff --git a/R/Intlr.R b/R/Intlr.R index 0d4f10d..953031f 100644 --- a/R/Intlr.R +++ b/R/Intlr.R @@ -26,16 +26,20 @@ #'forecast experiment data. If the forecast is provided, it will be downscaled using the #'hindcast and observations; if not provided, the hindcast will be downscaled instead. The #'default value is NULL. -#'@param lr_method a character vector indicating the linear regression method to be applied -#'after the interpolation. Accepted methods are 'basic', 'large-scale' and '4nn'. The 'basic' -#'method fits a linear regression using high resolution observations as predictands and the -#'interpolated model data as predictor. Then, the regression equation is to the interpolated -#'model data to correct the interpolated values. The 'large-scale' method fits a linear -#'regression with large-scale predictors from the same model (e.g. teleconnection indices) -#'as predictors and high-resolution observations as predictands. This equation is then -#'applied to the interpolated model values. Finally, the '4nn' method uses a linear -#'regression with the four nearest neighbours as predictors and high-resolution observations -#'as predictands. It is then applied to model data to correct the interpolated values. +#'@param lr_method a character vector indicating the linear regression method to be applied. +#'Accepted methods are 'basic', 'large-scale' and '9nn'. The 'basic' method fits a +#'linear regression using high resolution observations as predictands and the +#'interpolated model data as predictor. Then, the regression equation is applied to the +#'interpolated model data to correct the interpolated values. The 'large-scale' method +#'fits a linear regression with large-scale predictors (e.g. teleconnection indices) as +#'predictors and high-resolution observations as predictands. This equation is then applied +#'to the interpolated model values. Finally, the '9nn' method uses a linear regression +#'with the nine nearest neighbours as predictors and high-resolution observations as +#'predictands. Instead of constructing a regression model using all nine predictors, +#'principal component analysis is applied to the data of neighboring grids to reduce the +#'dimension of the predictors. The linear regression model is then built using the principal +#'components that explain 95% of the variance. The '9nn' method does not require a +#'pre-interpolation process. #'@param target_grid a character vector indicating the target grid to be passed to CDO. #'It must be a grid recognised by CDO or a NetCDF file. #'@param points a list of two elements containing the point latitudes and longitudes @@ -107,7 +111,7 @@ CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, p exp_cor_aux <- NULL if (!is.null(exp_cor)) { - if (identical(lr_method, 'basic') | identical(lr_method, '4nn')) { + if (identical(lr_method, 'basic') | identical(lr_method, '9nn')) { if (!inherits(exp_cor, 's2dv_cube')) { stop("Parameter 'exp_cor' must be of the class 's2dv_cube'") } @@ -147,7 +151,7 @@ CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, p res_s2dv <- list(exp = exp, obs = obs) } else { - if (identical(lr_method, 'basic') | identical(lr_method, '4nn')) { + if (identical(lr_method, 'basic') | identical(lr_method, '9nn')) { exp_cor$data <- res$data exp_cor$dims <- dim(exp_cor$data) exp_cor$coords[[lon_dim]] <- res$lon @@ -206,16 +210,20 @@ CST_Intlr <- function(exp, obs, exp_cor = NULL, lr_method, target_grid = NULL, p #'range from -90 to 90. #'@param obs_lons a numeric vector containing the longitude values in 'obs'. Longitudes #'can range from -180 to 180 or from 0 to 360. -#'@param lr_method a character vector indicating the linear regression method to be applied -#'after the interpolation. Accepted methods are 'basic', 'large-scale' and '4nn'. The 'basic' -#'method fits a linear regression using high resolution observations as predictands and the -#'interpolated model data as predictor. Then, the regression equation is to the interpolated -#'model data to correct the interpolated values. The 'large-scale' method fits a linear -#'regression with large-scale predictors from the same model (e.g. teleconnection indices) -#'as predictors and high-resolution observations as predictands. This equation is then -#'applied to the interpolated model values. Finally, the '4nn' method uses a linear -#'regression with the four nearest neighbours as predictors and high-resolution observations -#'as predictands. It is then applied to model data to correct the interpolated values. +#'@param lr_method a character vector indicating the linear regression method to be applied. +#'Accepted methods are 'basic', 'large-scale' and '9nn'. The 'basic' method fits a +#'linear regression using high resolution observations as predictands and the +#'interpolated model data as predictor. Then, the regression equation is applied to the +#'interpolated model data to correct the interpolated values. The 'large-scale' method +#'fits a linear regression with large-scale predictors (e.g. teleconnection indices) as +#'predictors and high-resolution observations as predictands. This equation is then applied +#'to the interpolated model values. Finally, the '9nn' method uses a linear regression +#'with the nine nearest neighbours as predictors and high-resolution observations as +#'predictands. Instead of constructing a regression model using all nine predictors, +#'principal component analysis is applied to the data of neighboring grids to reduce the +#'dimension of the predictors. The linear regression model is then built using the principal +#'components that explain 95% of the variance. The '9nn' method does not require a +#'pre-interpolation process. #'@param target_grid a character vector indicating the target grid to be passed to CDO. #'It must be a grid recognised by CDO or a NetCDF file. #'@param points a list of two elements containing the point latitudes and longitudes @@ -340,7 +348,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo if (loocv) { # loocv equal to false to train with the whole hindcast and predict with the forecast loocv <- FALSE warning("Forecast data will be downscaled. 'loocv' is set to FALSE ", - "to train the mode with the whole hindcast and predict with the forecast.") + "to train the model with the whole hindcast and predict with the forecast.") } if (is.null(predictors)) { @@ -357,8 +365,8 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo } } - if (lr_method == '4nn') { - warning("4nn method skips the interpolation step since the method itself inherently downscale ", + if (lr_method == '9nn') { + warning("9nn method skips the interpolation step since the method itself inherently downscale ", "the coarse resolution data to the reference data resolution without the need ", "for interpolation.") } @@ -429,7 +437,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo #----------------------------------- # Interpolation #----------------------------------- - if (lr_method != '4nn') { + if (lr_method != '9nn') { if (is.null(int_method)) { stop("Parameter 'int_method' must be a character vector indicating the interpolation method. ", @@ -534,19 +542,19 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo # Multi-linear regression with the four nearest neighbours # Predictors: model data # Predictand: observations - else if (lr_method == '4nn') { + else if (lr_method == '9nn') { predictor <- .find_nn(coar = exp, lats_hres = obs_lats, lons_hres = obs_lons, lats_coar = exp_lats, lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, - nn = 4, ncores = ncores) + nn = 9, ncores = ncores) if (!is.null(exp_cor)) { aux_dim <- 'nn' forecast <- .find_nn(coar = exp_cor, lats_hres = obs_lats, lons_hres = obs_lons, lats_coar = exp_lats, lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, - nn = 4, ncores = ncores) + nn = 9, ncores = ncores) } if (is.null(points) | ("location" %in% names(dim(obs)))) { @@ -558,7 +566,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo lats <- obs_lats lons <- obs_lons } else { - # If the downscaling is to point locations: Once the 4 nearest neighbours have been found, + # If the downscaling is to point locations: Once the 9 nearest neighbours have been found, # interpolate to point locations predictor <- Interpolation(exp = predictor, lats = obs_lats, lons = obs_lons, @@ -625,7 +633,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo # Return a list of three elements res <- list(data = res, obs = predictand, lon = lons, lat = lats) - # for testing 4nn + # for testing 9nn #x <- predictor[10,10,1,1,1,,,] #x <- Reorder(x,c(2,3,1)) #y <- predictand[1,1,1,10,,10] @@ -670,12 +678,23 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo } } else { # training -# if (!is.null(f)) { ## too many number of ensemble members, - ## so Principle component analysis is needed. -# PR <- prcomp(tmp_df[, -ncol(tmp_df)], scale = TRUE) ## prcomp between predictors -# N_predictors <- which (summary(PR)$importance[ "Cumulative Proportion",] > .8)[1] -# tmp_df <- data.frame(x = as.vector(PR$x[, 1:N_predictors]), y = y) -# } + + ## apply principle components if method is 9nn + if (any(names(dim(x)) == "nn")) { + PR <- prcomp(~ x.1 + x.2 + x.3 + x.4 + x.5 + x.6 + x.7 + x.8 + x.9, + data = tmp_df[, -ncol(tmp_df)], + na.action = na.omit, scale = FALSE) ## prcomp between predictors + ## 95% of variablity should be explained by the principle components + N_predictors <- which (summary(PR)$importance[ "Cumulative Proportion",] > .95)[1] + PR_mat <- PR$x[, 1:N_predictors, drop = FALSE] + if (!is.null(PR$na.action)) { ## if there are NAs, relocate them to the correct order + dum <- array(NA, dim = c(nrow(tmp_df), ncol(PR_mat) )) + dum [as.numeric(rownames(PR_mat)), ] <- PR_mat + PR_mat <- dum + } + tmp_df <- data.frame(x = PR_mat, y = y) + colnames(tmp_df) <- c(paste0("x.",1:(ncol(tmp_df)-1)), "y") + } lm1 <- .train_lm(df = tmp_df, loocv = loocv) @@ -734,7 +753,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo # case to downscale forecasts } else { if (!is.null(aux_dim)) { - # 4nn & large-scale + # 9nn & large-scale fcst_df <- adply(f,.margins = 3, .id = NULL) colnames(fcst_df) <- paste0('x.',1:ncol(fcst_df)) pred_vals_ls <- lapply(lm1, predict, newdata = fcst_df) @@ -768,7 +787,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo # 'coar' is an array with named dimensions #----------------------------------- .find_nn <- function(coar, lats_hres, lons_hres, lats_coar, lons_coar, - lat_dim, lon_dim, sdate_dim, member_dim, nn = 4, ncores = NULL) { + lat_dim, lon_dim, sdate_dim, member_dim, nn = 9, ncores = NULL) { # order lats and lons of the high res and coarse res data in 2-dimension in order to calculate the distances to each other. hres_mat <- expand.grid(as.numeric(lats_hres), as.numeric(lons_hres)) @@ -781,7 +800,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo idh <- SplitDim(dist_id, split_dim = lon_dim, indices = rep(1:(length(lons_hres)),length(lats_hres)), - new_dim_name = c(lat_dim)) + new_dim_name = lat_dim) idh <- aperm (idh, c(1,3,2)) # dimension: nn, lat_dim, lon_dim @@ -790,7 +809,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo idc_lon_grid <- match( coar_mat[,2], as.numeric(lons_coar)) cgrid <- cbind(idc_lat_grid,idc_lon_grid) # For the coarse res data, include the lat and lon order each grids ordered in 2-dimension. - idh_dum <- array(NA, dim = c(4, 2, length(lats_hres), length(lons_hres))) + idh_dum <- array(NA, dim = c(nn, 2, length(lats_hres), length(lons_hres))) for (i in 1:length(lats_hres)) { for (j in 1:length(lons_hres)) { -- GitLab