From acd2bc2531cbf5382987e063504cbd1b220b4a96 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Tue, 29 Jun 2021 22:51:51 +0200
Subject: [PATCH 1/7] Let target dimensions be able to have different length
 among dataset

---
 R/ByChunks.R                                  |   3 +-
 .../test-Compute-inconsistent_target_dim.R    | 147 ++++++++++++++++++
 2 files changed, 149 insertions(+), 1 deletion(-)
 create mode 100644 tests/testthat/test-Compute-inconsistent_target_dim.R

diff --git a/R/ByChunks.R b/R/ByChunks.R
index dd10112..0088235 100644
--- a/R/ByChunks.R
+++ b/R/ByChunks.R
@@ -354,7 +354,8 @@ ByChunks <- function(step_fun, cube_headers, ..., chunks = 'auto',
   # Check all input headers have matching dimensions
   cube_index <- 1
   for (cube_header in cube_headers) {
-    if (!all(attr(cube_header, 'Dimensions') == all_dims_merged[names(attr(cube_header, 'Dimensions'))])) {
+    # Only check margin dimensions (i.e., chunked_dims)
+    if (!all(attr(cube_header, 'Dimensions')[chunked_dims] == all_dims_merged[names(attr(cube_header, 'Dimensions'))][chunked_dims])) {
       stop("All provided 'cube_headers' must have matching dimension lengths ",
            "with each other.")
     }
diff --git a/tests/testthat/test-Compute-inconsistent_target_dim.R b/tests/testthat/test-Compute-inconsistent_target_dim.R
new file mode 100644
index 0000000..3da65c0
--- /dev/null
+++ b/tests/testthat/test-Compute-inconsistent_target_dim.R
@@ -0,0 +1,147 @@
+context("Compute()/ByChunks(): dimension consistence check")
+# If dataset are more than 1 (e.g., exp and obs), ByChunks() checks if 
+# they have consistent dimensions in favor of Apply() computation. However,
+# only margin dimensions need to be identical. Target dimensions can have
+# different lengths.
+
+path.exp <- '/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc'
+path.obs <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$var$_$date$.nc'
+
+var <- 'tos'
+y1 <- 1981
+y2 <- 1982
+months <- 12:1
+leadtimemin <- 1; leadtimemax <- 3
+
+indices <- c('NINO34', 'ALT3')
+lons.min <- 190
+lons.max <- 240
+lats.min <- -5
+lats.max <- 5
+
+
+ini_month <- sprintf("%02d", months[1])
+sdate <- paste0(y1:y2, ini_month, '01')
+
+suppressWarnings(
+exp <- Start(data = path.exp,
+             var = var,
+             member = 'all',
+             sdate = sdate,
+             time = leadtimemin:leadtimemax,
+             latitude = values(list(lats.min, lats.max)),
+             latitude_reorder = Sort(decreasing = TRUE),
+             longitude = values(list(lons.min, lons.max)),
+             synonims = list(longitude = c('lon', 'longitude'),
+                             latitude = c('lat', 'latitude'),
+                             member = c('member', 'ensemble')),
+             return_vars = list(longitude = NULL,
+                                latitude = NULL,
+                                time = 'sdate'),
+             retrieve = FALSE)
+)
+lons.exp <- attr(exp, 'Variables')$common$longitude
+lats.exp <- attr(exp, 'Variables')$common$latitude
+dates.exp <- attr(exp, 'Variables')$common$time
+
+# Manually create date/ime
+dates.obs <- c(paste0(y1:y2, ini_month, '01'),
+               paste0((y1 + 1):(y2 + 1), '01', '01'),
+               paste0((y1 + 1):(y2 + 1), '02', '01'))
+time.obs <- as.POSIXct(dates.obs, "%Y%m%d",
+                       origin = "1981-12", tz = 'UTC')
+dim(time.obs) <- c(dim(dates.exp)['sdate'], dim(dates.exp)['time'])
+
+suppressWarnings(
+obs <- Start(data = path.obs,
+             var = var,
+             date = unique(format(time.obs, '%Y%m')),
+             time = values(time.obs),
+             latitude = values(list(lats.min, lats.max)),
+             latitude_reorder = Sort(decreasing = TRUE),
+             longitude = values(list(lons.min, lons.max)),
+             time_across = 'date',
+             merge_across_dims = TRUE,
+             split_multiselected_dims = TRUE,
+             synonims = list(longitude = c('lon', 'longitude'),
+                             latitude = c('lat', 'latitude')),
+             return_vars = list(longitude = NULL,
+                                latitude = NULL,
+                                time = 'date'),
+             retrieve = FALSE)
+)
+
+lons.obs <- attr(obs, 'Variables')$common$longitude
+lats.obs <- attr(obs, 'Variables')$common$latitude
+dates.obs <- attr(obs, 'Variables')$common$time
+
+lons_exp <- as.vector(lons.exp)
+lats_exp <- as.vector(lats.exp)
+lons_obs <- as.vector(lons.obs)
+lats_obs <- as.vector(lats.obs)
+
+
+fun <- function(exp, obs, path.output,  
+                lons_exp = lons_exp, lats_exp = lats_exp,
+                lons_obs = lons_obs, lats_obs = lats_obs) {
+
+        e <- s2dv::MeanDims(drop(exp), c('member', 'time'))
+        sst.e <- ClimProjDiags::WeightedMean(e, lons_exp, lats_exp,
+                              londim = which(names(dim(e)) == 'longitude'),
+                              latdim = which(names(dim(e)) == 'latitude'))
+        index.exp <- (sst.e - mean(sst.e))/sd(sst.e)
+
+        lons.obs <- as.vector(attr(obs, 'Variables')$common$longitude)
+        lats.obs <- as.vector(attr(obs, 'Variables')$common$latitude)
+        o <- s2dv::MeanDims(drop(obs), 'time')
+        sst.o <- ClimProjDiags::WeightedMean(o, lons_obs, lats_obs,
+                              londim = which(names(dim(o)) == 'longitude'),
+                              latdim = which(names(dim(o)) == 'latitude'))
+        index.obs <- (sst.o - mean(sst.o))/sd(sst.o)
+
+        # give dim name
+        dim(index.exp) <- c(sdate = length(index.exp))
+        dim(index.obs) <- c(sdate = length(index.obs))
+
+        return(list(ind_exp = index.exp, ind_obs = index.obs))
+
+}
+
+step <- Step(fun,
+             target_dims = list(exp = c('member', 'sdate', 'time', 'latitude', 'longitude'),
+                                obs = c('sdate', 'time', 'latitude', 'longitude')),
+             output_dims = list(ind_exp = 'sdate', ind_obs = 'sdate'))
+
+workflow <- AddStep(list(exp = exp, obs = obs), step,
+                    lons_exp = lons_exp, lats_exp = lats_exp,
+                    lons_obs = lons_obs, lats_obs = lats_obs)
+
+res <- Compute(workflow$ind_exp,
+               chunks = list(var = 1),
+               threads_load = 2,
+               threads_compute = 4)
+
+
+expect_equal(
+attr(exp, 'Dimensions'),
+c(data = 1, var = 1, member = 25, sdate = 2, time = 3, latitude = 11, longitude = 51)
+)
+expect_equal(
+attr(obs, 'Dimensions'),
+c(data = 1, var = 1, sdate = 2, time = 3, latitude = 41, longitude = 201)
+)
+expect_equal(
+names(res),
+c('ind_exp', 'ind_obs')
+)
+expect_equal(
+mean(res$ind_exp)*10^14,
+-1.376677,
+tolerance = 0.00001
+)
+expect_equal(
+mean(res$ind_obs)*10^14,
+-4.424239,
+tolerance = 0.00001
+)
+
-- 
GitLab


From 0fadf66763b7d27d505c8f1630e03e489ca14dae Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Tue, 29 Jun 2021 23:39:18 +0200
Subject: [PATCH 2/7] Fix unit test

---
 tests/testthat/test-Compute-inconsistent_target_dim.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/testthat/test-Compute-inconsistent_target_dim.R b/tests/testthat/test-Compute-inconsistent_target_dim.R
index 3da65c0..2080be6 100644
--- a/tests/testthat/test-Compute-inconsistent_target_dim.R
+++ b/tests/testthat/test-Compute-inconsistent_target_dim.R
@@ -141,7 +141,7 @@ tolerance = 0.00001
 )
 expect_equal(
 mean(res$ind_obs)*10^14,
--4.424239,
+-1.471046,
 tolerance = 0.00001
 )
 
-- 
GitLab


From 8dac0e7462ba05443623341a51e137c5d7a770bf Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Wed, 30 Jun 2021 01:33:13 +0200
Subject: [PATCH 3/7] Reduce the weight

---
 .../test-Compute-inconsistent_target_dim.R    | 29 ++++++++-----------
 1 file changed, 12 insertions(+), 17 deletions(-)

diff --git a/tests/testthat/test-Compute-inconsistent_target_dim.R b/tests/testthat/test-Compute-inconsistent_target_dim.R
index 2080be6..2e4d6ee 100644
--- a/tests/testthat/test-Compute-inconsistent_target_dim.R
+++ b/tests/testthat/test-Compute-inconsistent_target_dim.R
@@ -9,24 +9,21 @@ path.obs <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$va
 
 var <- 'tos'
 y1 <- 1981
-y2 <- 1982
-months <- 12:1
+y2 <- 1983
 leadtimemin <- 1; leadtimemax <- 3
 
 indices <- c('NINO34', 'ALT3')
-lons.min <- 190
+lons.min <- 220
 lons.max <- 240
 lats.min <- -5
 lats.max <- 5
 
-
-ini_month <- sprintf("%02d", months[1])
-sdate <- paste0(y1:y2, ini_month, '01')
+sdate <- paste0(y1:y2, '1201')
 
 suppressWarnings(
 exp <- Start(data = path.exp,
              var = var,
-             member = 'all',
+             member = indices(1:2),
              sdate = sdate,
              time = leadtimemin:leadtimemax,
              latitude = values(list(lats.min, lats.max)),
@@ -45,7 +42,7 @@ lats.exp <- attr(exp, 'Variables')$common$latitude
 dates.exp <- attr(exp, 'Variables')$common$time
 
 # Manually create date/ime
-dates.obs <- c(paste0(y1:y2, ini_month, '01'),
+dates.obs <- c(paste0(y1:y2, '1201'),
                paste0((y1 + 1):(y2 + 1), '01', '01'),
                paste0((y1 + 1):(y2 + 1), '02', '01'))
 time.obs <- as.POSIXct(dates.obs, "%Y%m%d",
@@ -117,31 +114,29 @@ workflow <- AddStep(list(exp = exp, obs = obs), step,
                     lons_obs = lons_obs, lats_obs = lats_obs)
 
 res <- Compute(workflow$ind_exp,
-               chunks = list(var = 1),
-               threads_load = 2,
-               threads_compute = 4)
+               chunks = list(var = 1))
 
 
 expect_equal(
 attr(exp, 'Dimensions'),
-c(data = 1, var = 1, member = 25, sdate = 2, time = 3, latitude = 11, longitude = 51)
+c(data = 1, var = 1, member = 2, sdate = 3, time = 3, latitude = 11, longitude = 21)
 )
 expect_equal(
 attr(obs, 'Dimensions'),
-c(data = 1, var = 1, sdate = 2, time = 3, latitude = 41, longitude = 201)
+c(data = 1, var = 1, sdate = 3, time = 3, latitude = 41, longitude = 81)
 )
 expect_equal(
 names(res),
 c('ind_exp', 'ind_obs')
 )
 expect_equal(
-mean(res$ind_exp)*10^14,
--1.376677,
+mean(res$ind_exp)*10^18,
+-9.251859,
 tolerance = 0.00001
 )
 expect_equal(
-mean(res$ind_obs)*10^14,
--1.471046,
+mean(res$ind_obs)*10^15,
+-9.584944,
 tolerance = 0.00001
 )
 
-- 
GitLab


From 648d807a188ccde4b27097027c527782d226ea4a Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 1 Jul 2021 13:52:18 +0200
Subject: [PATCH 4/7] Create use case for different lengths of target dim.
 AAAdd check for inconsistent margin dim

---
 R/ByChunks.R                                  |   9 ++
 .../ex2_11_two_dat_inconsistent_target_dim.R  | 149 ++++++++++++++++++
 .../test-Compute-inconsistent_target_dim.R    |   6 +-
 3 files changed, 159 insertions(+), 5 deletions(-)
 create mode 100644 inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R

diff --git a/R/ByChunks.R b/R/ByChunks.R
index 0088235..8185763 100644
--- a/R/ByChunks.R
+++ b/R/ByChunks.R
@@ -354,6 +354,15 @@ ByChunks <- function(step_fun, cube_headers, ..., chunks = 'auto',
   # Check all input headers have matching dimensions
   cube_index <- 1
   for (cube_header in cube_headers) {
+
+    # Check if all the margin dims are consistent among datasets
+    if (!all(chunked_dims %in% names(attr(cube_header, "Dimensions")))) {
+      trouble_dim_name <- chunked_dims[which(!chunked_dims %in%              
+                                             names(attr(cube_header, "Dimensions")))]
+      stop(paste0("Found margin dimension, ", toString(trouble_dim_name),
+           ", is not in input data ", cube_index, "."))
+    }
+
     # Only check margin dimensions (i.e., chunked_dims)
     if (!all(attr(cube_header, 'Dimensions')[chunked_dims] == all_dims_merged[names(attr(cube_header, 'Dimensions'))][chunked_dims])) {
       stop("All provided 'cube_headers' must have matching dimension lengths ",
diff --git a/inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R b/inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R
new file mode 100644
index 0000000..521a248
--- /dev/null
+++ b/inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R
@@ -0,0 +1,149 @@
+# Author: Chihchung Chou, An-Chi Ho
+# Date: 1st July 2021
+# ------------------------------------------------------------------
+# This use case uses experimental and the corresponding observational data to calculate 
+# the temporal mean and spatial weighted mean. Notice that the spatial resolutions of the
+# two datasets are different, but it still works because lat and lon are target dimensions.
+# ------------------------------------------------------------------
+library(startR)
+
+# exp
+  repos <- paste0('/esarchive/exp/ecmwf/system4_m1/6hourly/',
+                  '$var$/$var$_$sdate$.nc')
+  sdate <- sapply(1994, function(x) paste0(x, sprintf('%02d', 9:12), '01'))
+
+
+  exp <- Start(dat = repos,
+               var = 'tas',
+               sdate = sdate,
+               time = indices(seq(1, 124, 4)),  #first time step per day
+               ensemble = 'all',
+               latitude = values(list(lats.min, lats.max)),
+               latitude_reorder = Sort(),
+               longitude = values(list(lons.min, lons.max)),
+               longitude_reorder = CircularSort(0, 360),
+               return_vars = list(latitude = NULL,
+                                  longitude = NULL,
+                                  time = 'sdate'))
+
+
+  dates <- attr(system4, 'Variables')$common$time
+
+  dates_file <- sort(unique(gsub('-', '', sapply(as.character(dates),
+  substr, 1, 7))))
+
+
+repos_exp <- paste0('/esarchive/exp/ecearth/a1tr/cmorfiles/CMIP/EC-Earth-Consortium/',
+                    'EC-Earth3/historical/r24i1p1f1/Amon/$var$/gr/v20190312/',
+                    '$var$_Amon_EC-Earth3_historical_r24i1p1f1_gr_$sdate$01-$sdate$12.nc')
+
+exp <- Start(dat = repos_exp,
+             var = 'tas',
+             sdate = as.character(c(2005:2008)),
+             time = indices(1:3),
+             lat = 'all',
+             lon = 'all',
+             synonims = list(lat = c('lat', 'latitude'),
+                             lon = c('lon', 'longitude')),
+             return_vars = list(lon = NULL,
+                                lat = NULL,
+                                time = 'sdate'),
+             retrieve = FALSE)
+
+lons_exp <- as.vector(attr(exp, 'Variables')$common$lon)
+lats_exp <- as.vector(attr(exp, 'Variables')$common$lat)
+dates_exp <- attr(exp, 'Variables')$common$time
+
+attr(exp, 'Dimensions')
+#  dat   var sdate  time   lat   lon 
+#    1     1     4     3   256   512 
+dates_exp
+# [1] "2005-01-16 12:00:00 UTC" "2006-01-16 12:00:00 UTC"
+# [3] "2007-01-16 12:00:00 UTC" "2008-01-16 12:00:00 UTC"
+# [5] "2005-02-15 00:00:00 UTC" "2006-02-15 00:00:00 UTC"
+# [7] "2007-02-15 00:00:00 UTC" "2008-02-15 12:00:00 UTC"
+# [9] "2005-03-16 12:00:00 UTC" "2006-03-16 12:00:00 UTC"
+#[11] "2007-03-16 12:00:00 UTC" "2008-03-16 12:00:00 UTC"
+
+
+# obs
+  path.obs <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$var$_$date$.nc'
+  obs <- Start(dat = path.obs,
+               var = 'tas',
+               date = unique(format(dates, '%Y%m')),
+               time = values(dates), 
+               time_across = 'date',
+               merge_across_dims = TRUE,
+               split_multiselected_dims = TRUE,
+               lat = 'all',
+               lon = 'all',
+               synonims = list(lon = c('lon', 'longitude'),
+                               lat = c('lat', 'latitude')),
+               return_vars = list(lon = NULL,
+                                  lat = NULL,
+                                  time = 'date'),
+               retrieve = FALSE)
+
+lons_obs <- as.vector(attr(obs, 'Variables')$common$lon)
+lats_obs <- as.vector(attr(obs, 'Variables')$common$lat)
+dates_obs <- attr(obs, 'Variables')$common$time
+
+
+attr(obs, 'Dimensions')
+#     data       var     sdate      time     lat      lon 
+#        1         1         4         3     721     1440 
+dates_obs
+# [1] "2005-01-16 11:30:00 UTC" "2006-01-16 11:30:00 UTC"
+# [3] "2007-01-16 11:30:00 UTC" "2008-01-16 11:30:00 UTC"
+# [5] "2005-02-14 23:30:00 UTC" "2006-02-14 23:30:00 UTC"
+# [7] "2007-02-14 23:30:00 UTC" "2008-02-15 11:30:00 UTC"
+# [9] "2005-03-16 11:30:00 UTC" "2006-03-16 11:30:00 UTC"
+#[11] "2007-03-16 11:30:00 UTC" "2008-03-16 11:30:00 UTC"
+
+
+fun <- function(exp, obs, 
+                lons_exp = lons_exp, lats_exp = lats_exp,
+                lons_obs = lons_obs, lats_obs = lats_obs) {
+  # exp
+  e <- s2dv::MeanDims(drop(exp), 'time')
+  sst.e <- ClimProjDiags::WeightedMean(e, lons_exp, lats_exp,
+                                       londim = which(names(dim(e)) == 'lon'),
+                                       latdim = which(names(dim(e)) == 'lat'))
+  index.exp <- (sst.e - mean(sst.e))/sd(sst.e)
+  
+  # obs         
+  o <- s2dv::MeanDims(drop(obs), 'time')
+  sst.o <- ClimProjDiags::WeightedMean(o, lons_obs, lats_obs,
+                                       londim = which(names(dim(o)) == 'lon'),
+                                       latdim = which(names(dim(o)) == 'lat'))
+  index.obs <- (sst.o - mean(sst.o))/sd(sst.o)
+
+  # give dim name
+  dim(index.exp) <- c(sdate = length(index.exp))
+  dim(index.obs) <- c(sdate = length(index.obs))
+
+  return(list(ind_exp = index.exp, ind_obs = index.obs))
+
+}
+
+# If ClimProjDiags::WeightedMean accepts two-dim input, 'sdate' can be margin dimension.
+step <- Step(fun,
+             target_dims = list(exp = c('sdate', 'time', 'lat', 'lon'),
+                                obs = c('sdate', 'time', 'lat', 'lon')),
+             output_dims = list(ind_exp = 'sdate', ind_obs = 'sdate'))
+
+workflow <- AddStep(list(exp = exp, obs = obs), step,
+                    lons_exp = lons_exp, lats_exp = lats_exp,
+                    lons_obs = lons_obs, lats_obs = lats_obs)
+
+res <- Compute(workflow$ind_exp,
+               chunks = list(var = 1))
+
+str(res)
+#List of 2
+# $ ind_exp: num [1:4, 1, 1] 1.195 0.422 -0.6 -1.017
+# $ ind_obs: num [1:4, 1, 1] 0.4642 0.0206 0.9123 -1.3971
+# ...
+# ...
+
+
diff --git a/tests/testthat/test-Compute-inconsistent_target_dim.R b/tests/testthat/test-Compute-inconsistent_target_dim.R
index 2e4d6ee..241ab6e 100644
--- a/tests/testthat/test-Compute-inconsistent_target_dim.R
+++ b/tests/testthat/test-Compute-inconsistent_target_dim.R
@@ -10,9 +10,7 @@ path.obs <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$va
 var <- 'tos'
 y1 <- 1981
 y2 <- 1983
-leadtimemin <- 1; leadtimemax <- 3
 
-indices <- c('NINO34', 'ALT3')
 lons.min <- 220
 lons.max <- 240
 lats.min <- -5
@@ -25,7 +23,7 @@ exp <- Start(data = path.exp,
              var = var,
              member = indices(1:2),
              sdate = sdate,
-             time = leadtimemin:leadtimemax,
+             time = 1:3,
              latitude = values(list(lats.min, lats.max)),
              latitude_reorder = Sort(decreasing = TRUE),
              longitude = values(list(lons.min, lons.max)),
@@ -88,8 +86,6 @@ fun <- function(exp, obs, path.output,
                               latdim = which(names(dim(e)) == 'latitude'))
         index.exp <- (sst.e - mean(sst.e))/sd(sst.e)
 
-        lons.obs <- as.vector(attr(obs, 'Variables')$common$longitude)
-        lats.obs <- as.vector(attr(obs, 'Variables')$common$latitude)
         o <- s2dv::MeanDims(drop(obs), 'time')
         sst.o <- ClimProjDiags::WeightedMean(o, lons_obs, lats_obs,
                               londim = which(names(dim(o)) == 'longitude'),
-- 
GitLab


From b18a1cecd2ced0e5dc184329c7ce07339a1b14e5 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 1 Jul 2021 14:10:45 +0200
Subject: [PATCH 5/7] Add new usecase to the document

---
 inst/doc/usecase.md | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/inst/doc/usecase.md b/inst/doc/usecase.md
index 06f2c04..2d1ab9d 100644
--- a/inst/doc/usecase.md
+++ b/inst/doc/usecase.md
@@ -72,8 +72,9 @@ You can find more explanation in FAQ [How-to-20](inst/doc/faq.md#20-use-metadata
    10. [Apply an existing mask on data](inst/doc/usecase/ex2_10_existing_mask.R)  
       This use case shows you how to apply the existing mask file on your data.
 If you need to create the mask file on your own, go to ex2_9_mask.R.
-
-
-
+   11. [Two datasets with different length of target dimensions](inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R)
+      This use case uses experimental and the corresponding observational data to calculate 
+the temporal mean and spatial weighted mean. Notice that the spatial resolutions of the two 
+datasets are different, but it still works because lat and lon are target dimensions.
 
  
-- 
GitLab


From 098ad85e7dfb7ae5a7395167b9fe4ccd435d6729 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 1 Jul 2021 16:28:54 +0200
Subject: [PATCH 6/7] Correct use case

---
 inst/doc/usecase/ex2_4_two_func.R             | 27 +++++++++----------
 inst/doc/usecase/ex2_6_ext_param_func.R       |  8 +++---
 .../ex2_7_seasonal_forecast_verification.R    |  4 +--
 inst/doc/usecase/ex2_9_mask.R                 | 18 +++++++++++--
 4 files changed, 34 insertions(+), 23 deletions(-)

diff --git a/inst/doc/usecase/ex2_4_two_func.R b/inst/doc/usecase/ex2_4_two_func.R
index 4f3da94..d4d1b5f 100644
--- a/inst/doc/usecase/ex2_4_two_func.R
+++ b/inst/doc/usecase/ex2_4_two_func.R
@@ -14,38 +14,35 @@ library(startR)
                 retrieve = FALSE)
   
   fun_deb3 <- function(x) {
-               source("/esarchive/scratch/nperez/Season_v2.R")
                lons_data = as.vector(attr(x, 'Variables')$dat1$longitude)
                lats_data = as.vector(attr(x, 'Variables')$dat1$latitude)
-               resgrid = "r360x180"   # prlr
-               y = Season_v2(x, posdim = 'time', monini = 1, moninf = 1, monsup = 3)
-               r <- s2dverification::CDORemap(y, lons_data, lats_data, resgrid,
-                                        'bil', crop = FALSE, force_remap = TRUE)[[1]]
+               resgrid = "r360x180"
+#               y <- s2dv::Season(x, time_dim = 'time', monini = 1, moninf = 1, monsup = 3)
+               y <- apply(x, c(1, 2), s2dv:::.Season, monini = 1, moninf = 1, monsup = 3)
+               r <- s2dv::CDORemap(y, lons_data, lats_data, resgrid,
+                                  'bil', crop = FALSE, force_remap = TRUE)[[1]]
                return(r)
   }
   
   step4 <- Step(fun = fun_deb3,
                target_dims = c('latitude','longitude', 'time'),
-               output_dims = c('latitude', 'longitude', 'time'),
+               output_dims = c('latitude', 'longitude'),
                use_attributes = list(data = "Variables"))
-  wf4 <- AddStep(data, step4)
+  wf4 <- AddStep(list(data = data), step4)
   
 ## locally
   res4 <- Compute(workflow = wf4,
               chunks = list(ensemble = 2,
                             sdate = 2))
-  dim(res4$output1)
-  head(res4$output1)
-  summary(res4$output1)
 
 # ------------------------------------------------------------------
 # Output:
-#> dim(res4$output1)
-# latitude longitude      time       dat       var     sdate  ensemble
-#      180       360         1         1         1         2         2
-#> head(res4$output1)
+dim(res4$output1)
+# latitude longitude       dat       var     sdate  ensemble
+#      180       360         1         1         2         2
+head(res4$output1)
 #[1] 237.1389 237.2601 238.0882 238.0312 237.7883 238.4835
-#> summary(res4$output1)
+summary(res4$output1)
 #   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 #  227.3   259.6   280.8   277.1   296.2   306.7
 # ------------------------------------------------------------------
diff --git a/inst/doc/usecase/ex2_6_ext_param_func.R b/inst/doc/usecase/ex2_6_ext_param_func.R
index 302b52a..8bd128f 100644
--- a/inst/doc/usecase/ex2_6_ext_param_func.R
+++ b/inst/doc/usecase/ex2_6_ext_param_func.R
@@ -88,15 +88,15 @@
 
 # Notice that the function uses rnorm() inside. So the results will be different.
 # -----------------------------------------------------------
-#names(res)
+names(res)
 #[1] "strat"  "t_test"
-#> dim(res$strat)
+dim(res$strat)
 #    phase longitude  latitude 
 #        8        30        20 
-#> summary(res$strat)
+summary(res$strat)
 #     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
 #-0.133300 -0.032530 -0.001822 -0.005715  0.031700  0.094220 
-#> res$strat[1:5, 1:2, 1]
+res$strat[1:5, 1:2, 1]
 #            [,1]        [,2]
 #[1,] -0.04661354 -0.04661539
 #[2,] -0.01058483 -0.01053589
diff --git a/inst/doc/usecase/ex2_7_seasonal_forecast_verification.R b/inst/doc/usecase/ex2_7_seasonal_forecast_verification.R
index f24fa42..8b74607 100644
--- a/inst/doc/usecase/ex2_7_seasonal_forecast_verification.R
+++ b/inst/doc/usecase/ex2_7_seasonal_forecast_verification.R
@@ -33,7 +33,7 @@
 
 
   func <- function(x, y) {
-    crps <- mean(SpecsVerification::EnsCrps(x, y, R.new = Inf))
+    crps <- mean(SpecsVerification::EnsCrps(x, y, R.new = Inf), na.rm = T)
     return(crps)
   }
   step <- Step(func, target_dims = list(c('sdate', 'ensemble'), c('sdate')),
@@ -100,7 +100,7 @@
 #        1         1         1       256       512 
  summary(res$output1)
 #   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-#0.09941 0.37760 0.71640 0.83570 1.20300 6.23400 
+#0.09882 0.37815 0.71648 0.83638 1.20353 6.23452 
 
 
 # Plotting
diff --git a/inst/doc/usecase/ex2_9_mask.R b/inst/doc/usecase/ex2_9_mask.R
index 13be878..aca6162 100644
--- a/inst/doc/usecase/ex2_9_mask.R
+++ b/inst/doc/usecase/ex2_9_mask.R
@@ -107,15 +107,20 @@ wf_mask <- AddStep(list(data, mask), stepMask)
 res <- Compute(workflow = wf_mask,
                 chunks = list(latitude = 2,
                               longitude = 2))
+
 #####################################################################
 # --------------------------------------------------------------------------------
 #####################################################################
 # Extra lines for output verification:
 # Output verification:
-
-summary(res$output1)
 dim(res$output1)
+#      var  latitude longitude 
+#        1        42        40 
+summary(res$output1)
+#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
+#  295.4   300.0   300.2   300.0   300.5   301.0     840 
 head(res$output1)
+#[1] 300.0644       NA       NA       NA       NA       NA
 
 mask_loaded <-  Start(dat = path,
                       var = 'mask',
@@ -125,8 +130,13 @@ mask_loaded <-  Start(dat = path,
                                          longitude = 'dat'),
                       retrieve = TRUE)
 summary(res$output1[mask_loaded == 0]) # All are NA's: correct
+#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
+#     NA      NA      NA     NaN      NA      NA     840 
 summary(res$output1[mask_loaded == 1]) # There is no NA's: correct
+#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
+#  295.4   300.0   300.2   300.0   300.5   301.0 
 sum(mask_loaded == 0) # The number of NA's are 840: correct
+#[1] 840
 
 # compare values:
 data_loaded <- Start(dat = repos,
@@ -141,9 +151,13 @@ data_loaded <- Start(dat = repos,
                                         time = 'sdate'),
                      retrieve = TRUE)
 mean(data_loaded[1,1, , , ,1,1])
+#[1] 300.0644
 res$output1[1,1,1]
+#[1] 300.0644
 mean(data_loaded[1,1, , , ,1,2])
+#[1] 300.3169
 res$output1[1,1,2]
+#[1] 300.3169
 
 out <- mask_loaded
 for (i in 1:(dim(data_loaded)['latitude'])) {
-- 
GitLab


From eb429e4299af1a803a5dea309236a28da6135235 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Mon, 5 Jul 2021 10:34:21 +0200
Subject: [PATCH 7/7] Correct the script

---
 .../ex2_11_two_dat_inconsistent_target_dim.R  | 30 ++-----------------
 1 file changed, 3 insertions(+), 27 deletions(-)

diff --git a/inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R b/inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R
index 521a248..69267b0 100644
--- a/inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R
+++ b/inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R
@@ -7,32 +7,8 @@
 # ------------------------------------------------------------------
 library(startR)
 
-# exp
-  repos <- paste0('/esarchive/exp/ecmwf/system4_m1/6hourly/',
-                  '$var$/$var$_$sdate$.nc')
-  sdate <- sapply(1994, function(x) paste0(x, sprintf('%02d', 9:12), '01'))
-
-
-  exp <- Start(dat = repos,
-               var = 'tas',
-               sdate = sdate,
-               time = indices(seq(1, 124, 4)),  #first time step per day
-               ensemble = 'all',
-               latitude = values(list(lats.min, lats.max)),
-               latitude_reorder = Sort(),
-               longitude = values(list(lons.min, lons.max)),
-               longitude_reorder = CircularSort(0, 360),
-               return_vars = list(latitude = NULL,
-                                  longitude = NULL,
-                                  time = 'sdate'))
-
-
-  dates <- attr(system4, 'Variables')$common$time
-
-  dates_file <- sort(unique(gsub('-', '', sapply(as.character(dates),
-  substr, 1, 7))))
-
 
+# exp
 repos_exp <- paste0('/esarchive/exp/ecearth/a1tr/cmorfiles/CMIP/EC-Earth-Consortium/',
                     'EC-Earth3/historical/r24i1p1f1/Amon/$var$/gr/v20190312/',
                     '$var$_Amon_EC-Earth3_historical_r24i1p1f1_gr_$sdate$01-$sdate$12.nc')
@@ -70,8 +46,8 @@ dates_exp
   path.obs <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h-r1440x721cds/$var$_$date$.nc'
   obs <- Start(dat = path.obs,
                var = 'tas',
-               date = unique(format(dates, '%Y%m')),
-               time = values(dates), 
+               date = unique(format(dates_exp, '%Y%m')),
+               time = values(dates_exp), 
                time_across = 'date',
                merge_across_dims = TRUE,
                split_multiselected_dims = TRUE,
-- 
GitLab