From e578e2798eeeceba25525cf93b7bda95347ddf9d Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 30 May 2023 14:28:10 +0200 Subject: [PATCH 1/2] Modify obs loading for multivar; add metadata_dims = 'var' --- modules/Loading/Loading_decadal.R | 94 +++++++++++++++++++++---------- 1 file changed, 63 insertions(+), 31 deletions(-) diff --git a/modules/Loading/Loading_decadal.R b/modules/Loading/Loading_decadal.R index 43b3c54f..626afeef 100644 --- a/modules/Loading/Loading_decadal.R +++ b/modules/Loading/Loading_decadal.R @@ -133,7 +133,8 @@ load_datasets <- function(recipe) { if (length(variable) > 1) { Start_default_arg_list <- c(Start_default_arg_list, list(table = table, grid = grid, version = version, - table_depends = 'var', grid_depends = 'var', version_depends = 'var')) + table_depends = 'var', grid_depends = 'var', version_depends = 'var', + metadata_dims = 'var')) } if (!multi_path) { @@ -277,9 +278,15 @@ load_datasets <- function(recipe) { #------------------------------------------- # Step 3. Load the reference #------------------------------------------- - obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src, store.freq, - paste0(variable, archive$Reference[[ref.name]][[store.freq]][[variable]])) - obs.files <- paste0('$var$_$file_date$.nc') +#-------NEW-------- + obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src, + store.freq, "$var$$var_dir$", "$var$_$file_date$.nc") + var_dir_obs <- archive$Reference[[ref.name]][[store.freq]][variable] # list(tas = "_f1h-r1440x721cds", tos = "_f1h-r1440x721cds") + +# obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src, store.freq, +# paste0(variable, archive$Reference[[ref.name]][[store.freq]][[variable]])) +# obs.files <- paste0('$var$_$file_date$.nc') +#--------NEW_END------------- # Get from startR_cube # dates <- attr(hcst, 'Variables')$common$time @@ -299,8 +306,12 @@ load_datasets <- function(recipe) { # Restore correct dimensions dim(dates) <- dim(dates_file) - obs <- Start(dat = file.path(obs.path, obs.files), +#--------NEW---------- + obs <- Start(dat = obs.path, var = variable, + var_dir = var_dir_obs, + var_dir_depends = 'var', +#--------NEW_END------------ file_date = unique(format(dates, '%Y%m')), time = dates, # [sday, sweek, syear, time] time_across = 'file_date', @@ -327,8 +338,12 @@ load_datasets <- function(recipe) { # Method 2: reshape hcst time attr's date into an array with time dim then as obs date selector #//////////////// - obs <- Start(dat = file.path(obs.path, obs.files), +#-----------NEW------------ + obs <- Start(dat = obs.path, var = variable, + var_dir = var_dir_obs, + var_dir_depends = 'var', +#--------NEW_END------------ file_date = dates_file, #dates_arr, # [sday, sweek, syear, time] split_multiselected_dims = TRUE, latitude = values(list(lats.min, lats.max)), @@ -344,6 +359,9 @@ load_datasets <- function(recipe) { longitude = c('lon','longitude')), return_vars = list(latitude = NULL, longitude = NULL, time = 'file_date'), +#-------NEW---------- + metadata_dims = 'var', +#-------NEW_END--------- silent = !DEBUG, retrieve = TRUE) } @@ -353,6 +371,9 @@ load_datasets <- function(recipe) { # sday sweek syear time # 1 1 2 14 + # Remove var_dir dimension + obs <- Subset(obs, along = "var_dir", indices = 1, drop = "selected") + # Only ensemble dim could be different if (!identical(dim(obs), dim(hcst$data)[-9])) { error(recipe$Run$logger, @@ -480,36 +501,47 @@ load_datasets <- function(recipe) { #------------------------------------------- # Remove negative values in accumulative variables dictionary <- read_yaml("conf/variable-dictionary.yml") - if (dictionary$vars[[variable]]$accum) { - info(recipe$Run$logger, - " Accumulated variable: setting negative values to zero.") - obs$data[obs$data < 0] <- 0 - hcst$data[hcst$data < 0] <- 0 - if (!is.null(fcst)) { - fcst$data[fcst$data < 0] <- 0 - } - } - # Convert precipitation to mm/day - ## TODO: Make a function? - if (variable == "prlr") { - # Verify that the units are m/s and the same in obs and hcst - if (((attr(obs$Variable, "variable")$units == "m s-1") || - (attr(obs$Variable, "variable")$units == "m s**-1")) && - ((attr(hcst$Variable, "variable")$units == "m s-1") || - (attr(hcst$Variable, "variable")$units == "m s**-1"))) { - +#----------NEW---------- + for (var_idx in 1:length(variable)) { + var_name <- variable[var_idx] + if (dictionary$vars[[var_name]]$accum) { info(recipe$Run$logger, - "Converting precipitation from m/s to mm/day.") - obs$data <- obs$data*86400*1000 - attr(obs$Variable, "variable")$units <- "mm/day" - hcst$data <- hcst$data*86400*1000 - attr(hcst$Variable, "variable")$units <- "mm/day" + paste0("Accumulated variable ", var_name, + ": setting negative values to zero.")) + # obs$data[, var_idx, , , , , , , ] <- pmax(Subset(obs$data, + # along = "var", + # indices = var_idx, F), 0) + obs$data[, var_idx, , , , , , , ][obs$data[, var_idx, , , , , , , ] < 0] <- 0 + hcst$data[, var_idx, , , , , , , ][hcst$data[, var_idx, , , , , , , ] < 0] <- 0 if (!is.null(fcst)) { - fcst$data <- fcst$data*86400*1000 - attr(fcst$Variable, "variable")$units <- "mm/day" + fcst$data[, var_idx, , , , , , , ][fcst$data[, var_idx, , , , , , , ] < 0] <- 0 + } + } + + # Convert prlr from m/s to mm/day + ## TODO: Make a unit conversion function + if (variable[[var_idx]] == "prlr") { + # Verify that the units are m/s and the same in obs and hcst + if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") || + (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) && + ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") || + (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) { + info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.") + obs$data[, var_idx, , , , , , , ] <- + obs$data[, var_idx, , , , , , , ]*86400*1000 + obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day" + hcst$data[, var_idx, , , , , , , ] <- + hcst$data[, var_idx, , , , , , , ]*86400*1000 + hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" + if (!is.null(fcst)) { + fcst$data[, var_idx, , , , , , , ] <- + fcst$data[, var_idx, , , , , , , ]*86400*1000 + fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" + } } } } +#---------NEW_END--------- #------------------------------------------- # Step 6. Print summary -- GitLab From d437f9d0bfe6cc1f847600f8f26d1d8787648fe4 Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 30 May 2023 14:43:52 +0200 Subject: [PATCH 2/2] Remove comment lines --- modules/Loading/Loading_decadal.R | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/modules/Loading/Loading_decadal.R b/modules/Loading/Loading_decadal.R index 626afeef..74c97e29 100644 --- a/modules/Loading/Loading_decadal.R +++ b/modules/Loading/Loading_decadal.R @@ -278,7 +278,6 @@ load_datasets <- function(recipe) { #------------------------------------------- # Step 3. Load the reference #------------------------------------------- -#-------NEW-------- obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src, store.freq, "$var$$var_dir$", "$var$_$file_date$.nc") var_dir_obs <- archive$Reference[[ref.name]][[store.freq]][variable] # list(tas = "_f1h-r1440x721cds", tos = "_f1h-r1440x721cds") @@ -286,7 +285,6 @@ load_datasets <- function(recipe) { # obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src, store.freq, # paste0(variable, archive$Reference[[ref.name]][[store.freq]][[variable]])) # obs.files <- paste0('$var$_$file_date$.nc') -#--------NEW_END------------- # Get from startR_cube # dates <- attr(hcst, 'Variables')$common$time @@ -306,12 +304,10 @@ load_datasets <- function(recipe) { # Restore correct dimensions dim(dates) <- dim(dates_file) -#--------NEW---------- obs <- Start(dat = obs.path, var = variable, var_dir = var_dir_obs, var_dir_depends = 'var', -#--------NEW_END------------ file_date = unique(format(dates, '%Y%m')), time = dates, # [sday, sweek, syear, time] time_across = 'file_date', @@ -338,12 +334,10 @@ load_datasets <- function(recipe) { # Method 2: reshape hcst time attr's date into an array with time dim then as obs date selector #//////////////// -#-----------NEW------------ obs <- Start(dat = obs.path, var = variable, var_dir = var_dir_obs, var_dir_depends = 'var', -#--------NEW_END------------ file_date = dates_file, #dates_arr, # [sday, sweek, syear, time] split_multiselected_dims = TRUE, latitude = values(list(lats.min, lats.max)), @@ -359,9 +353,7 @@ load_datasets <- function(recipe) { longitude = c('lon','longitude')), return_vars = list(latitude = NULL, longitude = NULL, time = 'file_date'), -#-------NEW---------- metadata_dims = 'var', -#-------NEW_END--------- silent = !DEBUG, retrieve = TRUE) } @@ -501,7 +493,6 @@ load_datasets <- function(recipe) { #------------------------------------------- # Remove negative values in accumulative variables dictionary <- read_yaml("conf/variable-dictionary.yml") -#----------NEW---------- for (var_idx in 1:length(variable)) { var_name <- variable[var_idx] if (dictionary$vars[[var_name]]$accum) { @@ -541,7 +532,6 @@ load_datasets <- function(recipe) { } } } -#---------NEW_END--------- #------------------------------------------- # Step 6. Print summary -- GitLab