Newer
Older
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", {
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)
resT <- eval(exp)
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
)
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
})
#####################################################################
#####################################################################
#####################################################################
#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", {
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',
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)
resT <- eval(exp)
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)
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),
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
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", {
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)
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)
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
)
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),
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)
)
expect_equal(
as.vector(res4$output1),
as.vector(resT)
)
as.vector(drop(res1$output1)[6:1, ]),
as.vector(drop(res4$output1))