From 05f8c82432f8f6219c9f63d341f3ce657243327f Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Wed, 16 Nov 2022 13:36:15 +0100
Subject: [PATCH] Reorder the largest other dim to the last dim as unlimited
 dim

---
 R/CDORemap.R                   | 12 +++++++---
 tests/testthat/test-CDORemap.R | 44 ++++++++++++++++++++++++++++++++++
 2 files changed, 53 insertions(+), 3 deletions(-)

diff --git a/R/CDORemap.R b/R/CDORemap.R
index f0044cb..927b107 100644
--- a/R/CDORemap.R
+++ b/R/CDORemap.R
@@ -666,7 +666,10 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
     total_slices <- 1
     other_dims_per_chunk <- ifelse(is_irregular, 1, 2)  # 4 (the maximum accepted by CDO) - 2 (lon, lat) = 2.
     if (length(other_dims) > 1 || (length(other_dims) > 0 && (is_irregular))) {
-      if (!(length(dim(data_array)) %in% other_dims)) {
+      # If lat/lon is the last dimension OR the largest other_dims is not the last one, 
+      # reorder the largest other dimension to the last as unlimited dim.
+      if (!(length(dim(data_array)) %in% other_dims) |
+          which.max(dim(data_array)[other_dims]) != length(other_dims)) {
         if (avoid_writes || is_irregular) {
           dims_mod <- dim(data_array)
           dims_mod[which(names(dim(data_array)) %in%
@@ -675,9 +678,9 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
           permutation <- (1:length(dim(data_array)))[-dim_to_move]
           permutation <- c(permutation, dim_to_move)
           permutation_back <- sort(permutation, index.return = TRUE)$ix
-          dim_backup <- dim(data_array)
+#          dim_backup <- dim(data_array)
           data_array <- aperm(data_array, permutation)
-          dim(data_array) <- dim_backup[permutation]
+#          dim(data_array) <- dim_backup[permutation]
           other_dims <- which(!(names(dim(data_array)) %in% c(lon_dim, lat_dim)))
         } else {
           # We allow only lon, lat and 1 more dimension per chunk, so 
@@ -695,6 +698,9 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
         total_slices <- prod(dim(slices_to_iterate))
       }
       if ((other_dims_per_chunk > 1) || (other_dims_per_chunk > 0 && is_irregular)) {
+        #NOTE: Why don't we use the second line here? In history, that line was never used.
+        #      The first line sort() can cause problems. If the largest other_dims is always
+        #      the last dim, tail(other_dims) is enough.
         unlimited_dim <- tail(sort(tail(other_dims_ordered_by_size, other_dims_per_chunk)), 1)
         #unlimited_dim <- tail(other_dims)
       }
diff --git a/tests/testthat/test-CDORemap.R b/tests/testthat/test-CDORemap.R
index 16da1b0..1ace086 100644
--- a/tests/testthat/test-CDORemap.R
+++ b/tests/testthat/test-CDORemap.R
@@ -39,6 +39,10 @@ data4_1 <- drop(data4)
 data4_2 <- ClimProjDiags::Subset(data4, c(1,2,3,4,8), list(1,1,1,1,1), drop = 'selected')
 data4_3 <- .aperm2(data4_2, c(1, 3, 2))
 
+ # data5: regular grid, more dimensions  
+ data5 <- array(1:(4*10*8*3*2), dim = c(dat = 1, var = 1, memb = 4, lon = 10, lat = 8, sdate = 3, sweek = 2))
+ data5_1 <- aperm(data5, c(1,2,3,6,4,5,7))
+
 ##############################################
 
 test_that("1. Input checks", {
@@ -269,3 +273,43 @@ c(-88.2, -84.6, -81.0, -77.4, -73.8)
 )
 
 })
+
+
+############################################################
+
+test_that("6. dat5: regular regrid, more dimensions", {
+suppressWarnings(
+res5 <- CDORemap(data5, lons = seq(1, 5.5, by = 0.5), lats = seq(10, 13.5, by = 0.5), 
+                 grid = 'r360x181', method = 'bil', crop = T)
+)
+expect_equal(
+as.vector(res5$lons),
+1:5
+)
+expect_equal(
+as.vector(res5$lats),
+10:13
+)
+expect_equal(
+as.vector(res5$data_array[1, 1, 1, , , 2, 1]),
+c(seq(321, 353, by = 8), seq(401, 433, by = 8), seq(481, 513, by = 8), seq(561, 593, by = 8))
+)
+
+suppressWarnings(
+res5_1 <- CDORemap(data5_1, lons = seq(1, 5.5, by = 0.5), lats = seq(10, 13.5, by = 0.5),
+                 grid = 'r360x181', method = 'bil', crop = T)
+)
+expect_equal(
+aperm(res5_1$data_array, c(1, 2, 3, 5, 6, 4, 7)),
+res5$data_array
+)
+expect_equal(
+as.vector(res5_1$lons),
+1:5
+)
+expect_equal(
+as.vector(res5_1$lats),
+10:13
+)
+})
+
-- 
GitLab