diff --git a/example_scripts/example_downscaling.R b/example_scripts/example_downscaling.R index 4e7e6294a5519b263cc3473b9d749435ef867dd3..1e9f1d6ba13ae328eae149624fb70df9c75bbb4e 100644 --- a/example_scripts/example_downscaling.R +++ b/example_scripts/example_downscaling.R @@ -6,6 +6,7 @@ # Load modules source("modules/Loading/Loading.R") source("modules/Units/Units.R") +source("modules/Anomalies/Anomalies.R") source("modules/Downscaling/Downscaling.R") source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") @@ -19,6 +20,8 @@ recipe <- prepare_outputs(recipe_file) data <- Loading(recipe) # Change units data <- Units(recipe, data) +# Compute anomalies +data <- Anomalies(recipe, data) # Downscale datasets data <- Downscaling(recipe, data) # Compute skill metrics diff --git a/modules/Anomalies/Anomalies.R b/modules/Anomalies/Anomalies.R index a2f0c0053791cde810f52a87379508ee227bf0b5..0eb9fbc7f5ece4fb5e99e4a98262dccb2561c469 100644 --- a/modules/Anomalies/Anomalies.R +++ b/modules/Anomalies/Anomalies.R @@ -12,7 +12,6 @@ Anomalies <- function(recipe, data) { "'Workflow:Anomalies:compute' is missing from the recipe.")) stop() } - if (recipe$Analysis$Workflow$Anomalies$compute) { if (recipe$Analysis$Workflow$Anomalies$cross_validation) { cross <- TRUE @@ -23,28 +22,51 @@ Anomalies <- function(recipe, data) { } original_dims <- data$hcst$dim - # Compute anomalies - anom <- CST_Anomaly(data$hcst, data$obs, - cross = cross, - memb = TRUE, - memb_dim = 'ensemble', - dim_anom = 'syear', - dat_dim = c('dat', 'ensemble'), - ftime_dim = 'time', - ncores = recipe$Analysis$ncores) - # Reorder dims - anom$exp$data <- Reorder(anom$exp$data, names(original_dims)) - anom$obs$data <- Reorder(anom$obs$data, names(original_dims)) - # Save full fields hcst_fullvalue <- data$hcst obs_fullvalue <- data$obs - - # Hindcast climatology - - data$hcst <- anom$exp - data$obs <- anom$obs - remove(anom) + # Compute anomalies + if (isTRUE(all.equal(as.vector(data$hcst$coords$latitude), + as.vector(data$obs$coords$latitude))) && + isTRUE(all.equal(as.vector(data$hcst$coords$longitude), + as.vector(data$obs$coord$longitude)))) { + anom <- CST_Anomaly(data$hcst, data$obs, + cross = cross, + memb = TRUE, + memb_dim = 'ensemble', + dim_anom = 'syear', + dat_dim = c('dat', 'ensemble'), + ftime_dim = 'time', + ncores = recipe$Analysis$ncores) + # Reorder dims + anom$exp$data <- Reorder(anom$exp$data, names(original_dims)) + anom$obs$data <- Reorder(anom$obs$data, names(original_dims)) + # Hindcast climatology + data$hcst <- anom$exp + data$obs <- anom$obs + remove(anom) + } else { + ## TODO: Remove when cross-validation is implemented for this use case + if (cross) { + warn(recipe$Run$logger, + paste("Anomaly computation in cross-validation has not been", + "implemented yet for the case where hcst and obs have", + "different grids. The climatology will be computed as", + "a simple mean.")) + } + clim_hcst <- Apply(data$hcst$data, + target_dims = c('syear', 'ensemble'), + mean, + na.rm = recipe$Analysis$remove_NAs, + ncores = recipe$Analysis$ncores)$output1 + clim_obs <- Apply(data$obs$data, + target_dims = c('syear', 'ensemble'), + mean, + na.rm = recipe$Analysis$remove_NAs, + ncores = recipe$Anaysis$ncores)$output1 + data$hcst$data <- Ano(data = data$hcst$data, clim = clim_hcst) + data$obs$data <- Ano(data = data$obs$data, clim = clim_obs) + } # Change variable metadata for (var in data$hcst$attrs$Variable$varName) { # Change hcst longname @@ -57,15 +79,20 @@ Anomalies <- function(recipe, data) { # Compute forecast anomaly field if (!is.null(data$fcst)) { # Compute hindcast climatology ensemble mean - clim <- s2dv::Clim(hcst_fullvalue$data, obs_fullvalue$data, - time_dim = "syear", - dat_dim = c("dat", "ensemble"), - memb = FALSE, - memb_dim = "ensemble", - ftime_dim = "time", - ncores = recipe$Analysis$ncores) - clim_hcst <- InsertDim(clim$clim_exp, posdim = 1, lendim = 1, - name = "syear") + if (isTRUE(all.equal(as.vector(data$hcst$coords$latitude), + as.vector(data$obs$coords$latitude))) && + isTRUE(all.equal(as.vector(data$hcst$coords$longitude), + as.vector(data$obs$coord$longitude)))) { + clim <- s2dv::Clim(hcst_fullvalue$data, obs_fullvalue$data, + time_dim = "syear", + dat_dim = c("dat", "ensemble"), + memb = FALSE, + memb_dim = "ensemble", + ftime_dim = "time", + ncores = recipe$Analysis$ncores) + clim_hcst <- InsertDim(clim$clim_exp, posdim = 1, lendim = 1, + name = "syear") + } # Store original dimensions dims <- dim(clim_hcst) # Repeat the array as many times as ensemble members @@ -113,7 +140,7 @@ Anomalies <- function(recipe, data) { } else { warn(recipe$Run$logger, paste("The Anomalies module has been called, but", - "recipe parameter Analysis:Variables:anomaly is set to FALSE.", + "recipe parameter Workflow:anomalies:compute is set to FALSE.", "The full fields will be returned.")) hcst_fullvalue <- NULL obs_fullvalue <- NULL diff --git a/recipes/atomic_recipes/recipe_seasonal_downscaling.yml b/recipes/atomic_recipes/recipe_seasonal_downscaling.yml index 219c5d6168f299153f83ce5fed2b7b9c957a927d..c526140d4cb91a3fd15a3e52a6d86b5178581fe6 100644 --- a/recipes/atomic_recipes/recipe_seasonal_downscaling.yml +++ b/recipes/atomic_recipes/recipe_seasonal_downscaling.yml @@ -29,7 +29,7 @@ Analysis: type: none Workflow: Anomalies: - compute: no # yes/no, default yes + compute: yes # yes/no, default yes cross_validation: no # yes/no, default yes save: 'none' # 'all'/'none'/'exp_only'/'fcst_only' Calibration: diff --git a/tests/recipes/recipe-seasonal_downscaling.yml b/tests/recipes/recipe-seasonal_downscaling.yml index 41df9e82ce41e02630f4bd13a4a7042f650e7301..994f0bd8181ea47b3659d7fc287b7e6fae3beda5 100644 --- a/tests/recipes/recipe-seasonal_downscaling.yml +++ b/tests/recipes/recipe-seasonal_downscaling.yml @@ -25,8 +25,8 @@ Analysis: type: 'none' Workflow: Anomalies: - compute: no - cross_validation: + compute: TRUE + cross_validation: FALSE save: 'none' Calibration: method: 'none' @@ -46,8 +46,11 @@ Analysis: target_grid: /esarchive/recon/ecmwf/era5/daily_mean/tas_f1h/tas_199301.nc save: 'all' Output_format: S2S4E + ncores: 1 + remove_NAs: FALSE Run: Loglevel: INFO Terminal: yes output_dir: ./out-logs/ code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/ + diff --git a/tests/testthat/test-seasonal_downscaling.R b/tests/testthat/test-seasonal_downscaling.R index fde0756d4408935d15f6176fa82eebe13624273d..8a52657a8b4dfd3665d82b27af46af6a60375b63 100644 --- a/tests/testthat/test-seasonal_downscaling.R +++ b/tests/testthat/test-seasonal_downscaling.R @@ -5,6 +5,7 @@ source("modules/Loading/Loading.R") source("modules/Skill/Skill.R") source("modules/Downscaling/Downscaling.R") source("modules/Saving/Saving.R") +source("modules/Anomalies/Anomalies.R") recipe_file <- "tests/recipes/recipe-seasonal_downscaling.yml" recipe <- prepare_outputs(recipe_file, disable_checks = F) @@ -14,9 +15,14 @@ suppressWarnings({invisible(capture.output( data <- Loading(recipe) ))}) +# Compute anomalies +suppressWarnings({invisible(capture.output( +ano_data <- Anomalies(recipe, data) +))}) + # Downscale the data suppressWarnings({invisible(capture.output( -downscaled_data <- Downscaling(recipe, data) +downscaled_data <- Downscaling(recipe, ano_data) ))}) # Compute skill metrics @@ -24,6 +30,7 @@ suppressWarnings({invisible(capture.output( skill_metrics <- Skill(recipe, downscaled_data) ))}) +#====================================== test_that("1. Loading", { expect_equal( @@ -101,7 +108,48 @@ as.POSIXct("1994-12-03 11:30:00 UTC", tz = 'UTC') }) #====================================== -test_that("2. Downscaling", { +test_that("2. Anomalies", { + +expect_equal( +is.list(ano_data), +TRUE +) +expect_equal( +names(ano_data), +c("hcst", "obs", "fcst", "hcst.full_val", "obs.full_val") +) +expect_equal( +class(ano_data$hcst), +"s2dv_cube" +) +expect_equal( +class(ano_data$fcst), +"NULL" +) +expect_equal( +dim(ano_data$hcst$data), +c(dat = 1, var = 1, sday = 1, sweek = 1, syear = 4, time = 31, latitude = 4, longitude = 4, ensemble = 25) +) +expect_equal( +mean(ano_data$hcst$data), +-3.466088e-16, +tolerance = 0.0001 +) +expect_equal( +as.vector(drop(ano_data$hcst$data)[3, 26:28, 3, 2, 12]), +c(-1.5654303, -0.6506540, -0.4674155), +tolerance = 0.0001 +) +expect_equal( +range(ano_data$hcst$data), +c(-6.265505, 9.440247), +tolerance = 0.0001 +) + +}) + +#====================================== +test_that("3. Downscaling", { expect_equal( is.list(downscaled_data), @@ -141,24 +189,24 @@ c(ensemble = 1, dat = 1, var = 1, sday = 1, sweek = 1, syear = 4, time = 31, lat ) expect_equal( as.vector(drop(downscaled_data$hcst$data)[1:2,1:2,1,2,3]), -c(293.8468, 293.5830, 294.0713, 293.7104), +c(1.482848, 1.503154, 1.599051, 1.591909), tolerance = 0.0001 ) expect_equal( mean(downscaled_data$hcst$data), -288.5359, +-0.04462092, tolerance = 0.0001 ) expect_equal( range(downscaled_data$hcst$data), -c(284.0950, 297.4405), +c(-4.611440, 5.771648), tolerance = 0.0001 ) }) #====================================== -test_that("3. Metrics", { +test_that("4. Metrics", { expect_equal( is.list(skill_metrics), @@ -182,7 +230,7 @@ dim(skill_metrics$rpss) ) expect_equal( skill_metrics$rpss[1], --0.05325714, +-0.005942857, tolerance = 0.0001 ) expect_equal( @@ -192,7 +240,7 @@ FALSE }) -test_that("4. Check saved data", { +test_that("5. Check saved data", { outputs <- paste0(recipe$Run$output_dir, "/outputs/Downscaling/") expect_equal( @@ -206,3 +254,4 @@ TRUE }) unlink(recipe$Run$output_dir, recursive = TRUE) +