Commits (2)
......@@ -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)
......
......@@ -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
......
......@@ -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)
......
......@@ -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)
}
......
......@@ -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)
......