From 0c534fc3c8c2d74c08cb55007229b137152056ae Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 27 Jan 2021 14:49:28 +0100 Subject: [PATCH 1/3] Correct the time attribute calculation when units = month. --- R/NcDataReader.R | 72 +++++++++++++++++++++++++--- tests/testthat/test-Start-calendar.R | 30 +++++++++++- 2 files changed, 95 insertions(+), 7 deletions(-) diff --git a/R/NcDataReader.R b/R/NcDataReader.R index 222d11d..72abc09 100644 --- a/R/NcDataReader.R +++ b/R/NcDataReader.R @@ -197,6 +197,10 @@ NcDataReader <- function(file_path = NULL, file_object = NULL, result[] <- paste(result[], units) } else if (grepl(' since ', units)) { + # Find the calendar + calendar <- attr(result, 'variables')$time$calendar + if (calendar == 'standard') calendar <- 'gregorian' + parts <- strsplit(units, ' since ')[[1]] units <- parts[1] @@ -211,14 +215,70 @@ NcDataReader <- function(file_path = NULL, file_object = NULL, # units <- 'days' result <- result * 24 * 60 * 60 # day to sec } else if (units %in% c('month', 'months')) { - result <- result * 30.5 - result <- result * 24 * 60 * 60 # day to sec -# units <- 'days' + # define day in each month + leap_month_day <- c(31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) + no_leap_month_day <- c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) + # Origin year and month + ori_year <- as.numeric(substr(parts[2], 1, 4)) + ori_month <- as.numeric(substr(parts[2], 6, 7)) + + if (calendar == 'gregorian') { + # Find how many years + months + yr_num <- floor(result / 12) + month_left <- result - yr_num * 12 + # Find the leap years we care + if (ori_month <= 2) { + leap_num <- length(which(sapply(ori_year:(ori_year + yr_num - 1), s2dv::LeapYear))) + } else { + leap_num <- length(which(sapply((ori_year + 1):(ori_year + yr_num), s2dv::LeapYear))) + } + total_days <- leap_num * 366 + (yr_num - leap_num) * 365 # not include month_left yet + + + if (month_left != 0) { + if ((ori_month + month_left) <= 12) { # the last month is still in the same last yr + # Is the last year a leap year? + last_leap <- s2dv::LeapYear(ori_year + yr_num) + if (last_leap) { + total_days <- total_days + sum(leap_month_day[ori_month:(ori_month + month_left - 1)]) + } else { + total_days <- total_days + sum(no_leap_month_day[ori_month:(ori_month + month_left - 1)]) + } + } else { # the last month ends in the next yr + if (ori_month == 2) { # e.g., 2005-02-16 + 11mth = 2006-01-16 + last_leap <- s2dv::LeapYear(ori_year + yr_num) # still consider 2005 + if (last_leap) { + total_days <- total_days + sum(leap_month_day[2:12]) + } else { + total_days <- total_days + sum(no_leap_month_day[2:12]) + } + } else { # e.g., 2005-04-16 + 11mth = 2006-03-16 + last_leap <- s2dv::LeapYear(ori_year + yr_num + 1) + needed_month <- c(ori_month:12, 1:(ori_month + month_left - 12 - 1)) + if (last_leap) { + total_days <- total_days + sum(leap_month_day[needed_month]) + } else { + total_days <- total_days + sum(no_leap_month_day[needed_month]) + } + } + } + } + result <- total_days * 24 * 60 * 60 # day to sec + } else if (calendar %in% c('365_day',' 365', 'noleap')) { + yr_num <- floor(result / 12) + month_left <- result - yr_num * 12 + total_days <- 365 * yr_num + sum(no_leap_month_day[ori_month:(month_left - 1)]) + result <- total_days * 24 * 60 * 60 # day to sec + + } else if (calendar %in% c('360_day', '360')) { + result <- result * 30 * 24 * 60 * 60 # day to sec + + } else { #old code. The calendar is not in any of the above. + result <- result * 30.5 + result <- result * 24 * 60 * 60 # day to sec + } } - # Find the calendar - calendar <- attr(result, 'variables')$time$calendar - if (calendar == 'standard') calendar <- 'gregorian' new_array <- PCICt::as.PCICt(result, cal = calendar, origin = parts[2])[] new_array <- suppressWarnings(PCICt::as.POSIXct.PCICt(new_array, tz = "UTC")) diff --git a/tests/testthat/test-Start-calendar.R b/tests/testthat/test-Start-calendar.R index 19ac0df..654f8b8 100644 --- a/tests/testthat/test-Start-calendar.R +++ b/tests/testthat/test-Start-calendar.R @@ -147,7 +147,7 @@ expect_equal( }) -test_that("4. gregorian/standard, 6hrly", { +test_that("5. gregorian/standard, 6hrly", { repos_obs <- paste0('/esarchive/recon/ecmwf/erainterim/6hourly/', '$var$/$var$_199405.nc') date <- paste0('1994-05-', sprintf('%02d', 1:31), ' 00:00:00') @@ -178,3 +178,31 @@ expect_equal( tolerance = 0.0001 ) }) + +test_that("6. gregorian/standard, monthly, unit = month", { + + repos_obs <- '/esarchive/obs/ukmo/hadisst_v1.1/monthly_mean/$var$/$var$_$date$.nc' + obs <- Start(dat = repos_obs, + var = 'tos', + date = '200505', #dates_file, + time = 'all', + lat = indices(1:10), + lon = indices(1:10), + time_across = 'date', + #combine time and file_date dims + merge_across_dims = TRUE, + #exclude the additional NAs generated by merge_across_dims + merge_across_dims_narm = TRUE, + synonims = list(lat = c('lat', 'latitude'), + lon = c('lon', 'longitude')), + return_vars = list(lat = 'dat', + lon = 'dat', + time = 'date'), + retrieve = FALSE) + +expect_equal( + attr(obs, 'Variables')$common$time + as.POSIXct('2000-05-16 00:00:00', tz = 'UTC') +) + +}) -- GitLab From 829c26805f3181944421c95913b96b665583036a Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 27 Jan 2021 15:16:26 +0100 Subject: [PATCH 2/3] Correct new unit test. --- tests/testthat/test-Start-calendar.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-Start-calendar.R b/tests/testthat/test-Start-calendar.R index 654f8b8..3a667bd 100644 --- a/tests/testthat/test-Start-calendar.R +++ b/tests/testthat/test-Start-calendar.R @@ -201,8 +201,8 @@ test_that("6. gregorian/standard, monthly, unit = month", { retrieve = FALSE) expect_equal( - attr(obs, 'Variables')$common$time - as.POSIXct('2000-05-16 00:00:00', tz = 'UTC') + attr(obs, 'Variables')$common$time[1, 1], + as.POSIXct('2005-05-16 12:00:00', tz = 'UTC') ) }) -- GitLab From 2be260e9cd7c43ee080c2a35a6a2b4edaaaa01c4 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 27 Jan 2021 18:13:48 +0100 Subject: [PATCH 3/3] Add more unit tests. --- tests/testthat/test-Start-calendar.R | 75 +++++++++++++++++++++++++--- 1 file changed, 69 insertions(+), 6 deletions(-) diff --git a/tests/testthat/test-Start-calendar.R b/tests/testthat/test-Start-calendar.R index 3a667bd..328aa0e 100644 --- a/tests/testthat/test-Start-calendar.R +++ b/tests/testthat/test-Start-calendar.R @@ -1,6 +1,6 @@ context("Start() different calendar") -test_that("1. 360_day, daily", { +test_that("1. 360_day, daily, unit = 'days since 1850-01-01'", { path_hadgem3 <- paste0('/esarchive/exp/CMIP6/dcppA-hindcast/hadgem3-gc31-mm/', 'cmip6-dcppA-hindcast_i1p1/DCPP/MOHC/HadGEM3-GC31-MM/', @@ -47,7 +47,7 @@ expect_equal( }) -test_that("2. 365_day, daily", { +test_that("2. 365_day, daily, unit = 'days since 1984-01-01'", { path_bcc_csm2 <- '/esarchive/exp/CMIP6/dcppA-hindcast/bcc-csm2-mr/cmip6-dcppA-hindcast_i1p1/DCPP/BCC/BCC-CSM2-MR/dcppA-hindcast/r1i1p1f1/day/$var$/gn/v20200408/$var$_day_BCC-CSM2-MR_dcppA-hindcast_s$sdate$-r1i1p1f1_gn_19800101-19891231.nc' data <- Start(dat = path_bcc_csm2, @@ -76,7 +76,7 @@ expect_equal( }) -test_that("3. gregorian/standard, daily", { +test_that("3. standard, daily, unit = 'days since 1850-1-1 00:00:00'", { path_mpi_esm <- paste0('/esarchive/exp/CMIP6/dcppA-hindcast/mpi-esm1-2-hr/', 'cmip6-dcppA-hindcast_i1p1/DCPP/MPI-M/MPI-ESM1-2-HR/', @@ -114,7 +114,7 @@ expect_equal( }) -test_that("4. gregorian/standard, monthly", { +test_that("4. standard, monthly, unit = 'days since 1850-1-1 00:00:00'", { path_mpi_esm <- paste0('/esarchive/exp/CMIP6/dcppA-hindcast/mpi-esm1-2-hr/', 'cmip6-dcppA-hindcast_i1p1/DCPP/MPI-M/MPI-ESM1-2-HR/', @@ -147,7 +147,7 @@ expect_equal( }) -test_that("5. gregorian/standard, 6hrly", { +test_that("5. proleptic_gregorian, 6hrly, unit = 'hours since 2000-11-01 00:00:00'", { repos_obs <- paste0('/esarchive/recon/ecmwf/erainterim/6hourly/', '$var$/$var$_199405.nc') date <- paste0('1994-05-', sprintf('%02d', 1:31), ' 00:00:00') @@ -179,7 +179,7 @@ expect_equal( ) }) -test_that("6. gregorian/standard, monthly, unit = month", { +test_that("6. standard, monthly, unit = 'months since 1870-01-16 12:00:00'", { repos_obs <- '/esarchive/obs/ukmo/hadisst_v1.1/monthly_mean/$var$/$var$_$date$.nc' obs <- Start(dat = repos_obs, @@ -206,3 +206,66 @@ expect_equal( ) }) + +test_that("7. proleptic_gregorian, monthly, unit = 'days since 1850-1-1 00:00:00'", { + + repos <- '/esarchive/exp/mpi-esm-lr/cmip5-historical_i0p1/monthly_mean/$var$/$var$_$sdate$.nc' + data <- Start(dat = repos, + var = 'tas', + sdate = '20000101', + time = indices(1:3), + ensemble = indices(1), + latitude = indices(1:4), + longitude = indices(1:3), + return_vars = list(time = NULL)) + + time <- c(as.POSIXct('2000-01-16 12:00:00', tz = 'UTC'), + as.POSIXct('2000-02-15 12:00:00', tz = 'UTC'), + as.POSIXct('2000-03-16 12:00:00', tz = 'UTC')) + attr(time, "tzone") <- "UTC" + + expect_equal( + attr(data, 'Variables')$common$time[1], + time[1] + ) + expect_equal( + attr(data, 'Variables')$common$time[2], + time[2] + ) + expect_equal( + attr(data, 'Variables')$common$time[3], + time[3] + ) + +}) + +test_that("8. gregorian, 3hrly, unit = 'days since 1850-1-1'", { + repos <- '/esarchive/exp/CMIP5/historical/ecearth/cmip5-historical_i0p1/$var$_3hr_EC-EARTH_historical_r6i1p1_$period$.nc' + data <- Start(dat = repos, + var = 'vas', + period = '200501010300-200601010000', + time = indices(1:3), + lat = indices(1:4), + lon = indices(1:3), + return_vars = list(time = NULL)) + + time <- c(as.POSIXct('2005-01-01 03:00:00', tz = 'UTC'), + as.POSIXct('2005-01-01 06:00:00', tz = 'UTC'), + as.POSIXct('2005-01-01 09:00:00', tz = 'UTC')) + attr(time, "tzone") <- "UTC" + + expect_equal( + attr(data, 'Variables')$common$time[1], + time[1] + ) + expect_equal( + attr(data, 'Variables')$common$time[2], + time[2] + ) + expect_equal( + attr(data, 'Variables')$common$time[3], + time[3] + ) + +}) + -- GitLab