context("Transform with indices")
# Using indinces() to assign lat and lon, and transform the data. 
# Also test transform + chunk along lat/lon.

#----------------------------------------------------------
# cdo result
#library(easyNCDF)
#pathh <- "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc"
#file <- NcOpen(pathh)
#arr <- NcToArray(file,
#                 dim_indices = list(time = 1, ensemble = 1,
#                                    latitude = 1:640, longitude = 1:1296),
#                 vars_to_read = 'tas')
#lats <- NcToArray(file,
#                  dim_indices = list(latitude = 1:640), vars_to_read = 'latitude')
#lons <- NcToArray(file,
#                  dim_indices = list(longitude = 1:1296), vars_to_read = 'longitude')
#NcClose(file)
#
#arr2 <- s2dv::CDORemap(arr, lons = as.vector(lons), lats = as.vector(lats),
#                       grid = 'r100x50', method = 'con', crop = FALSE)


test_that("1. global", {
#skip_on_cran()

path <- '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc'

#-----------------------------------
# crop = region
suppressWarnings(
exp <- Start(dat = path,
              var = 'tas',
              sdate = '20000101',
              ensemble = indices(1),
              time = indices(1),
              lat = indices(1:640),
              lon = indices(1:1296),
              lat_reorder = Sort(),
              lon_reorder = CircularSort(0, 360),
              transform = CDORemapper,
              transform_extra_cells = 8,
              transform_params = list(grid = 'r100x50', method = 'conservative'),
#                                      crop = c(0, 360, -90, 90)),
              transform_vars = c('lat','lon'),
              synonims = list(lon = c('lon', 'longitude'), lat = c('lat', 'latitude')),
              return_vars = list(lat = 'dat', lon = 'dat', time = 'sdate'),
              retrieve = F)
)
func <- function(x) {
  return(x)
}
step <- Step(func, target_dims = 'time', output_dims = 'time')
wf <- AddStep(exp, step)

#---
suppressWarnings(
resT <- eval(exp)
)
suppressWarnings(
res1 <- Compute(wf, chunks = list(lon = 2))
)
#suppressWarnings(
#res2 <- Compute(wf, chunks = list(ensemble = 1))
#)
suppressWarnings(
res3 <- Compute(wf, chunks = list(lon = 3))
)


expect_equal(
as.vector(res1$output1),
as.vector(resT)
)
expect_equal(
res1$output1,
res3$output1
)

#-----------------------------------

# crop = region, selector is indices(list(, ))
suppressWarnings(
exp <- Start(dat = path,
              var = 'tas',
              sdate = '20000101',
              ensemble = indices(1),
              time = indices(1),
              lat = indices(list(1, 640)),
              lon = indices(list(1, 1296)),
              lat_reorder = Sort(),
              lon_reorder = CircularSort(0, 360),
              transform = CDORemapper,
              transform_extra_cells = 8,
              transform_params = list(grid = 'r100x50', method = 'conservative'),
#                                      crop = c(0, 360, -90, 90)),
              transform_vars = c('lat','lon'),
              synonims = list(lon = c('lon', 'longitude'), lat = c('lat', 'latitude')),
              return_vars = list(lat = 'dat', lon = 'dat', time = 'sdate'),
              retrieve = F)
)
func <- function(x) {
  return(x)
}
step <- Step(func, target_dims = 'time', output_dims = 'time')
wf <- AddStep(exp, step)

#---
suppressWarnings(
res1_list <- Compute(wf, chunks = list(lon = 2))
)
expect_equal(
res1$output1,
res1_list$output1
)

})


#####################################################################
#####################################################################
#####################################################################

#NOTE: The numbers in the unit test are testified by the following code. First, we subset
#      the desired region plus the extra cells (that is, the desired lon is (19:65), so we
#      subset (19-8:65+8)); then, we use CDORemap() to transform and crop like the Start()
#      call. Actually, the crop region is not correct. The lat region should be (-90, 67)
#      rather than (-90, -60). It causes wrong values at the end of lat because the selected
#      region is not big enough to do the interpolation at -60. But anyway, the startR
#      result is identical to arr2, and that's what we expect.

#pathh <- "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc"
#file <- NcOpen(pathh)
#arr <- NcToArray(file,
#                 dim_indices = list(time = 1, ensemble = 1,
#                                    latitude = 553:640, longitude = 11:83),
#                 vars_to_read = 'tas')
#lats <- NcToArray(file,
#                  dim_indices = list(latitude = 553:640), vars_to_read = 'latitude')
#lons <- NcToArray(file,
#                  dim_indices = list(longitude = 11:83), vars_to_read = 'longitude')
#NcClose(file)
#
#arr2 <- s2dv::CDORemap(arr, lons = as.vector(lons), lats = as.vector(lats),
#                       grid = 'r100x50', method = 'con', crop = c(0, 22, -90, -60))


test_that("2. regional, no border", {
#skip_on_cran()

path <- '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc'

# crop = region
suppressWarnings(
exp <- Start(dat = path,
              var = 'tas',
              sdate = '20000101',
              ensemble = indices(1),
              time = indices(1),
              lat = indices(1:80), # 1:80 = -89.78488:-67.58778
              lon = indices(19:65),  # 19:65 = 5.00000:17.7777778
              lat_reorder = Sort(),
              lon_reorder = CircularSort(0, 360),
              transform = CDORemapper,
              transform_extra_cells = 8,
              transform_params = list(grid = 'r100x50', 
                                      method = 'conservative'), 
              transform_vars = c('lat','lon'),
              synonims = list(lon = c('lon', 'longitude'), lat = c('lat', 'latitude')),
              return_vars = list(lat = 'dat', lon = 'dat', time = 'sdate'),
              retrieve = F)
)
func <- function(x) {
  return(x)
}
step <- Step(func, target_dims = 'time', output_dims = 'time')
wf <- AddStep(exp, step)

#---
suppressWarnings(
resT <- eval(exp)
)
suppressWarnings(
res1 <- Compute(wf, chunks = list(lon = 2))
)
#suppressWarnings(
#res2 <- Compute(wf, chunks = list(ensemble = 1))
#)
suppressWarnings(
res3 <- Compute(wf, chunks = list(lon = 3))
)

expect_equal(
as.vector(res1$output1),
as.vector(resT)
)
expect_equal(
res1$output1,
res3$output1
)

expect_equal(
drop(res1$output1)[, 1],
c(241.4042, 242.5804, 246.8507, 245.8008, 246.4318, 267.0983),
tolerance = 0.001
)
expect_equal(
drop(res1$output1)[, 2],
c(241.2223, 242.2564, 245.9863, 244.5377, 244.8937, 266.5749),
tolerance = 0.001
)
expect_equal(
drop(res1$output1)[, 3],
c(241.0895, 242.1896, 245.3183, 243.1169, 243.9446, 266.4386),
tolerance = 0.001
)


})

#####################################################################
#####################################################################
#####################################################################

#NOTE: The numbers in the unit test below is identical to the result from the following
#      code. Unlike unit test 2 above, we need to retrieve the global grids for 
#      transformation here because lon is at the border and the extra cells at the other
#      side (i.e., 360, 359, etc.) are needed.

#library(easyNCDF)
#pathh <- "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc"
#file <- NcOpen(pathh)
#arr <- NcToArray(file,
#                 dim_indices = list(time = 1, ensemble = 1,
#                                    latitude = 1:640, longitude = 1:1296),
#                 vars_to_read = 'tas')
#lats <- NcToArray(file,
#                  dim_indices = list(latitude = 1:640), vars_to_read = 'latitude')
#lons <- NcToArray(file,
#                  dim_indices = list(longitude = 1:1296), vars_to_read = 'longitude')
#NcClose(file)
#
#arr2 <- s2dv::CDORemap(arr, lons = as.vector(lons), lats = as.vector(lats),
#                       grid = 'r100x50', method = 'con', crop = c(0, 18, -90, -67))


test_that("3. regional, at lon border", {
#skip_on_cran()

path <- '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc'

# crop = region
suppressWarnings(
exp <- Start(dat = path,
              var = 'tas',
              sdate = '20000101',
              ensemble = indices(1),
              time = indices(1),
              lat = indices(1:80), # 1:80 = -89.78488:-67.58778
              lon = indices(1:65),# 1:65 = 0.00000:17.7777778
              lat_reorder = Sort(),
              lon_reorder = CircularSort(0, 360),
              transform = CDORemapper,
              transform_extra_cells = 8,
              transform_params = list(grid = 'r100x50', method = 'conservative'), 
              transform_vars = c('lat','lon'),
              synonims = list(lon = c('lon', 'longitude'), lat = c('lat', 'latitude')),
              return_vars = list(lat = 'dat', lon = 'dat', time = 'sdate'),
              retrieve = F)
)

func <- function(x) {
  return(x)
}
step <- Step(func, target_dims = 'time', output_dims = 'time')
wf <- AddStep(exp, step)

suppressWarnings(
resT <- eval(exp)
)
suppressWarnings(
res1 <- Compute(wf, chunks = list(lon = 2))
)
#suppressWarnings(
#res2 <- Compute(wf, chunks = list(ensemble = 1))
#)
suppressWarnings(
res3 <- Compute(wf, chunks = list(lon = 3))
)

expect_equal(
as.vector(res1$output1),
as.vector(resT)
)
expect_equal(
res1$output1,
res3$output1
)
expect_equal(
drop(res1$output1)[, 1],
c(241.8592, 243.7243, 248.7337, 247.9308, 252.0744, 268.5533),
tolerance = 0.001
)
expect_equal(
drop(res1$output1)[, 2],
c(241.6231, 243.0969, 247.8179, 246.8879, 249.1226, 267.8804),
tolerance = 0.001
)
expect_equal(
drop(res1$output1)[, 3],
c(241.4042, 242.5804, 246.8507, 245.8008, 246.4318, 267.0983),
tolerance = 0.001
)
expect_equal(
drop(res1$output1)[, 4],
c(241.2223, 242.2564, 245.9863, 244.5377, 244.8937, 266.5749),
tolerance = 0.001
)
expect_equal(
drop(res1$output1)[, 5],
c(241.0894, 242.1896, 245.3183, 243.1169, 243.9446, 266.4386),
tolerance = 0.001
)

#---------------------------------
# lat indices is reversed

suppressWarnings(
exp <- Start(dat = path,
              var = 'tas',
              sdate = '20000101',
              ensemble = indices(1),
              time = indices(1),
              lat = indices(80:1), # 1:80 = -89.78488:-67.58778
              lon = indices(1:65),# 1:65 = 0.00000:17.7777778
              lat_reorder = Sort(),
              lon_reorder = CircularSort(0, 360),
              transform = CDORemapper,
              transform_extra_cells = 8,
              transform_params = list(grid = 'r100x50', method = 'conservative'),
              transform_vars = c('lat','lon'),
              synonims = list(lon = c('lon', 'longitude'), lat = c('lat', 'latitude')),
              return_vars = list(lat = 'dat', lon = 'dat', time = 'sdate'),
              retrieve = F)
)

func <- function(x) {
  return(x)
}
step <- Step(func, target_dims = 'time', output_dims = 'time')
wf <- AddStep(exp, step)

#---
suppressWarnings(
resT <- eval(exp)
)
suppressWarnings(
res4 <- Compute(wf, chunks = list(lon = 2))
)
suppressWarnings(
res5 <- Compute(wf, chunks = list(lat = 2))
)
expect_equal(
as.vector(res4$output1),
as.vector(resT)
)
expect_equal(
res4$output1,
res5$output1
)
expect_equal(
as.vector(drop(res1$output1)[6:1, ]),
as.vector(drop(res4$output1))
)

})