diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 8cd9e963bf92843ac7de669222149d1c4dbddfc4..ef540fe63c5bdfb2b3ee3bfb3a3e9e6513f1b2a5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -7,5 +7,4 @@ build: - module load CDO/1.9.8-foss-2015a - R CMD build --resave-data . - R CMD check --as-cran --no-manual --run-donttest startR_*.tar.gz - - R -e 'covr::package_coverage()' diff --git a/R/Start.R b/R/Start.R index b6e32637d9bee67cecbd387b765c9421cea5a305..2e2258a50d0bbb4262e12412f9ecec001053e713 100644 --- a/R/Start.R +++ b/R/Start.R @@ -1330,13 +1330,21 @@ Start <- function(..., # dim = indices/selectors, "the Start call.")) } } - if (attr(dat_selectors[[dim_name]], 'indices') & (dim_name %in% transform_vars) & - !(dim_name %in% names(var_params))) { - var_params <- c(var_params, setNames(list(dim_name), dim_name)) - .warning(paste0("Found dimension '", dim_name, "' is required to transform but no '", + + if (attr(dat_selectors[[dim_name]], 'indices') & !(dim_name %in% names(var_params))) { + if (dim_name %in% transform_vars) { + var_params <- c(var_params, setNames(list(dim_name), dim_name)) + .warning(paste0("Found dimension '", dim_name, "' is required to transform but no '", + dim_name, "_var' provided. ", '"', dim_name, "_var = '", + dim_name, "'", '"', " has been automatically added to ", + "the Start call.")) + } else if (dim_name %in% names(dim_reorder_params)) { + var_params <- c(var_params, setNames(list(dim_name), dim_name)) + .warning(paste0("Found dimension '", dim_name, "' is required to reorder but no '", dim_name, "_var' provided. ", '"', dim_name, "_var = '", dim_name, "'", '"', " has been automatically added to ", "the Start call.")) + } } } @@ -1725,7 +1733,26 @@ Start <- function(..., # dim = indices/selectors, } #//////////////////////////////////////////// - # Create 'picked_common_vars' + # Change the structure of 'dat'. If the selector is not list or it is list of 2 that represents + # range, make it as list. The dimensions that go across files will later be extended to have + # lists of lists/vectors of selectors. + for (i in 1:length(dat)) { + if (dataset_has_files[i]) { + for (inner_dim in expected_inner_dims[[i]]) { + if (!is.list(dat[[i]][['selectors']][[inner_dim]]) || # not list, or + (is.list(dat[[i]][['selectors']][[inner_dim]]) && # list of 2 that represents range + length(dat[[i]][['selectors']][[inner_dim]]) == 2 && + is.null(names(dat[[i]][['selectors']][[inner_dim]])))) { + dat[[i]][['selectors']][[inner_dim]] <- list(dat[[i]][['selectors']][[inner_dim]]) + } + } + } + } + + + # Use 'common_return_vars' and 'return_vars' to generate the initial picked(_common)_vars, + # picked(_common)_vars_ordered, and picked(_common)_vars_unorder_indices. + ## Create 'picked_common_vars' if (length(common_return_vars) > 0) { picked_common_vars <- vector('list', length = length(common_return_vars)) names(picked_common_vars) <- names(common_return_vars) @@ -1735,33 +1762,19 @@ Start <- function(..., # dim = indices/selectors, picked_common_vars_ordered <- picked_common_vars picked_common_vars_unorder_indices <- picked_common_vars - # Create 'picked_vars' + ## Create 'picked_vars' picked_vars <- vector('list', length = length(dat)) names(picked_vars) <- dat_names + if (length(return_vars) > 0) { + picked_vars <- lapply(picked_vars, function(x) { + x <- vector('list', length = length(return_vars))} ) + picked_vars <- lapply(picked_vars, setNames, names(return_vars)) + } picked_vars_ordered <- picked_vars picked_vars_unorder_indices <- picked_vars for (i in 1:length(dat)) { if (dataset_has_files[i]) { - # Put all selectors in a list of a single list/vector of selectors. - # The dimensions that go across files will later be extended to have - # lists of lists/vectors of selectors. - for (inner_dim in expected_inner_dims[[i]]) { - if (!is.list(dat[[i]][['selectors']][[inner_dim]]) || # not list, or - (is.list(dat[[i]][['selectors']][[inner_dim]]) && # list of 2 that represents range - length(dat[[i]][['selectors']][[inner_dim]]) == 2 && - is.null(names(dat[[i]][['selectors']][[inner_dim]])))) { - dat[[i]][['selectors']][[inner_dim]] <- list(dat[[i]][['selectors']][[inner_dim]]) - } - } - - if (length(return_vars) > 0) { - picked_vars[[i]] <- vector('list', length = length(return_vars)) - names(picked_vars[[i]]) <- names(return_vars) - picked_vars_ordered[[i]] <- picked_vars[[i]] - picked_vars_unorder_indices[[i]] <- picked_vars[[i]] - } - indices_of_first_file <- as.list(indices_of_first_files_with_data[[i]]) array_file_dims <- sapply(dat[[i]][['selectors']][found_file_dims[[i]]], function(x) length(x[[1]])) names(array_file_dims) <- found_file_dims[[i]] @@ -1809,179 +1822,42 @@ Start <- function(..., # dim = indices/selectors, # Need to translate accoridng to synonims: names(var_dims) <- replace_with_synonmins(var_dims, synonims) if (!is.null(var_dims)) { - var_file_dims <- NULL + + ## (1) common_return_vars if (var_to_read %in% names(common_return_vars)) { var_to_check <- common_return_vars[[var_to_read]] + list_picked_var_of_read <- generate_picked_var_of_read( + var_to_read, var_to_check, array_of_files_to_load, var_dims, + array_of_var_files = array_of_var_files[j], file_var_reader, + file_object, synonims, associated_dim_name, dim_reorder_params, + aiat, current_indices, var_params, + either_picked_vars = picked_common_vars[[var_to_read]], + either_picked_vars_ordered = picked_common_vars_ordered[[var_to_read]], + either_picked_vars_unorder_indices = picked_common_vars_unorder_indices[[var_to_read]] + ) + picked_common_vars[[var_to_read]] <- list_picked_var_of_read$either_picked_vars + picked_common_vars_ordered[[var_to_read]] <- + list_picked_var_of_read$either_picked_vars_ordered + picked_common_vars_unorder_indices[[var_to_read]] <- + list_picked_var_of_read$either_picked_vars_unorder_indices + + ## (2) return_vars } else { var_to_check <- return_vars[[var_to_read]] - } - if (any(names(dim(array_of_files_to_load)) %in% var_to_check)) { - var_file_dims <- dim(array_of_files_to_load)[which(names(dim(array_of_files_to_load)) %in% - var_to_check)] - } - if (((var_to_read %in% names(common_return_vars)) && - is.null(picked_common_vars[[var_to_read]])) || - ((var_to_read %in% names(return_vars)) && - is.null(picked_vars[[i]][[var_to_read]]))) { - if (any(names(var_file_dims) %in% names(var_dims))) { - stop("Found a requested var in 'return_var' requested for a ", - "file dimension which also appears in the dimensions of ", - "the variable inside the file.\n", array_of_var_files[j]) - } - first_sample <- file_var_reader(NULL, file_object, NULL, - var_to_read, synonims) - if (any(class(first_sample) %in% names(time_special_types()))) { - array_size <- prod(c(var_file_dims, var_dims)) - new_array <- rep(time_special_types()[[class(first_sample)[1]]](NA), array_size) - dim(new_array) <- c(var_file_dims, var_dims) - } else { - new_array <- array(dim = c(var_file_dims, var_dims)) - } - attr(new_array, 'variables') <- attr(first_sample, 'variables') - if (var_to_read %in% names(common_return_vars)) { - picked_common_vars[[var_to_read]] <- new_array - pick_ordered <- FALSE - if (var_to_read %in% unlist(var_params)) { - if (associated_dim_name %in% names(dim_reorder_params) && !aiat) { - picked_common_vars_ordered[[var_to_read]] <- new_array - pick_ordered <- TRUE - } - } - if (!pick_ordered) { - picked_common_vars_ordered[[var_to_read]] <- NULL - } - } else { - picked_vars[[i]][[var_to_read]] <- new_array - pick_ordered <- FALSE - if (var_to_read %in% unlist(var_params)) { - if (associated_dim_name %in% names(dim_reorder_params) && !aiat) { - picked_vars_ordered[[i]][[var_to_read]] <- new_array - pick_ordered <- TRUE - } - } - if (!pick_ordered) { - picked_vars_ordered[[i]][[var_to_read]] <- NULL - } - } - } else { - if (var_to_read %in% names(common_return_vars)) { - array_var_dims <- dim(picked_common_vars[[var_to_read]]) - } else { - array_var_dims <- dim(picked_vars[[i]][[var_to_read]]) - } - full_array_var_dims <- array_var_dims - if (any(names(array_var_dims) %in% names(var_file_dims))) { - array_var_dims <- array_var_dims[-which(names(array_var_dims) %in% names(var_file_dims))] - } - if (names(array_var_dims) != names(var_dims)) { - stop("Error while reading the variable '", var_to_read, "' from ", - "the file. Dimensions do not match.\nExpected ", - paste(paste0("'", names(array_var_dims), "'"), - collapse = ', '), " but found ", - paste(paste0("'", names(var_dims), "'"), - collapse = ', '), ".\n", array_of_var_files[j]) - } - if (any(var_dims > array_var_dims)) { - longer_dims <- which(var_dims > array_var_dims) - if (length(longer_dims) == 1) { - longer_dims_in_full_array <- longer_dims - if (any(names(full_array_var_dims) %in% names(var_file_dims))) { - candidates <- (1:length(full_array_var_dims))[-which(names(full_array_var_dims) %in% names(var_file_dims))] - longer_dims_in_full_array <- candidates[longer_dims] - } - padding_dims <- full_array_var_dims - padding_dims[longer_dims_in_full_array] <- var_dims[longer_dims] - - array_var_dims[longer_dims] - if (var_to_read %in% names(common_return_vars)) { - var_class <- class(picked_common_vars[[var_to_read]]) - } else { - var_class <- class(picked_vars[[i]][[var_to_read]]) - } - if (any(var_class %in% names(time_special_types()))) { - padding_size <- prod(padding_dims) - padding <- rep(time_special_types()[[var_class[1]]](NA), padding_size) - dim(padding) <- padding_dims - } else { - padding <- array(dim = padding_dims) - } - if (var_to_read %in% names(common_return_vars)) { - picked_common_vars[[var_to_read]] <- .abind2( - picked_common_vars[[var_to_read]], - padding, - names(full_array_var_dims)[longer_dims_in_full_array] - ) - } else { - picked_vars[[i]][[var_to_read]] <- .abind2( - picked_vars[[i]][[var_to_read]], - padding, - names(full_array_var_dims)[longer_dims_in_full_array] - ) - } - } else { - stop("Error while reading the variable '", var_to_read, "' from ", - "the file. Found size (", paste(var_dims, collapse = ' x '), - ") is greater than expected maximum size (", - array_var_dims, ").") - } - } - } - var_store_indices <- c(as.list(current_indices[names(var_file_dims)]), lapply(var_dims, function(x) 1:x)) - var_values <- file_var_reader(NULL, file_object, NULL, var_to_read, synonims) - if (var_to_read %in% unlist(var_params)) { - if ((associated_dim_name %in% names(dim_reorder_params)) && !aiat) { - ## Is this check really needed? - if (length(dim(var_values)) > 1) { - stop("Requested a '", associated_dim_name, "_reorder' for a dimension ", - "whose coordinate variable that has more than 1 dimension. This is ", - "not supported.") - } - ordered_var_values <- dim_reorder_params[[associated_dim_name]](var_values) - attr(ordered_var_values$x, 'variables') <- attr(var_values, 'variables') - if (!all(c('x', 'ix') %in% names(ordered_var_values))) { - stop("All the dimension reorder functions must return a list with the components 'x' and 'ix'.") - } - # Save the indices to reorder back the ordered variable values. - # This will be used to define the first round indices. - unorder <- sort(ordered_var_values$ix, index.return = TRUE)$ix - if (var_to_read %in% names(common_return_vars)) { - picked_common_vars_ordered[[var_to_read]] <- do.call('[<-', - c(list(x = picked_common_vars_ordered[[var_to_read]]), - var_store_indices, - list(value = ordered_var_values$x))) - picked_common_vars_unorder_indices[[var_to_read]] <- do.call('[<-', - c(list(x = picked_common_vars_unorder_indices[[var_to_read]]), - var_store_indices, - list(value = unorder))) - } else { - picked_vars_ordered[[i]][[var_to_read]] <- do.call('[<-', - c(list(x = picked_vars_ordered[[i]][[var_to_read]]), - var_store_indices, - list(value = ordered_var_values$x))) - picked_vars_unorder_indices[[i]][[var_to_read]] <- do.call('[<-', - c(list(x = picked_vars_unorder_indices[[i]][[var_to_read]]), - var_store_indices, - list(value = unorder))) - } - } - } - if (var_to_read %in% names(common_return_vars)) { - picked_common_vars[[var_to_read]] <- do.call('[<-', - c(list(x = picked_common_vars[[var_to_read]]), - var_store_indices, - list(value = var_values))) - # Turn time zone back to UTC if this var_to_read is 'time' - if (all(class(picked_common_vars[[var_to_read]]) == names(time_special_types))) { - attr(picked_common_vars[[var_to_read]], "tzone") <- 'UTC' - } - } else { - picked_vars[[i]][[var_to_read]] <- do.call('[<-', - c(list(x = picked_vars[[i]][[var_to_read]]), - var_store_indices, - list(value = var_values))) - # Turn time zone back to UTC if this var_to_read is 'time' - if (all(class(picked_vars[[i]][[var_to_read]]) == names(time_special_types))) { - attr(picked_vars[[i]][[var_to_read]], "tzone") <- 'UTC' - } + list_picked_var_of_read <- generate_picked_var_of_read( + var_to_read, var_to_check, array_of_files_to_load, var_dims, + array_of_var_files = array_of_var_files[j], file_var_reader, + file_object, synonims, associated_dim_name, dim_reorder_params, + aiat, current_indices, var_params, + either_picked_vars = picked_vars[[i]][[var_to_read]], + either_picked_vars_ordered = picked_vars_ordered[[i]][[var_to_read]], + either_picked_vars_unorder_indices = picked_vars_unorder_indices[[i]][[var_to_read]] + ) + picked_vars[[i]][[var_to_read]] <- list_picked_var_of_read$either_picked_vars + picked_vars_ordered[[i]][[var_to_read]] <- + list_picked_var_of_read$either_picked_vars_ordered + picked_vars_unorder_indices[[i]][[var_to_read]] <- + list_picked_var_of_read$either_picked_vars_unorder_indices } if (var_to_read %in% names(first_found_file)) { first_found_file[var_to_read] <- TRUE @@ -1990,7 +1866,7 @@ Start <- function(..., # dim = indices/selectors, common_first_found_file[var_to_read] <- TRUE } } else { - stop("Could not find variable '", var_to_read, + stop("Could not find variable '", var_to_read, "' in the file ", array_of_var_files[j]) } } @@ -2003,21 +1879,7 @@ Start <- function(..., # dim = indices/selectors, } # Once we have the variable values, we can work out the indices # for the implicitly defined selectors. - # - # Trnasforms a vector of indices v expressed in a world of - # length N from 1 to N, into a world of length M, from - # 1 to M. Repeated adjacent indices are collapsed. - transform_indices <- function(v, n, m) { - #unique2 turns e.g. 1 1 2 2 2 3 3 1 1 1 into 1 2 3 1 - unique2 <- function(v) { - if (length(v) < 2) { - v - } else { - v[c(1, v[2:length(v)] - v[1:(length(v) - 1)]) != 0] - } - } - unique2(round(((v - 1) / (n - 1)) * (m - 1))) + 1 # this rounding may generate 0s. what then? - } + beta <- transform_extra_cells dims_to_crop <- vector('list') transformed_vars <- vector('list', length = length(dat)) @@ -2032,7 +1894,10 @@ Start <- function(..., # dim = indices/selectors, if (dataset_has_files[i]) { indices <- indices_of_first_files_with_data[[i]] if (!is.null(indices)) { - if (largest_dims_length == FALSE | is.numeric(largest_dims_length)) { #old code. use the 1st valid file to determine the dims + #////////////////////////////////////////////////// + # Find data_dims + ## old code. use the 1st valid file to determine the dims + if (!largest_dims_length | is.numeric(largest_dims_length)) { file_path <- do.call("[", c(list(array_of_files_to_load), as.list(indices_of_first_files_with_data[[i]]))) # The following 5 lines should go several lines below, but were moved # here for better performance. @@ -2060,7 +1925,8 @@ Start <- function(..., # dim = indices/selectors, } } - } else { # largest_dims_length == TRUE + ## largest_dims_length = TRUE + } else { data_dims <- find_largest_dims_length( selectors_total_list[[i]], array_of_files_to_load, selector_indices_save[[i]], dat[[i]], expected_inner_dims[[i]], @@ -2070,8 +1936,11 @@ Start <- function(..., # dim = indices/selectors, # Need to translate accoridng to synonims: names(data_dims) <- replace_with_synonmins(data_dims, synonims) - } # end of if (largest_dims_length == TRUE) + } # end if (largest_dims_length == TRUE) + #////////////////////////////////////////////////// + + #/////////////////////////////////////////////////////////////////// # Transform the variables if needed and keep them apart. if (!is.null(transform) && (length(transform_vars) > 0)) { if (!all(transform_vars %in% c(names(picked_vars[[i]]), names(picked_common_vars)))) { @@ -2105,7 +1974,7 @@ Start <- function(..., # dim = indices/selectors, for (var_to_read in names(transformed_data$variables)) { if (var_to_read %in% unlist(var_params)) { associated_dim_name <- names(var_params)[which(unlist(var_params) == var_to_read)] - if ((associated_dim_name %in% names(dim_reorder_params)) && aiat) { + if ((associated_dim_name %in% names(dim_reorder_params))) { ## Is this check really needed? if (length(dim(transformed_data$variables[[associated_dim_name]])) > 1) { stop("Requested a '", associated_dim_name, "_reorder' for a dimension ", @@ -2118,7 +1987,8 @@ Start <- function(..., # dim = indices/selectors, stop("All the dimension reorder functions must return a list with the components 'x' and 'ix'.") } # Save the indices to reorder back the ordered variable values. - # This will be used to define the first round indices. + # This will be used to define the first round indices (if aiat) or second round + # indices (if !aiat). unorder <- sort(ordered_var_values$ix, index.return = TRUE)$ix if (var_to_read %in% names(picked_common_vars)) { transformed_common_vars_ordered[[var_to_read]] <- ordered_var_values$x @@ -2130,16 +2000,16 @@ Start <- function(..., # dim = indices/selectors, } } } - transformed_picked_vars <- which(names(picked_vars[[i]]) %in% names(transformed_data$variables)) - if (length(transformed_picked_vars) > 0) { - transformed_picked_vars <- names(picked_vars[[i]])[transformed_picked_vars] - transformed_vars[[i]][transformed_picked_vars] <- transformed_data$variables[transformed_picked_vars] + transformed_picked_vars_names <- which(names(picked_vars[[i]]) %in% names(transformed_data$variables)) + if (length(transformed_picked_vars_names) > 0) { + transformed_picked_vars_names <- names(picked_vars[[i]])[transformed_picked_vars_names] + transformed_vars[[i]][transformed_picked_vars_names] <- transformed_data$variables[transformed_picked_vars_names] } if (is.null(transformed_common_vars)) { - transformed_picked_common_vars <- which(names(picked_common_vars) %in% names(transformed_data$variables)) - if (length(transformed_picked_common_vars) > 0) { - transformed_picked_common_vars <- names(picked_common_vars)[transformed_picked_common_vars] - transformed_common_vars <- transformed_data$variables[transformed_picked_common_vars] + transformed_picked_common_vars_names <- which(names(picked_common_vars) %in% names(transformed_data$variables)) + if (length(transformed_picked_common_vars_names) > 0) { + transformed_picked_common_vars_names <- names(picked_common_vars)[transformed_picked_common_vars_names] + transformed_common_vars <- transformed_data$variables[transformed_picked_common_vars_names] } } } @@ -2255,6 +2125,7 @@ Start <- function(..., # dim = indices/selectors, print(str(transform)) } } + # For fri if (var_with_selectors_name %in% names(picked_vars[[i]])) { var_with_selectors <- picked_vars[[i]][[var_with_selectors_name]] var_ordered <- picked_vars_ordered[[i]][[var_with_selectors_name]] @@ -2265,10 +2136,13 @@ Start <- function(..., # dim = indices/selectors, var_unorder_indices <- picked_common_vars_unorder_indices[[var_with_selectors_name]] } n <- prod(dim(var_with_selectors)) + # if no _reorder, var_unorder_indices is NULL if (is.null(var_unorder_indices)) { var_unorder_indices <- 1:n } + # For sri if (with_transform) { + ## var in 'dat' if (var_with_selectors_name %in% names(transformed_vars[[i]])) { m <- prod(dim(transformed_vars[[i]][[var_with_selectors_name]])) if (aiat) { @@ -2276,6 +2150,22 @@ Start <- function(..., # dim = indices/selectors, var_ordered <- transformed_vars_ordered[[i]][[var_with_selectors_name]] var_unorder_indices <- transformed_vars_unorder_indices[[i]][[var_with_selectors_name]] } + # For making sri ordered later + transformed_var_unordered_indices <- transformed_vars_unorder_indices[[i]][[var_with_selectors_name]] + if (is.null(transformed_var_unordered_indices)) { + transformed_var_unordered_indices <- 1:m + } + transformed_var_with_selectors <- transformed_vars[[i]][transformed_picked_vars_names][[var_with_selectors_name]][transformed_var_unordered_indices] + # Sorting the transformed variable and working out the indices again after transform. + if (!is.null(dim_reorder_params[[var_with_selectors_name]])) { + transformed_var_with_selectors_reorder <- dim_reorder_params[[var_with_selectors_name]](transformed_var_with_selectors) + transformed_var_with_selectors <- transformed_var_with_selectors_reorder$x + transformed_var_with_selectors_unorder <- transformed_var_with_selectors_reorder$ix + } else { + transformed_var_with_selectors_unorder <- 1:length(transformed_var_with_selectors) + } + + ## var in common } else if (var_with_selectors_name %in% names(transformed_common_vars)) { m <- prod(dim(transformed_common_vars[[var_with_selectors_name]])) if (aiat) { @@ -2283,6 +2173,20 @@ Start <- function(..., # dim = indices/selectors, var_ordered <- transformed_common_vars_ordered[[var_with_selectors_name]] var_unorder_indices <- transformed_common_vars_unorder_indices[[var_with_selectors_name]] } + # For making sri ordered later + transformed_var_unordered_indices <- transformed_common_vars_unorder_indices[[var_with_selectors_name]] + if (is.null(transformed_var_unordered_indices)) { + transformed_var_unordered_indices <- 1:m + } + transformed_var_with_selectors <- transformed_common_vars[[var_with_selectors_name]][transformed_var_unordered_indices] + # Sorting the transformed variable and working out the indices again after transform. + if (!is.null(dim_reorder_params[[var_with_selectors_name]])) { + transformed_var_with_selectors_reorder <- dim_reorder_params[[var_with_selectors_name]](transformed_var_with_selectors) + transformed_var_with_selectors <- transformed_var_with_selectors_reorder$x + transformed_var_with_selectors_unorder <- transformed_var_with_selectors_reorder$ix + } else { + transformed_var_with_selectors_unorder <- 1:length(transformed_var_with_selectors) + } } if (is.null(var_unorder_indices)) { var_unorder_indices <- 1:m @@ -2433,7 +2337,11 @@ Start <- function(..., # dim = indices/selectors, #sri <- NULL } else { ordered_sri[] <- replicate(prod(var_file_dims), list(1:m)) - sri[] <- replicate(prod(var_file_dims), list(1:m)) + if (inner_dim %in% names(dim_reorder_params)) { + sri[] <- replicate(prod(var_file_dims), list(transformed_var_unordered_indices[1:m])) + } else { + sri[] <- replicate(prod(var_file_dims), list(1:m)) + } ## var_file_dims instead?? #if (!aiat) { #fri[] <- replicate(prod(var_file_dims), list(1:n)) @@ -2618,16 +2526,30 @@ Start <- function(..., # dim = indices/selectors, # The selector_checker will return either a vector of indices or a list # with the first and last desired indices. - #NOTE: goes_across_prime_meridian may be TRUE only if the selector is list + #NOTE: goes_across_prime_meridian may be TRUE only if the selector is list of values goes_across_prime_meridian <- FALSE is_circular_dim <- FALSE + # If selectors are indices and _reorder = CircularSort() is used, change + # is_circular_dim to TRUE. + if (!is.null(var_ordered) & selectors_are_indices & + !is.null(dim_reorder_params[[inner_dim]])) { + if ('circular' %in% names(attributes(dim_reorder_params[[inner_dim]]))) { + is_circular_dim <- attr(dim_reorder_params[[inner_dim]], "circular") + if (is_circular_dim & is.list(sub_array_of_selectors)) { + tmp <- dim_reorder_params[[inner_dim]](unlist(sub_array_of_selectors))$ix + goes_across_prime_meridian <- tmp[1] > tmp[2] + } + } + } + + # If selectors are values and _reorder is defined. if (!is.null(var_ordered) && !selectors_are_indices) { if (!is.null(dim_reorder_params[[inner_dim]])) { + if ('circular' %in% names(attributes(dim_reorder_params[[inner_dim]]))) { + is_circular_dim <- attr(dim_reorder_params[[inner_dim]], "circular") + } if (is.list(sub_array_of_selectors)) { ## NOTE: The check of 'goes_across_prime_meridian' is moved forward to here. - if ('circular' %in% names(attributes(dim_reorder_params[[inner_dim]]))) { - is_circular_dim <- attr(dim_reorder_params[[inner_dim]], "circular") - } if (is_circular_dim) { # NOTE: Use CircularSort() to put the values in the assigned range, and get the order. # For example, [-10, 20] in CircularSort(0, 360) is [350, 20]. The $ix list is [2, 1]. @@ -2639,8 +2561,12 @@ Start <- function(..., # dim = indices/selectors, goes_across_prime_meridian <- tmp[1] > tmp[2] } - # HERE change to the same code as below (under 'else'). Not sure why originally - #it uses additional lines, which make reorder not work. + #NOTE: HERE change to the same code as below (under 'else'). Not sure why originally + # it uses additional lines, which make reorder not work. + # If "_reorder" is used, here 'sub_array_of_selectors' is adjusted to + # follow the reorder rule. E.g., if lat = values(list(-90, 90)) and + # lat_reorder = Sort(decreasing = T), 'sub_array_of_selectors' changes + # from list(-90, 90) to list(90, -90). sub_array_of_selectors <- as.list(dim_reorder_params[[inner_dim]](unlist(sub_array_of_selectors))$x) #sub_array_reordered <- dim_reorder_params[[inner_dim]](unlist(sub_array_of_selectors)) #sub_array_unorder <- sort(sub_array_reordered$ix, index.return = TRUE)$ix @@ -2696,7 +2622,7 @@ Start <- function(..., # dim = indices/selectors, } else { # Add warning if the boundary is out of range - if (is.list(sub_array_of_selectors)) { + if (is.list(sub_array_of_selectors) & !selectors_are_indices) { if (min(unlist(sub_array_of_selectors)) < min(sub_array_of_values)) { show_out_of_range_warning(inner_dim, range = range(sub_array_of_values), bound = 'lower') @@ -2707,6 +2633,9 @@ Start <- function(..., # dim = indices/selectors, } } + # sub_array_of_values here is NULL if selectors are indices, and + # 'sub_array_of_indices' will be sub_array_of_selectors, i.e., the indices + # assigned (but rounded). sub_array_of_indices <- selector_checker(sub_array_of_selectors, sub_array_of_values, tolerance = if (aiat) { NULL @@ -2722,28 +2651,66 @@ Start <- function(..., # dim = indices/selectors, } } - ## This 'if' runs in both Start() and Compute(). In Start(), it doesn't have any effect (no chunk). - ## In Compute(), it creates the indices for each chunk. For example, if 'sub_array_of_indices' - ## is c(5:10) and chunked into 2, 'sub_array_of_indices' becomes c(5:7) for chunk = 1, c(8:10) - ## for chunk = 2. If 'sub_array_of_indices' is list(55, 62) and chunked into 2, it becomes - ## list(55, 58) for chunk = 1 and list(59, 62) for chunk = 2. - ## TODO: The list can be turned into vector here? So afterward no need to judge if it is list - ## or vector. - if (!is.list(sub_array_of_indices)) { - sub_array_of_indices <- sub_array_of_indices[chunk_indices(length(sub_array_of_indices), - chunks[[inner_dim]]["chunk"], - chunks[[inner_dim]]["n_chunks"], - inner_dim)] - } else { - tmp <- chunk_indices(length(sub_array_of_indices[[1]]:sub_array_of_indices[[2]]), - chunks[[inner_dim]]["chunk"], chunks[[inner_dim]]["n_chunks"], - inner_dim) - vect <- sub_array_of_indices[[1]]:sub_array_of_indices[[2]] - sub_array_of_indices[[1]] <- vect[tmp[1]] - sub_array_of_indices[[2]] <- vect[tmp[length(tmp)]] - } + # If chunking along this inner dim, this part creates the indices for each chunk. + + # For example, if 'sub_array_of_indices' is c(5:10) and chunked into 2, + # 'sub_array_of_indices' becomes c(5:7) for chunk = 1, c(8:10) for chunk = 2. + # If 'sub_array_of_indices' is list(55, 62) and chunked into 2, it becomes + # list(55, 58) for chunk = 1 and list(59, 62) for chunk = 2. + #TODO: The list can be turned into vector here? So afterward no need to judge if + # it is list or vector. + #NOTE: chunking cannot be done if goes_across_prime_meridian = TRUE. + #TODO: Change the algorithm to make chunking works for goes_across_prime_meridian = TRUE. + # If goes_across_prime_meridian = TRUE, "sub_array_of_indices" are not + # continuous numbers. For example, list(37, 1243) means sub_array_of_fri + # that will be generated based on sub_array_of_indices later is c(1:37, 1243:1296). + # the longitude are separated into 2 parts, therefore, cannot be chunked here. + if (chunks[[inner_dim]]["n_chunks"] > 1) { + if (goes_across_prime_meridian) { + stop(paste0("Chunking over ", inner_dim, " that goes across the circular border assigned by '", inner_dim, "_reorder' is not supported by startR now. Adjust the ", inner_dim, " selector to be within the border or change the borders." )) + } + if (!is.list(sub_array_of_indices)) { + sub_array_of_indices <- + sub_array_of_indices[chunk_indices(length(sub_array_of_indices), + chunks[[inner_dim]]["chunk"], + chunks[[inner_dim]]["n_chunks"], + inner_dim)] + } else { + tmp <- + chunk_indices(length(sub_array_of_indices[[1]]:sub_array_of_indices[[2]]), + chunks[[inner_dim]]["chunk"], chunks[[inner_dim]]["n_chunks"], + inner_dim) + vect <- sub_array_of_indices[[1]]:sub_array_of_indices[[2]] + sub_array_of_indices[[1]] <- vect[tmp[1]] + sub_array_of_indices[[2]] <- vect[tmp[length(tmp)]] + } + } # The sub_array_of_indices now contains numeric indices of the values to be taken by each chunk. - + + #---------------------------------------------------------- + # 'sub_sub_array_of_values' is for sri chunking. If this inner dim is chunked, + # the sri has to follow the chunking of fri. Therefore, we save the original + # value of this chunk here for later use. We'll find the corresponding + # transformed value within 'sub_sub_array_of_values' and chunk sri. + if (with_transform & chunks[[inner_dim]]["n_chunks"] > 1) { + if (!is.null(var_ordered)) { #var_ordered + input_array_of_values <- var_ordered + } else { + if (is.null(sub_array_of_values)) { # selectors are indices + #NOTE: Not sure if 'vars_to_transform' is the correct one to use. + input_array_of_values <- vars_to_transform[[var_with_selectors_name]] + } else { + input_array_of_values <- sub_array_of_values + } + } + tmp <- generate_sub_sub_array_of_values( + input_array_of_values, sub_array_of_indices, + number_of_chunk = chunks[[inner_dim]]["chunk"]) + sub_sub_array_of_values <- tmp$sub_sub_array_of_values + previous_sub_sub_array_of_values <- tmp$previous_sub_sub_array_of_values + } + #---------------------------------------------------------- + if (debug) { if (inner_dim %in% dims_to_check) { print("-> TRANSFORMATION REQUESTED?") @@ -2768,6 +2735,12 @@ Start <- function(..., # dim = indices/selectors, with_transform, goes_across_prime_meridian, sub_array_of_indices, n, beta, is_circular_dim) + # May be useful for crop = T. 'subset_vars_to_transform' may not need + # to include extra cells, but currently it shows mistake if not include. + sub_array_of_fri_no_beta <- generate_sub_array_of_fri( + with_transform, goes_across_prime_meridian, sub_array_of_indices, n, beta, + is_circular_dim, add_beta = FALSE) + subset_vars_to_transform <- vars_to_transform if (!is.null(var_ordered)) { @@ -2780,16 +2753,23 @@ Start <- function(..., # dim = indices/selectors, subset_vars_to_transform[[var_with_selectors_name]] <- Subset(var_ordered, inner_dim, sub_array_of_fri) } else { - ##NOTE: It should be redundant because without reordering the var should remain array + if (!selectors_are_indices) { # selectors are values + #NOTE: It should be redundant because without reordering the var should remain array ## But just stay same with above... if (!is.array(sub_array_of_values)) { sub_array_of_values <- as.array(sub_array_of_values) names(dim(sub_array_of_values)) <- inner_dim } - + subset_vars_to_transform[[var_with_selectors_name]] <- Subset(sub_array_of_values, inner_dim, sub_array_of_fri) + + } else { # selectors are indices + subset_vars_to_transform[[var_with_selectors_name]] <- + Subset(subset_vars_to_transform[[var_with_selectors_name]], + inner_dim, sub_array_of_fri) + } } - + # Change the order of longitude crop if no reorder + from big to small. # cdo -sellonlatbox, the lon is west, east (while lat can be north # to south or opposite) @@ -2827,12 +2807,36 @@ Start <- function(..., # dim = indices/selectors, } else { transformed_subset_var_unorder <- 1:length(transformed_subset_var) } - sub_array_of_sri <- selector_checker(sub_array_of_selectors, transformed_subset_var, - tolerance = if (aiat) { - tolerance_params[[inner_dim]] - } else { - NULL - }) + if (!selectors_are_indices) { # selectors are values + sub_array_of_sri <- selector_checker( + sub_array_of_selectors, transformed_subset_var, + tolerance = if (aiat) { + tolerance_params[[inner_dim]] + } else { + NULL + }) + if (!is.list(sub_array_of_sri)) { + sub_array_of_sri <- unique(sub_array_of_sri) + } + } else { # selectors are indices + # Need to transfer to values first, then use the values to get the new + # indices in transformed_subset_var. + if (is.list(sub_array_of_selectors)) { + ori_values <- vars_to_transform[[var_with_selectors_name]][sub_array_of_selectors[[1]]:sub_array_of_selectors[[2]]] + } else { + ori_values <- vars_to_transform[[var_with_selectors_name]][sub_array_of_selectors] + } + sub_array_of_sri <- selector_checker( + ori_values, transformed_subset_var, + tolerance = if (aiat) { + tolerance_params[[inner_dim]] + } else { + NULL + }) + # Here may need to further modify considering aiat. If aiat = FALSE, + # (i.e., indices are taken before transform), unique() is needed. + sub_array_of_sri <- unique(sub_array_of_sri) + } # Check if selectors fall out of the range of the transform grid # It may happen when original lon is [-180, 180] while want to regrid to @@ -2852,17 +2856,89 @@ Start <- function(..., # dim = indices/selectors, sub_array_of_sri <- c(1:length(transformed_subset_var)) } else { # the common case, i.e., non-global - # NOTE: Because sub_array_of_sri order is exchanged due to - # previous development, here [[1]] and [[2]] should exchange - sub_array_of_sri <- c(1:sub_array_of_sri[[1]], - sub_array_of_sri[[2]]:length(transformed_subset_var)) +# # NOTE: Because sub_array_of_sri order is exchanged due to +# # previous development, here [[1]] and [[2]] should exchange +# sub_array_of_sri <- c(1:sub_array_of_sri[[1]], +# sub_array_of_sri[[2]]:length(transformed_subset_var)) + #NOTE: the old code above is not suitable for all the possible cases. + # If sub_array_of_selectors is not exactly the value in transformed_subset_var, sub_array_of_sri[[1]] will be larger than sub_array_of_sri[[2]]. + # Though here is not global case, we already have transformed_subset_var cropped as the desired region, so it is okay to use the whole length. Not sure if it will cause other problems... + sub_array_of_sri <- 1:length(transformed_subset_var) } } else if (is.list(sub_array_of_sri)) { sub_array_of_sri <- sub_array_of_sri[[1]]:sub_array_of_sri[[2]] } + + # Chunk sub_array_of_sri if this inner_dim needs to be chunked + #TODO: Potential problem: the transformed_subset_var value falls between + # the end of sub_sub_array_of_values of the 1st chunk and the beginning + # of sub_sub_array_of_values of the 2nd chunk. Then, one sub_array_of_sri + # will miss. 'previous_sri' is checked and will be included if this + # situation happens, but don't know if the transformed result is + # correct or not. + if (chunks[[inner_dim]]["n_chunks"] > 1) { + if (is.list(sub_sub_array_of_values)) { # list + sub_array_of_sri <- + which(transformed_subset_var >= min(unlist(sub_sub_array_of_values)) & + transformed_subset_var <= max(unlist(sub_sub_array_of_values))) + # Check if sub_array_of_sri perfectly connects to the previous sri. + # If not, inlclude the previous sri. + #NOTE 1: don't know if the transform for the previous sri is + # correct or not. + #NOTE 2: If crop = T, sub_array_of_sri always starts from 1. + # Don't know if the cropping will miss some sri or not. + if (sub_array_of_sri[1] != 1) { + if (!is.null(previous_sub_sub_array_of_values)) { + previous_sri <- max(which(transformed_subset_var <= previous_sub_sub_array_of_values)) + if (previous_sri + 1 != sub_array_of_sri[1]) { + sub_array_of_sri <- (previous_sri + 1):sub_array_of_sri[length(sub_array_of_sri)] + } + } + } + + } else { # is vector + tmp <- which(transformed_subset_var >= min(sub_sub_array_of_values) & + transformed_subset_var <= max(sub_sub_array_of_values)) + # Include first or last sri if tmp doesn't have. It's only for + # ""vectors"" because vectors look for the closest value. + #NOTE: The condition here is not correct. The criteria should be + # 'vector' instead of indices. + if (chunks[[inner_dim]]["chunk"] == 1) { + sub_array_of_sri <- unique(c(sub_array_of_sri[1], tmp)) + } else if (chunks[[inner_dim]]["chunk"] == + chunks[[inner_dim]]["n_chunks"]) { # last chunk + sub_array_of_sri <- unique(c(tmp, sub_array_of_sri[length(sub_array_of_sri)])) + } else { + sub_array_of_sri <- tmp + } + # Check if sub_array_of_sri perfectly connects to the previous sri. + # If not, inlclude the previous sri. + #NOTE 1: don't know if the transform for the previous sri is + # correct or not. + #NOTE 2: If crop = T, sub_array_of_sri always starts from 1. + # Don't know if the cropping will miss some sri or not. + if (sub_array_of_sri[1] != 1) { + if (!is.null(previous_sub_sub_array_of_values)) { + previous_sri <- max(which(transformed_subset_var <= previous_sub_sub_array_of_values)) + if (previous_sri + 1 != sub_array_of_sri[1]) { + sub_array_of_sri <- (previous_sri + 1):sub_array_of_sri[length(sub_array_of_sri)] + } + } + } + } + } ordered_sri <- sub_array_of_sri sub_array_of_sri <- transformed_subset_var_unorder[sub_array_of_sri] + +###########################old################################## +# if (chunks[[inner_dim]]["n_chunks"] > 1) { +# tmp <- which(transformed_subset_var >= min(sub_sub_array_of_values) & +# transformed_subset_var <= max(sub_sub_array_of_values)) +# sub_array_of_sri <- sub_array_of_sri[tmp] +# } +################################################################ + # In this case, the tvi are not defined and the 'transformed_subset_var' # will be taken instead of the var transformed before in the code. if (debug) { @@ -2918,6 +2994,9 @@ Start <- function(..., # dim = indices/selectors, with_transform, goes_across_prime_meridian, sub_array_of_indices, n, beta, is_circular_dim) } + + # Reorder sub_array_of_fri if reordering function is used. + # It was index in the assigned order (e.g., in [-180, 180] if CircularSort(-180, 180)), and here is changed to the index in the original order. if (!is.null(var_unorder_indices)) { if (is.null(ordered_fri)) { ordered_fri <- sub_array_of_fri @@ -3153,6 +3232,9 @@ Start <- function(..., # dim = indices/selectors, if (!(length(selector_array) == 1 & selector_array %in% c('all', 'first', 'last'))) { vars_to_crop[[var_to_crop]] <- Subset(transformed_subset_var, inner_dim, crop_indices) + } else { + vars_to_crop[[var_to_crop]] <- + Subset(transformed_var_with_selectors, inner_dim, crop_indices) } } else { vars_to_crop[[var_to_crop]] <- Subset(vars_to_crop[[var_to_crop]], inner_dim, crop_indices) @@ -3172,6 +3254,9 @@ Start <- function(..., # dim = indices/selectors, selector_array %in% c('all', 'first', 'last'))) { common_vars_to_crop[[common_var_to_crop]] <- Subset(transformed_subset_var, inner_dim, crop_indices) + } else { + common_vars_to_crop[[common_var_to_crop]] <- + Subset(transformed_var_with_selectors, inner_dim, crop_indices) } } else { #old code common_vars_to_crop[[common_var_to_crop]] <- Subset(common_vars_to_crop[[common_var_to_crop]], inner_dim, crop_indices) diff --git a/R/zzz.R b/R/zzz.R index b045e85e0d275a1fbbe0e1b932944626c5a009be..3eeed1a5c667e6deb88d861f5e100fd6593fa4e4 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -568,10 +568,37 @@ show_out_of_range_warning <- function(inner_dim, range, bound) { "Check if the desired range is all included.")) } +# 'sub_sub_array_of_values' is for sri chunking. If this inner dim is chunked, +# the sri has to follow the chunking of fri. Therefore, we save the original +# value of this chunk here for later use. We'll find the corresponding +# transformed value within 'sub_sub_array_of_values' and chunk sri. This +# function also returns 'previous_sub_subarray_of_values', which is used for +# checking if there is sri being skipped. +generate_sub_sub_array_of_values <- function(input_array_of_values, sub_array_of_indices, + number_of_chunk) { + previous_sub_sub_array_of_values <- NULL + + if (is.list(sub_array_of_indices)) { + sub_sub_array_of_values <- list(input_array_of_values[sub_array_of_indices[[1]]], + input_array_of_values[sub_array_of_indices[[2]]]) + if (number_of_chunk > 1) { + previous_sub_sub_array_of_values <- input_array_of_values[sub_array_of_indices[[1]] - 1] + } + } else { # is vector + sub_sub_array_of_values <- input_array_of_values[sub_array_of_indices] + if (number_of_chunk > 1) { + previous_sub_sub_array_of_values <- input_array_of_values[sub_array_of_indices[1] - 1] + } + } + + return(list(sub_sub_array_of_values = sub_sub_array_of_values, + previous_sub_sub_array_of_values = previous_sub_sub_array_of_values)) +} + # Generate sub_array_of_fri generate_sub_array_of_fri <- function(with_transform, goes_across_prime_meridian, sub_array_of_indices, n, beta, - is_circular_dim) { + is_circular_dim, add_beta = TRUE) { print_warning <- FALSE if (goes_across_prime_meridian) { #NOTE: The potential problem here is, if it is global longitude, @@ -584,7 +611,7 @@ generate_sub_array_of_fri <- function(with_transform, goes_across_prime_meridian # global longitude sub_array_of_fri <- 1:n # n = prod(dim(var_with_selectors)) - if (with_transform & beta != 0) { + if (with_transform & beta != 0 & add_beta) { # Warning if transform_extra_cell != 0 print_warning <- TRUE } @@ -593,7 +620,7 @@ generate_sub_array_of_fri <- function(with_transform, goes_across_prime_meridian # normal case, i.e., not global first_index <- min(unlist(sub_array_of_indices)) last_index <- max(unlist(sub_array_of_indices)) - if (with_transform) { + if (with_transform & add_beta) { gap_width <- last_index - first_index - 1 actual_beta <- min(gap_width, beta) sub_array_of_fri <- c(1:(first_index + actual_beta), @@ -612,7 +639,7 @@ generate_sub_array_of_fri <- function(with_transform, goes_across_prime_meridian # sub_array_of_indices <- sub_array_of_indices[[1]]:sub_array_of_indices[[2]] # } #NOTE: sub_array_of_indices may be vector or list - if (with_transform) { + if (with_transform & add_beta) { first_index <- min(unlist(sub_array_of_indices)) last_index <- max(unlist(sub_array_of_indices)) start_padding <- min(beta, first_index - 1) @@ -624,15 +651,18 @@ generate_sub_array_of_fri <- function(with_transform, goes_across_prime_meridian print_warning <- TRUE } } else { #longitude - # if (((last_index + beta) - (first_index - beta) + 1) >= n) { - if (start_padding <= beta & end_padding <= beta) { + if (start_padding == beta & end_padding == beta) { + # normal regional situation + sub_array_of_fri <- (first_index - start_padding):(last_index + end_padding) + } else if (start_padding < beta & end_padding < beta) { + # global sub_array_of_fri <- 1:n - } else if (start_padding < beta) { # left side too close to border, need to go to right side + } else if (start_padding < beta) { + # left side too close to border, need to go to right side sub_array_of_fri <- c((first_index - start_padding):(last_index + end_padding), (n - (beta - start_padding - 1)):n) - } else if (end_padding < beta) { # right side too close to border, need to go to left side + } else if (end_padding < beta) { + # right side too close to border, need to go to left side sub_array_of_fri <- c(1: (beta - end_padding), (first_index - start_padding):(last_index + end_padding)) - } else { #normal - sub_array_of_fri <- (first_index - start_padding):(last_index + end_padding) } } @@ -1172,3 +1202,158 @@ combine_metadata_picked_vars <- function(return_metadata, picked_vars, picked_co return(Variables_list) } + +# This function generates a list of 3, containing picked(_common)_vars, +# picked(_common)_vars_ordered, and picked(_common)_vars_unorder_indices for the 'var_to_read' +# of this dataset (i) and file (j). +generate_picked_var_of_read <- function(var_to_read, var_to_check, array_of_files_to_load, + var_dims, array_of_var_files, file_var_reader, + file_object, synonims, associated_dim_name, + dim_reorder_params, aiat, current_indices, var_params, + either_picked_vars, + either_picked_vars_ordered, + either_picked_vars_unorder_indices) { + var_file_dims <- NULL + + if (any(names(dim(array_of_files_to_load)) %in% var_to_check)) { + var_file_dims <- dim(array_of_files_to_load)[which(names(dim(array_of_files_to_load)) %in% + var_to_check)] + } + if (is.null(either_picked_vars)) { + + if (any(names(var_file_dims) %in% names(var_dims))) { + stop("Found a requested var in 'return_var' requested for a ", + "file dimension which also appears in the dimensions of ", + "the variable inside the file.\n", array_of_var_files) + } + first_sample <- file_var_reader(NULL, file_object, NULL, + var_to_read, synonims) + if (any(class(first_sample) %in% names(time_special_types()))) { + array_size <- prod(c(var_file_dims, var_dims)) + new_array <- rep(time_special_types()[[class(first_sample)[1]]](NA), array_size) + dim(new_array) <- c(var_file_dims, var_dims) + } else { + new_array <- array(dim = c(var_file_dims, var_dims)) + } + attr(new_array, 'variables') <- attr(first_sample, 'variables') + + either_picked_vars <- new_array + pick_ordered <- FALSE + if (var_to_read %in% unlist(var_params)) { + if (associated_dim_name %in% names(dim_reorder_params) && !aiat) { + either_picked_vars_ordered <- new_array + pick_ordered <- TRUE + } + } + if (!pick_ordered) { + either_picked_vars_ordered <- NULL + } + + } else { + array_var_dims <- dim(either_picked_vars) + full_array_var_dims <- array_var_dims + if (any(names(array_var_dims) %in% names(var_file_dims))) { + array_var_dims <- array_var_dims[-which(names(array_var_dims) %in% names(var_file_dims))] + } + if (names(array_var_dims) != names(var_dims)) { + stop("Error while reading the variable '", var_to_read, "' from ", + "the file. Dimensions do not match.\nExpected ", + paste(paste0("'", names(array_var_dims), "'"), collapse = ', '), + " but found ", + paste(paste0("'", names(var_dims), "'"), collapse = ', '), + ".\n", array_of_var_files) + } + if (any(var_dims > array_var_dims)) { + longer_dims <- which(var_dims > array_var_dims) + if (length(longer_dims) == 1) { + longer_dims_in_full_array <- longer_dims + if (any(names(full_array_var_dims) %in% names(var_file_dims))) { + candidates <- (1:length(full_array_var_dims))[-which(names(full_array_var_dims) %in% names(var_file_dims))] + longer_dims_in_full_array <- candidates[longer_dims] + } + padding_dims <- full_array_var_dims + padding_dims[longer_dims_in_full_array] <- + var_dims[longer_dims] - array_var_dims[longer_dims] + + var_class <- class(either_picked_vars) + if (any(var_class %in% names(time_special_types()))) { + padding_size <- prod(padding_dims) + padding <- rep(time_special_types()[[var_class[1]]](NA), padding_size) + dim(padding) <- padding_dims + } else { + padding <- array(dim = padding_dims) + } + either_picked_vars <- .abind2(either_picked_vars, padding, + names(full_array_var_dims)[longer_dims_in_full_array]) + } else { + stop("Error while reading the variable '", var_to_read, "' from ", + "the file. Found size (", paste(var_dims, collapse = ' x '), + ") is greater than expected maximum size (", array_var_dims, ").") + } + } + } + + var_store_indices <- c(as.list(current_indices[names(var_file_dims)]), + lapply(var_dims, function(x) 1:x)) + var_values <- file_var_reader(NULL, file_object, NULL, var_to_read, synonims) + if (var_to_read %in% unlist(var_params)) { + if ((associated_dim_name %in% names(dim_reorder_params)) && !aiat) { + ## Is this check really needed? + if (length(dim(var_values)) > 1) { + stop("Requested a '", associated_dim_name, "_reorder' for a dimension ", + "whose coordinate variable that has more than 1 dimension. This is ", + "not supported.") + } + ordered_var_values <- dim_reorder_params[[associated_dim_name]](var_values) + attr(ordered_var_values$x, 'variables') <- attr(var_values, 'variables') + if (!all(c('x', 'ix') %in% names(ordered_var_values))) { + stop("All the dimension reorder functions must return a list with the components 'x' and 'ix'.") + } + # Save the indices to reorder the ordered variable values back to original order. + # 'unorder' refers to the indices of 'ordered_var_values' if it is unordered. + # This will be used to define the first round indices. + unorder <- sort(ordered_var_values$ix, index.return = TRUE)$ix + either_picked_vars_ordered <- do.call('[<-', + c(list(x = either_picked_vars_ordered), + var_store_indices, + list(value = ordered_var_values$x))) + either_picked_vars_unorder_indices <- do.call('[<-', + c(list(x = either_picked_vars_unorder_indices), + var_store_indices, + list(value = unorder))) + + + } + } + + either_picked_vars <- do.call('[<-', + c(list(x = either_picked_vars), + var_store_indices, + list(value = var_values))) + # Turn time zone back to UTC if this var_to_read is 'time' + if (all(class(either_picked_vars) == names(time_special_types))) { + attr(either_picked_vars, "tzone") <- 'UTC' + } + + + return(list(either_picked_vars = either_picked_vars, + either_picked_vars_ordered = either_picked_vars_ordered, + either_picked_vars_unorder_indices = either_picked_vars_unorder_indices)) +} + + +# Trnasforms a vector of indices v expressed in a world of +# length N from 1 to N, into a world of length M, from +# 1 to M. Repeated adjacent indices are collapsed. +transform_indices <- function(v, n, m) { + #unique2 turns e.g. 1 1 2 2 2 3 3 1 1 1 into 1 2 3 1 + unique2 <- function(v) { + if (length(v) < 2) { + v + } else { + v[c(1, v[2:length(v)] - v[1:(length(v) - 1)]) != 0] + } + } + unique2(round(((v - 1) / (n - 1)) * (m - 1))) + 1 # this rounding may generate 0s. what then? +} + diff --git a/inst/doc/data_check.md b/inst/doc/data_check.md index ae73f50ea2e665f11ebbaf4152f4024c2475e859..156a2531ce091bac515310c240e297c282b7619f 100644 --- a/inst/doc/data_check.md +++ b/inst/doc/data_check.md @@ -21,6 +21,7 @@ Here we list some tips recommended to pay attention, and some tools for data com - Extra examination (5) [Compute()](inst/doc/data_check.md#5-compute) + (6) [Regridding](inst/doc/data_check.md#6-regridding) ## Tips @@ -383,5 +384,35 @@ res[1:3, 1, 1, 1, 1, 1:2] ``` +### (5) Regridding +If `transform = CDORemapper` is used in Start(), you can use other regridding tool to +verify the result, like cdo or s2dv::CDORemap. Here is an example using easyNCDF and +CDORemap() to get the transformed data of file "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc". +```r +library(easyNCDF) +file <- NcOpen("/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc") +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) + +res <- s2dv::CDORemap(arr, lons = as.vector(lons), lats = as.vector(lats), + grid = 'r100x50', method = 'con', crop = FALSE) +dim(res$data_array) +# var time ensemble latitude longitude +# 1 1 1 50 100 + +head(res$lons) +#[1] 0.0 3.6 7.2 10.8 14.4 18.0 + +head(res$lats) +#[1] -88.2 -84.6 -81.0 -77.4 -73.8 -70.2 + +``` diff --git a/inst/doc/faq.md b/inst/doc/faq.md index 4cd91243666772ed3a42b3a6dc506fbe9628862d..fbf2efbc098210a04abd1d5a35cd207638e8521c 100644 --- a/inst/doc/faq.md +++ b/inst/doc/faq.md @@ -26,6 +26,8 @@ This document intends to be the first reference for any doubts that you may have 20. [Use 'metadata_dims' to retrieve variable metadata](#20-use-metadata_dims-to-retrieve-variable-metadata) 21. [Retrieve the complete data when the dimension length varies among files](#21-retrieve-the-complete-data-when-the-dimension-length-varies-among-files) 22. [Define the selector when the indices in the files are not aligned](#22-define-the-selector-when-the-indices-in-the-files-are-not-aligned) + 23. [The best practice of using vector and list for selectors](#23-the-best-practice-of-using-vector-and-list-for-selectors) + 24. [Do both interpolation and chunking on spatial dimensions](#24-do-both-interpolation-and-chunking-on-spatial-dimensions) 2. **Something goes wrong...** @@ -322,7 +324,7 @@ If you want to do the interpolation within Start(), you can use the following fo 4. **`transform_extra_cells`**: A numeric indicating the number of grid cell to extend from the borders if the interpolating region is a subset of the whole region. 2 as default, which is consistent with the method in s2dverification::Load(). You can find an example script here [ex1_1_tranform.R](/inst/doc/usecase/ex1_1_tranform.R) -You can see more information in s2dverification::CDORemap documentation [here](https://earth.bsc.es/gitlab/es/s2dverification/blob/master/man/CDORemap.Rd). +You can see more information in s2dverification::CDORemap documentation [here](https://earth.bsc.es/gitlab/es/s2dverification/blob/master/man/CDORemap.Rd). ### 6. Get data attributes without retrieving data to workstation @@ -920,6 +922,46 @@ data <- Start(dat = path, retrieve = T) ``` +### 23. The best practice of using vector and list for selectors +There are three ways to define the selectors in Start(): `indices()`, `values()`, and character string +like 'all', 'first', and 'last'. For `indices()` and `values()`, we can put either a vector or a list +in them (here we talk about the common cases, not including the dependency case mentioned in how-to-22 above.) + +For file dimensions, it is common to simply define the selectors by a vector of character string +(which belongs to `values()` but `values()` can be ommitted), e.g., `sdate = c('200001', '200002')`; `var = 'tas'`. +You can also use a vector of indices, but you cannot gurantee the files you get is the desired one +since the file order in the repository may change. + +For inner dimensions, it is recommended using "list of 2 values" or "vector of indices". +The main difference between vector and list is that the vector looks for the exact or closest +(could be larger or smaller) value in the data while the list looks for the data falling between the two numbers in the list. +You can assign all the indices needed by a vector, e.g., `time = indices(1:12)`, or give a range +that covers all the data needed by a list of 2, e.g., `lon = values(list(0, 30))`. +Note that `lon = values(list(0, 30))` means the data between 0 degE and 30 degE is taken; on the +other hand, `lon = indices(list(0, 30))` means that index 0 to index 30 of lon is taken (and it +will return an error in this case because there is no index 0.) + +In conclusion, if you know the exact values or indices of the selector, you can use vector of values or indices; if not, usually for longitude and latitude, it is better to use list of 2 values instead. + + +### 24. Do both interpolation and chunking on spatial dimensions +If all other dimensions are used as target dimensions in the operation, it would be necessary to +to chunk the spatial dimensions. The chunking can be done even if regridding is also required in +Start() (See those transform arguments at [how-to-5](#5-do-interpolation-in-start-using-parameter-transform), and the script has no difference with chunking other dimensions. +However, there are some things you need to bear in mind when using startR in this way. + +The regridding function provided by startR is CDORemapper(), which is a wrapper function of s2dv::CDORemap; +and CDORemap() uses cdo inside. Therefore, the regridding of startR has the same performance as cdo. +The errors due to transformation at borders may increase by chunking because it produces more +borders. For example, if `longitude = indices(1:20)` is chunked by 2, the first chunk will be indices(1:10) and the second chunk will be indices(11:20). Therefore, we have borders at 0, 10, 11, and 20. +In most cases, the border errors can be eliminated by increasing the number of extra cells (argument `transform_extra_cells` in Start()). With enough extra cells, the result will be identical as +global regridding. + +However, there are many factors that may impact the results of regridding, like the `crop` option, +the way to define the longitude/latitude selectors, etc. It is important to know how CDO works and +the usage of those parameters to avoid unecessary errors. +We provide some [use cases](inst/doc/usecase/ex2_12_transform_and_chunk.R) showing the secure ways of transformation + chunking. + # Something goes wrong... diff --git a/inst/doc/usecase/ex2_12_transform_and_chunk.R b/inst/doc/usecase/ex2_12_transform_and_chunk.R new file mode 100644 index 0000000000000000000000000000000000000000..8b2eb831878c0844b5262a3c501c7f8ccdbb5882 --- /dev/null +++ b/inst/doc/usecase/ex2_12_transform_and_chunk.R @@ -0,0 +1,163 @@ +# Author: An-Chi Ho +# Date: 10th September 2021 +# ------------------------------------------------------------------ +# This use case provides an example of transforming and chunking latitude and longitude +# dimensions. If all other dimensions are used as target dimensions in the operation, +# it would be good to have the option of chunking the spatial dimensions. However, the +# errors due to transformation at borders may increase because chunking produces more +# borders. There are many factors may impact the results of transformation or +# transformation + chunking. See FAQ How-to-24 for related information(https://earth.bsc.es/gitlab/es/startR/-/blob/master/inst/doc/faq.md#24-do-both-interpolation-and-chunking-on-spatial-dimensions) +# Here we provide some scripts that are more common and less error-prone. +# Common things to notice: +# 1. 'transform_extra_cells' is increased to 8 to avoid errors at borders. +# Fewer cells may be enough, depending on cases. +# 2. The 'crop' argument in 'transform_params' is defined by the borders of the region or FALSE. +# TRUE may return wrong values, depending on cases. +# 3. CircularSort() is required to use even if the longitude fully falls in the range because +# it tells startR that longitude dimension is circular and the extra cells should be got from +# the other side if the border is reached. In the scripts below, CircularSort(0, 360) is used, +# but it can also be replaced by CircularSort(-180, 180). +# ------------------------------------------------------------------ + +library(startR) + +lons.min <- 0 +lons.max <- 359.9 +lats.min <- -90 +lats.max <- 90 +sdates <- paste0(1981:2011, '0101') + +#--------------------------------- +# Method 1: +# - Use list of 2 values to define longitude and latitude. +#--------------------------------- + +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = sdates, + ensemble = 'all', + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('latitude', 'lat'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = F) + +func <- function(x) { + # x: [sdate, ensemble] + ensemble_mean <- s2dv::MeanDims(x, 2) + trend <- s2dv:::.Trend(ensemble_mean)$trend + return(trend) +} +step <- Step(func, + target_dims = c('sdate', 'ensemble'), output_dims = 'stats') +wf <- AddStep(exp, step) + +res1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 + +dim(res1) +# stats dat var time latitude longitude +# 2 1 1 1 50 100 +summary(res1) +# Min. 1st Qu. Median Mean 3rd Qu. Max. +# -0.09884 0.00879 117.19019 138.18788 279.98564 305.25010 + + + +#--------------------------------- +# Method 2: +# - Use vector of indices to define longitude and latitude. +# - The 'crop' argument in 'transform_params' is FALSE (but it can also be defined by the borders +# of the region.) +#--------------------------------- + +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = sdates, + ensemble = 'all', + time = indices(1), + latitude = indices(1:640), + latitude_reorder = Sort(), + longitude = indices(1:1296), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('latitude', 'lat'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = F) + +func <- function(x) { + # x: [sdate, ensemble] + ensemble_mean <- s2dv::MeanDims(x, 2) + trend <- s2dv:::.Trend(ensemble_mean)$trend + return(trend) +} +step <- Step(func, + target_dims = c('sdate', 'ensemble'), output_dims = 'stats') +wf <- AddStep(exp, step) + +res2 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 + +identical(res1, res2) +#[1] TRUE + +#--------------------------------- +# Method 3: +# - Use 'all' to define longitude and latitude. +# - The 'crop' argument in 'transform_params' is FALSE (but it can also be defined by the borders +# of the region.) +#--------------------------------- +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = sdates, + ensemble = 'all', + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('latitude', 'lat'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = F) + +func <- function(x) { + # x: [sdate, ensemble] + ensemble_mean <- s2dv::MeanDims(x, 2) + trend <- s2dv:::.Trend(ensemble_mean)$trend + return(trend) +} +step <- Step(func, + target_dims = c('sdate', 'ensemble'), output_dims = 'stats') +wf <- AddStep(exp, step) + +res3 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 + +identical(res2, res3) +#[1] TRUE diff --git a/tests/testthat.R b/tests/testthat.R index 5073b5e22eab429537b8a493043b4162806ccc3f..04a698ce3550b762582b7442237b9f3721350fb7 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -1,8 +1,9 @@ library(testthat) library(startR) library(SpecsVerification) -library(dplyr) library(plyr) +library(dplyr) +library(easyNCDF) library(s2dv) test_check("startR") diff --git a/tests/testthat/test-Compute-transform_all.R b/tests/testthat/test-Compute-transform_all.R index 67011ee9822b7be3b38679508966dd7e4bb5d99f..46430d7a2736f1946d6e369d32c4037179494aa8 100644 --- a/tests/testthat/test-Compute-transform_all.R +++ b/tests/testthat/test-Compute-transform_all.R @@ -1,73 +1,7 @@ context("Transform with 'all'") -test_that("1. Specify lat and lon with 'all'; retrieve = TRUE", { - -path <- '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/$member$/Amon/$var$/gr/v20190713/$var$_Amon_*_s$sdate$-$member$_gr_$fyear$.nc' -suppressWarnings( -data1 <- Start(dat = path, - var = 'tos', - sdate = paste0(1960), - time = indices(1:2), #'all', - lat = 'all', - lon = 'all', -# lat_reorder = Sort(decreasing = F), -# lon_reorder = CircularSort(0, 360), -# lat_var = 'lat', -# lon_var = 'lon', - fyear = indices(1), - member = indices(1), - transform = CDORemapper, - transform_extra_cells = 2, - transform_params = list(grid = 'r100x50', method = 'conservative', crop = FALSE), - 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 = T) -) - -suppressWarnings( -data2 <- Start(dat = path, - var = 'tos', - sdate = paste0(1960), - time = indices(1:2), #'all', - lat = values(list(-90, 90)), - lon = values(list(0, 359.9)), - lat_reorder = Sort(decreasing = F), - lon_reorder = CircularSort(0, 360), -# lat_var = 'lat', -# lon_var = 'lon', - fyear = indices(1), - member = indices(1), - transform = CDORemapper, - transform_extra_cells = 2, - transform_params = list(grid = 'r100x50', method = 'conservative', crop = FALSE), - transform_vars = c('lat','lon'), - synonims = list(lon = c('lon', 'longitude'), lat = c('lat', 'latitude')), - return_vars = list(lat = NULL, lon = NULL, time = 'sdate'), - retrieve = T) -) - -expect_equal( -dim(data1), -c(dat = 1, var = 1, sdate = 1, time = 2, lat = 50, lon = 100, fyear = 1, member = 1) -) -expect_equal( -dim(data1), -dim(data2) -) -expect_equal( -as.vector(data1), -as.vector(data2) -) -expect_equal( -data1[1, 1, 1, 2, 10:12, 20, 1, 1], -c(274.6942, 276.2658, 278.2566), -tolerance = 0.0001 -) - -}) - -test_that("2. Specify lat and lon with 'all'; retrieve = FALSE", { +test_that("1. Chunk along non-lat/lon dim", { +skip_on_cran() path <- '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/$member$/Amon/$var$/gr/v20190713/$var$_Amon_*_s$sdate$-$member$_gr_$fyear$.nc' suppressWarnings( @@ -110,5 +44,83 @@ c(274.2808, 275.8509, 277.7623), tolerance = 0.0001 ) +}) + +test_that("2. chunk along lon", { +skip_on_cran() + +#!!!!!!!!!!!!!!!!!!!NOTE: the results are not identical when exp has extra cells = 2!!!!!!!!!!!!!!!!!! +# But exp2 (retrieve = T) has the same results with extra_cells = 2 and 8. + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat',#NULL, + longitude = 'dat',#NULL, + time = 'sdate'), + retrieve = F) +) + + func <- function(x) { + return(x) + } + step <- Step(func, target_dims = 'time', output_dims = 'time') + wf <- AddStep(exp, step) +suppressWarnings( + res <- Compute(wf, chunks = list(longitude = 2)) +) +suppressWarnings( + res2 <- Compute(wf, chunks = list(ensemble = 1)) +) + +expect_equal( +res$output1, +res2$output1 +) + +# Check with retrieve = TRUE +suppressWarnings( +exp2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat',#NULL, + longitude = 'dat',#NULL, + time = 'sdate'), + retrieve = T) +) + +expect_equal( +as.vector(res$output1), +as.vector(exp2) +) }) diff --git a/tests/testthat/test-Compute-transform_indices.R b/tests/testthat/test-Compute-transform_indices.R new file mode 100644 index 0000000000000000000000000000000000000000..d9d65cb8ede93b2b55859db60bca3036805c4008 --- /dev/null +++ b/tests/testthat/test-Compute-transform_indices.R @@ -0,0 +1,650 @@ +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( +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( +res1$output1, +res2$output1 +) +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 +) + +#----------------------------------- +# crop = FALSE +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 = FALSE), + 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( +res_crop_F_1 <- Compute(wf, chunks = list(lon = 2)) +) +suppressWarnings( +res_crop_F_2 <- Compute(wf, chunks = list(ensemble = 1)) +) +suppressWarnings( +res_crop_F_3 <- Compute(wf, chunks = list(lon = 3)) +) + +expect_equal( +res1$output1, +res_crop_F_1$output1 +) +expect_equal( +res_crop_F_1$output1, +res_crop_F_2$output1 +) +expect_equal( +res_crop_F_1$output1, +res_crop_F_3$output1 +) + +#--------------------------------------------- +#!!!!!!!!!!!!!!!!!!!!Problem when global + crop = T + chunk along lon!!!!!!!!!!!!!!!! +# crop = TRUE +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 = TRUE), + 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) + +#WRONG!!!!!!!!!! +#suppressWarnings( +#res_crop_T_1 <- Compute(wf, chunks = list(lon = 2)) +#) +suppressWarnings( +res_crop_T_2 <- Compute(wf, chunks = list(ensemble = 1)) +) +#WRONG!!!!!!!!!! +#suppressWarnings( +#res_crop_T_3 <- Compute(wf, chunks = list(lon = 3)) +#) +suppressWarnings( +res_crop_T_4 <- Compute(wf, chunks = list(lat = 2)) +) + +#expect_equal( +#res1$output1, +#res_crop_T_1$output1 +#) +expect_equal( +res1$output1, +res_crop_T_2$output1 +) +#expect_equal( +#res1$output1, +#res_crop_T_3$output1 +#) +expect_equal( +res1$output1, +res_crop_T_4$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', + crop = c(0, 22, -90, -60)), + 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( +res <- Compute(wf, chunks = list(lon = 2)) +) +suppressWarnings( +res2 <- Compute(wf, chunks = list(ensemble = 1)) +) +suppressWarnings( +res3 <- Compute(wf, chunks = list(lon = 3)) +) + + +expect_equal( +res$output1, +res2$output1 +) +expect_equal( +res$output1, +res3$output1 +) + +expect_equal( +drop(res$output1)[, 1], +c(241.5952, 243.0271, 247.6998, 246.7727, 248.7175, 267.7744, 273.2705), +tolerance = 0.001 +) +expect_equal( +drop(res$output1)[, 2], +c(241.4042, 242.5804, 246.8507, 245.8008, 246.4318, 267.0983, 272.9651), +tolerance = 0.001 +) +expect_equal( +drop(res$output1)[, 3], +c(241.2223, 242.2564, 245.9863, 244.5377, 244.8937, 266.5749, 272.5154), +tolerance = 0.001 +) +expect_equal( +drop(res$output1)[, 4], +c(241.0894, 242.1896, 245.3183, 243.1169, 243.9446, 266.4386, 272.4731), +tolerance = 0.001 +) +expect_equal( +drop(res$output1)[, 5], +c(241.0217, 242.3326, 244.6789, 241.6538, 244.3845, 266.6960, 272.4390), +tolerance = 0.001 +) + +#------------------------------------------------------ +# crop = FALSE +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', + crop = FALSE), + 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( +res_crop_F_1 <- Compute(wf, chunks = list(lon = 2)) +) +suppressWarnings( +res_crop_F_2 <- Compute(wf, chunks = list(ensemble = 1)) +) +suppressWarnings( +res_crop_F_3 <- Compute(wf, chunks = list(lon = 3)) +) + +expect_equal( +res$output1, +res_crop_F_1$output1 +) +expect_equal( +res$output1, +res_crop_F_2$output1 +) +expect_equal( +res$output1, +res_crop_F_3$output1 +) + + + +#------------------------------------------------------ +# crop = TRUE +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', + crop = T), + 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( +res_crop_T_1 <- Compute(wf, chunks = list(lat = 2, lon = 2)) +) +suppressWarnings( +res_crop_T_2 <- Compute(wf, chunks = list(ensemble = 1)) +) +suppressWarnings( +res_crop_T_3 <- Compute(wf, chunks = list(lon = 3)) +) + + +expect_equal( +res$output1, +res_crop_T_1$output1 +) +expect_equal( +res$output1, +res_crop_T_2$output1 +) +expect_equal( +res$output1, +res_crop_T_3$output1 +) + + +}) + +##################################################################### +##################################################################### +##################################################################### + +#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', + crop = c(0, 18, -90, -67)), + 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 <- Compute(wf, chunks = list(lon = 2)) +) +suppressWarnings( +res2 <- Compute(wf, chunks = list(ensemble = 1)) +) +suppressWarnings( +res3 <- Compute(wf, chunks = list(lon = 3)) +) + +expect_equal( +res1$output1, +res2$output1 +) +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 +) +expect_equal( +drop(res1$output1)[, 6], +c(241.0217, 242.3326, 244.6789, 241.6538, 244.3845, 266.6960), +tolerance = 0.001 +) + +#------------------------------------------------------ +# crop = FALSE +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', + crop = F), + 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( +res_crop_F_1 <- Compute(wf, chunks = list(lon = 2)) +) +suppressWarnings( +res_crop_F_2 <- Compute(wf, chunks = list(ensemble = 1)) +) +suppressWarnings( +res_crop_F_3 <- Compute(wf, chunks = list(lon = 3)) +) + +expect_equal( +as.vector(res1$output1), +as.vector(drop(res_crop_F_1$output1)[1:6, ]) +) +expect_equal( +res_crop_F_1$output1, +res_crop_F_2$output1 +) +expect_equal( +res_crop_F_1$output1, +res_crop_F_3$output1 +) + +#---------------------------------------------- +# crop = TRUE +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', + crop = T), + 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( +res_crop_T_1 <- Compute(wf, chunks = list(lon = 2)) +) +suppressWarnings( +res_crop_T_2 <- Compute(wf, chunks = list(ensemble = 1)) +) +suppressWarnings( +res_crop_T_3 <- Compute(wf, chunks = list(lon = 3)) +) + +expect_equal( +res_crop_F_1$output1, +res_crop_T_1$output1 +) +expect_equal( +res_crop_T_1$output1, +res_crop_T_2$output1 +) +expect_equal( +res_crop_T_1$output1, +res_crop_T_3$output1 +) + +}) diff --git a/tests/testthat/test-Compute-transform_values.R b/tests/testthat/test-Compute-transform_values.R new file mode 100644 index 0000000000000000000000000000000000000000..74975f7ca0904b46da0ec1608462894e9268cef2 --- /dev/null +++ b/tests/testthat/test-Compute-transform_values.R @@ -0,0 +1,924 @@ +context("Compute: Transform and chunk values()") +# Using values() to assign lat and lon, and transform the data. +# Also test transform + chunk along lat/lon. + +##################################################################### +##################################################################### +##################################################################### + +test_that("1. Global", { +skip_on_cran() + +lons.min <- 0 +lons.max <- 359.9 +lats.min <- -90 +lats.max <- 90 + +# crop = region +#NOTE: res1 and res3 differ if extra_cells = 2. But if retrieve = T, extra_cells = 2 or 8 is equal. + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res2 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res3 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +dim(res1), +c(sdate = 1, dat = 1, var = 1, ensemble = 1, time = 1, latitude = 50, longitude = 100) +) +expect_equal( +res1, +res2 +) +expect_equal( +res1, +res3 +) + +expect_equal( +drop(res1)[1:5, 1], +c(241.8592, 243.7243, 248.7337, 247.9308, 252.0744), +tolerance = 0.001 +) +expect_equal( +drop(res1)[23:28, 2], +c(298.0772, 299.4716, 299.7746, 300.2744, 300.3914, 299.5223), +tolerance = 0.001 +) +expect_equal( +mean(res1), +276.3901, +tolerance = 0.001 +) +#-------------------------------------------------- + +# crop = region, selector is values(c()) +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) + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(as.vector(lats)), + latitude_reorder = Sort(), + longitude = values(as.vector(lons)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res1_vector <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +expect_equal( +res1, +res1_vector +) + +#----------------------------------------------------------------- + +# crop = region, CircularSort(-180, 180) +lons.min <- -180 +lons.max <- 179.9 +lats.min <- -90 +lats.max <- 90 +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(-180, 180), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res1_180 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res2_180 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res3_180 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +as.vector(drop(res1)[, c(51:100, 1:50)]), +as.vector(res1_180), +tolerance = 0.0001 +) +expect_equal( +res1_180, +res2_180 +) +expect_equal( +res1_180, +res3_180 +) + +#============================================================ + +# crop = FALSE +lons.min <- 0 +lons.max <- 359.9 +lats.min <- -90 +lats.max <- 90 + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res_crop_F_1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res_crop_F_2 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res_crop_F_3 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +res1, +res_crop_F_1 +) +expect_equal( +res_crop_F_1, +res_crop_F_2 +) +expect_equal( +res_crop_F_1, +res_crop_F_3 +) + +#------------------------------------------------------- +#!!!!!!!!!!!!!!!!!!!!Problem when global + crop = T + chunk along lon!!!!!!!!!!!!!!!! + +# crop = TRUE +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = TRUE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +#WRONG +#suppressWarnings( +#res_crop_T_1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +#) +suppressWarnings( +res_crop_T_2 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +#WRONG +#suppressWarnings( +#res_crop_T_3 <- Compute(wf, chunks = list(longitude = 3))$output1 +#) +suppressWarnings( +res_crop_T_4 <- Compute(wf, chunks = list(latitude = 3))$output1 +) + +#expect_equal( +#res1, +#res_crop_T_1 +#) +expect_equal( +res1, +res_crop_T_2 +) +#expect_equal( +#res1, +#res_crop_T_3 +#) +expect_equal( +res1, +res_crop_T_4 +) + +}) + +############################################################################ +############################################################################ +############################################################################ + +# The numbers below are consistent with the result of this script. +#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 = 171:257, longitude = 30:81), +# vars_to_read = 'tas') +#lats <- NcToArray(file, +# dim_indices = list(latitude = 171:257), vars_to_read = 'latitude') +#lons <- NcToArray(file, +# dim_indices = list(longitude = 30:81), vars_to_read = 'longitude') +#NcClose(file) +# +#arr2 <- s2dv::CDORemap(arr, lons = as.vector(lons), lats = as.vector(lats), +# grid = 'r100x50', method = 'con', crop = c(10, 20, 20, 40)) + +test_that("2. Regional, no border", { + +skip_on_cran() + +lons.min <- 10 +lons.max <- 20 +lats.min <- 20 +lats.max <- 40 + +# crop = region +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', #paste0(2000:2001, '0101'), + ensemble = indices(1), #'all', + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, #'dat', + longitude = NULL, #'dat', + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res2 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res3 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +dim(res1), +c(sdate = 1, dat = 1, var = 1, ensemble = 1, time = 1, latitude = 5, longitude = 3) +) +expect_equal( +res1, +res2 +) +expect_equal( +res1, +res3 +) + +expect_equal( +drop(res1)[, 1], +c(284.9907, 282.9883, 281.2574, 284.1387, 285.6547), +tolerance = 0.0001 +) +expect_equal( +drop(res1)[, 2], +c(285.4820, 282.9362, 282.6088, 287.3716, 285.0194), +tolerance = 0.0001 +) +expect_equal( +drop(res1)[, 3], +c(286.1208, 284.3523, 285.9198, 287.7389, 286.1099), +tolerance = 0.0001 +) + +#------------------------------------------------------- + +# crop = FALSE +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res_crop_F_1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res_crop_F_2 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res_crop_F_3 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +res1, +res_crop_F_1 +) +expect_equal( +res_crop_F_1, +res_crop_F_2 +) +expect_equal( +res_crop_F_1, +res_crop_F_3 +) + +#------------------------------------------------------- + +# crop = TRUE +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = TRUE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res_crop_T_1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res_crop_T_2 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res_crop_T_3 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +res1, +res_crop_T_1 +) +expect_equal( +res_crop_T_1, +res_crop_T_2 +) +expect_equal( +res_crop_T_1, +res_crop_T_3 +) + +}) + +############################################################################ +############################################################################ +############################################################################ + +# The numbers below are consistent with the result of this script. +#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 = 171:257, longitude = c(1:81, 1289:1296)), +# vars_to_read = 'tas') +#lats <- NcToArray(file, +# dim_indices = list(latitude = 171:257), vars_to_read = 'latitude') +#lons <- NcToArray(file, +# dim_indices = list(longitude = c(1:81, 1289: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, 20, 20, 40)) + + +test_that("3. Regional, at lon border", { +skip_on_cran() + +lons.min <- 0 +lons.max <- 20 +lats.min <- 20 +lats.max <- 40 + +# crop = region +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res2 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res3 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +dim(res1), +c(sdate = 1, dat = 1, var = 1, ensemble = 1, time = 1, latitude = 5, longitude = 6) +) +expect_equal( +res1, +res2 +) +expect_equal( +res1, +res3 +) + +expect_equal( +drop(res1)[, 1], +c(286.4231, 283.8847, 280.4234, 277.7688, 284.3575), +tolerance = 0.001 +) +expect_equal( +drop(res1)[, 2], +c(288.0916, 283.9386, 280.9974, 278.4432, 284.8728), +tolerance = 0.001 +) +expect_equal( +drop(res1)[, 3], +c(285.7436, 283.1867, 281.7465, 280.2615, 284.6408), +tolerance = 0.001 +) +expect_equal( +drop(res1)[, 4], +c(284.9907, 282.9883, 281.2574, 284.1387, 285.6547), +tolerance = 0.001 +) +expect_equal( +drop(res1)[, 5], +c(285.4820, 282.9362, 282.6088, 287.3716, 285.0194), +tolerance = 0.001 +) +expect_equal( +drop(res1)[, 6], +c(286.1208, 284.3523, 285.9198, 287.7389, 286.1099), +tolerance = 0.001 +) + + +#-------------------------------------------------------------- + +# crop = region, CircularSort(-180, 180) +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(-180, 180), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res1_180 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res2_180 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res3_180 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +res1, +res1_180 +) +expect_equal( +res1_180, +res2_180 +) +expect_equal( +res1_180, +res3_180 +) + +#================================================================ + +# crop = FALSE +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res_crop_F_1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res_crop_F_2 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res_crop_F_3 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +res1, +res_crop_F_1 +) +expect_equal( +res_crop_F_1, +res_crop_F_2 +) +expect_equal( +res_crop_F_1, +res_crop_F_3 +) + +#------------------------------------------------------- + +# crop = FALSE, CircularSort(-180, 180) +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(-180, 180), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res_crop_F_1_180 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res_crop_F_2_180 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res_crop_F_3_180 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +res1, +res_crop_F_1_180 +) +expect_equal( +res_crop_F_1, +res_crop_F_2_180 +) +expect_equal( +res_crop_F_1, +res_crop_F_3_180 +) + +#=========================================================== + +# crop = TRUE +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = TRUE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res_crop_T_1 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res_crop_T_2 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res_crop_T_3 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +res1, +res_crop_T_1 +) +expect_equal( +res_crop_T_1, +res_crop_T_2 +) +expect_equal( +res_crop_T_1, +res_crop_T_3 +) +#-------------------------------------------------- +# crop = TRUE, CircularSort(-180, 180) +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(-180, 180), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = TRUE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = F) +) +func <- function(exp) { + return(exp) +} +step <- Step(func, + target_dims = 'sdate', output_dims = 'sdate') +wf <- AddStep(exp, step) +suppressWarnings( +res_crop_T_1_180 <- Compute(wf, chunks = list(latitude = 2, longitude = 2))$output1 +) +suppressWarnings( +res_crop_T_2_180 <- Compute(wf, chunks = list(ensemble = 1))$output1 +) +suppressWarnings( +res_crop_T_3_180 <- Compute(wf, chunks = list(longitude = 3))$output1 +) + +expect_equal( +res1, +res_crop_T_1_180 +) +expect_equal( +res_crop_T_1, +res_crop_T_2_180 +) +expect_equal( +res_crop_T_1, +res_crop_T_3_180 +) + +}) diff --git a/tests/testthat/test-Start-indices_list_vector.R b/tests/testthat/test-Start-indices_list_vector.R new file mode 100644 index 0000000000000000000000000000000000000000..39ecb24f15edde38460a02d06095e59f612dae19 --- /dev/null +++ b/tests/testthat/test-Start-indices_list_vector.R @@ -0,0 +1,241 @@ +# This unit test tests the consistence between list of indices and vector of indices. +# 1. transform +# 2. no transform +# 3. transform, indices reversed +# 4. no transform, indices reversed + +context("List of indices and vector of indices") + + +test_that("1. transform", { + +# lat and lon are lists of indices +suppressWarnings( +exp1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(list(1, 30)), + latitude_reorder = Sort(), + longitude = indices(list(1, 40)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(0, 11, -90, -81)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +# lat and lon are vectors of indices +suppressWarnings( +exp2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 1:30, + latitude_reorder = Sort(), + longitude = 1:40, + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(0, 11, -90, -81)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +expect_equal( +as.vector(exp1), +as.vector(exp2) +) + +}) + +############################################################# +############################################################# +############################################################# + +test_that("2. no transform", { + +# lat and lon are lists of indices +suppressWarnings( +exp1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(list(1, 30)), + latitude_reorder = Sort(), + longitude = indices(list(1, 40)), + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +# lat and lon are vectors of indices +suppressWarnings( +exp2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 1:30, + latitude_reorder = Sort(), + longitude = 1:40, + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +expect_equal( +as.vector(exp1), +as.vector(exp2) +) + +}) + + +############################################################# +############################################################# +############################################################# +# PROBLEM: +# latitude is -81 to -88.2 now, but it should be 81 to 88.2 because the indices is retrieved first then do the transform (aiat = F); and it should be ascending. + +# .. ..$ latitude : num [1:3(1d)] -81 -84.6 -88.2 + +test_that("3. transform, indices reverse", { + +# lat and lon are lists of indices +suppressWarnings( +exp1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(list(30, 1)), + latitude_reorder = Sort(), + longitude = indices(list(1, 40)), # can't reverse. Different meaning + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(0, 11, -90, -81)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +# lat and lon are vectors of indices +suppressWarnings( +exp2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 30:1, + latitude_reorder = Sort(), + longitude = 40:1, + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(0, 11, -90, -81)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +expect_equal( +as.vector(drop(exp1)[, 4:1]), +as.vector(exp2) +) + +}) + +################################################################ +################################################################ +################################################################ + +test_that("4. no transform, indices reverse", { + +# lat and lon are lists of indices +suppressWarnings( +exp1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(list(30, 1)), + latitude_var = 'latitude', + latitude_reorder = Sort(), + longitude = indices(list(1, 40)), # can't reverse. different meaning + longitude_var = 'longitude', + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +# lat and lon are vectors of indices +suppressWarnings( +exp2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 30:1, + latitude_var = 'latitude', + latitude_reorder = Sort(), + longitude = 40:1, + longitude_var = 'longitude', + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +expect_equal( +as.vector(drop(exp1)[, 40:1]), +as.vector(exp2) +) + +}) diff --git a/tests/testthat/test-Start-metadata_dims.R b/tests/testthat/test-Start-metadata_dims.R index 841786dcef6c68df19554108965e801c97abf41b..3239e7de37b878030d776ea5019d15ceacf43fe6 100644 --- a/tests/testthat/test-Start-metadata_dims.R +++ b/tests/testthat/test-Start-metadata_dims.R @@ -406,7 +406,7 @@ suppressWarnings( ) expect_equal( attr(data, "Variables")$system4_m1, - NULL + list(lon = NULL, lat = NULL) ) expect_equal( length(attr(data, "Variables")$system5_m1$lon), diff --git a/tests/testthat/test-Start-reorder-lat.R b/tests/testthat/test-Start-reorder-lat.R index e2ef5d9d00a99305780d770dbfa43bcf2ff88482..4133cf01e4b357350bb85dfc58212c472afc8060 100644 --- a/tests/testthat/test-Start-reorder-lat.R +++ b/tests/testthat/test-Start-reorder-lat.R @@ -880,7 +880,175 @@ test_that("1-4. Selector type: indices(vector)", { }) ############################################## -test_that("1-4. Selector type: indices(vector)", { +test_that("4-x-2-12-123-2-1-x", { + +# 1-1. no Sort(), NULL +## lat should be descending +suppressWarnings( +exp1_1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', +# latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('latitude', 'lat'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +# 1-2. Sort(), NULL +## lat should be ascending +suppressWarnings( +exp1_2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('latitude', 'lat'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +# 1-3. Sort(drcreasing = T), NULL +## lat should be descending +suppressWarnings( +exp1_3 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(decreasing = T), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('latitude', 'lat'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +expect_equal( +drop(exp1_1)[1:5, 2], +c(250.8470, 251.0054, 251.1874, 251.3769, 251.5602), +tolerance = 0.0001 +) +expect_equal( +as.vector(drop(exp1_1)), +as.vector(drop(exp1_2)[640:1, ]) +) +expect_equal( +as.vector(drop(exp1_1)), +as.vector(drop(exp1_3)) +) +expect_equal( +as.vector(attr(exp1_1, 'Variables')$common$latitude)[1:5], +c(89.78488, 89.50620, 89.22588, 88.94519, 88.66436), +tolerance = 0.0001 +) +expect_equal( +as.vector(attr(exp1_2, 'Variables')$common$latitude)[1:5], +c(-89.78488, -89.50620, -89.22588, -88.94519, -88.66436), +tolerance = 0.0001 +) +expect_equal( +as.vector(attr(exp1_1, 'Variables')$common$latitude), +as.vector(attr(exp1_3, 'Variables')$common$latitude) +) -}) +# 2-1. no Sort(), 'dat' +## lat should be descending +suppressWarnings( +exp2_1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', +# latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('latitude', 'lat'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) +# 2-2. Sort(), 'dat' +## lat should be ascending +suppressWarnings( +exp2_2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('latitude', 'lat'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) +# 2-3. Sort(drcreasing = T), NULL +## lat should be descending +suppressWarnings( +exp2_3 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(decreasing = T), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('latitude', 'lat'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) +expect_equal( +as.vector(drop(exp1_1)), +as.vector(drop(exp2_1)) +) +expect_equal( +as.vector(drop(exp1_2)), +as.vector(drop(exp2_2)) +) +expect_equal( +as.vector(drop(exp1_3)), +as.vector(drop(exp2_3)) +) +expect_equal( +as.vector(attr(exp2_1, 'Variables')$dat1$latitude), +as.vector(attr(exp1_1, 'Variables')$common$latitude) +) +expect_equal( +as.vector(attr(exp2_2, 'Variables')$dat1$latitude), +as.vector(attr(exp1_2, 'Variables')$common$latitude) +) +expect_equal( +as.vector(attr(exp2_3, 'Variables')$dat1$latitude), +as.vector(attr(exp1_3, 'Variables')$common$latitude) +) + +}) +############################################## diff --git a/tests/testthat/test-Start-transform-border.R b/tests/testthat/test-Start-transform-border.R new file mode 100644 index 0000000000000000000000000000000000000000..a84f1eca315b5929cc57da45ba3f305bbc6cdbb9 --- /dev/null +++ b/tests/testthat/test-Start-transform-border.R @@ -0,0 +1,501 @@ +context("Transform: check with cdo") +# This unit test checks different border situations: normal regional that doesn't touch the borders, +# global situation that uses all the grids, or one side reaches the border. + +# Compare the results with cdo. The example script is as below: +#library(easyNCDF) +#path <- "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc" +#file <- NcOpen(path) +#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) +# +#dim(arr) +#dim(lats) +#arr2 <- s2dv::CDORemap(arr, lons = as.vector(lons), lats = as.vector(lats), +# grid = 'r100x50', method = 'con', crop = c(10,20,20,40)) +#dim(arr2$data_array) + +# The result of cdo is from CDO version 1.9.8. + +test_that("1. normal regional situation", { + +lons.min <- 10 +lons.max <- 20 +lats.min <- 20 +lats.max <- 40 + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +expect_equal( +dim(drop(exp)), +c(latitude = 5, longitude = 3) +) +expect_equal( +drop(exp)[, 1], +c(284.9907, 282.9883, 281.2574, 284.1387, 285.6547), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 2], +c(285.4820, 282.9362, 282.6088, 287.3716, 285.0194), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 3], +c(286.1208, 284.3523, 285.9198, 287.7389, 286.1099), +tolerance = 0.0001 +) + +}) + +# The result is consistent with cdo result, using all the grid points to transform +# and crop the region. +# drop(exp)[, ] +# [,1] [,2] [,3] +#[1,] 284.9907 285.4820 286.1208 +#[2,] 282.9883 282.9362 284.3523 +#[3,] 281.2574 282.6088 285.9198 +#[4,] 284.1387 287.3716 287.7389 +#[5,] 285.6547 285.0194 286.1099 + +#------------------------------------------------ + +test_that("2. global situation", { + +lons.min <- 0 +lons.max <- 359.9 +lats.min <- 20 +lats.max <- 40 + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +expect_equal( +dim(drop(exp)), +c(latitude = 5, longitude = 100) +) +expect_equal( +drop(exp)[, 1], +c(286.4231, 283.8847, 280.4234, 277.7688, 284.3575), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 2], +c(288.0916, 283.9386, 280.9974, 278.4432, 284.8728), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 100], +c(285.9373, 283.6340, 280.6685, 279.6016, 279.5081), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 99], +c(286.3033, 283.8651, 279.9846, 285.0679, 282.6013), +tolerance = 0.0001 +) + +}) + +#> drop(arr2$data_array)[,1:3] +# [,1] [,2] [,3] +#[1,] 286.4231 288.0916 285.7435 +#[2,] 283.8847 283.9386 283.1867 +#[3,] 280.4234 280.9974 281.7465 +#[4,] 277.7688 278.4432 280.2615 +#[5,] 284.3575 284.8728 284.6408 +#> drop(arr2$data_array)[,98:100] +# [,1] [,2] [,3] +#[1,] 286.4648 286.3033 285.9373 +#[2,] 285.5226 283.8651 283.6340 +#[3,] 287.8567 279.9846 280.6685 +#[4,] 288.6723 285.0679 279.6016 +#[5,] 286.8253 282.6013 279.5081 + +#----------------------------------------------- + +test_that("3. left side too close to border", { + +lons.min <- 0 +lons.max <- 20 +lats.min <- 20 +lats.max <- 40 + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +expect_equal( +dim(drop(exp)), +c(latitude = 5, longitude = 6) +) +expect_equal( +drop(exp)[, 1], +c(286.4231, 283.8847, 280.4234, 277.7688, 284.3575), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 2], +c(288.0916, 283.9386, 280.9974, 278.4432, 284.8728), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 3], +c(285.7436, 283.1867, 281.7465, 280.2615, 284.6408), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 4], +c(284.9907, 282.9883, 281.2574, 284.1387, 285.6547), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 5], +c(285.4820, 282.9362, 282.6088, 287.3716, 285.0194), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 6], +c(286.1208, 284.3523, 285.9198, 287.7389, 286.1099), +tolerance = 0.0001 +) + +}) + + +#drop(arr2$data_array)[,] +# [,1] [,2] [,3] [,4] [,5] [,6] +#[1,] 286.4231 288.0916 285.7435 284.9907 285.4820 286.1208 +#[2,] 283.8847 283.9386 283.1867 282.9882 282.9362 284.3523 +#[3,] 280.4234 280.9974 281.7465 281.2574 282.6088 285.9198 +#[4,] 277.7688 278.4432 280.2615 284.1387 287.3716 287.7389 +#[5,] 284.3575 284.8728 284.6408 285.6547 285.0194 286.1099 + + + +test_that("4. right side too close to border", { + +lons.min <- 350 +lons.max <- 359 +lats.min <- 20 +lats.max <- 40 + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +expect_equal( +dim(drop(exp)), +c(latitude = 5, longitude = 2) +) +expect_equal( +drop(exp)[, 1], +c(286.3033, 283.8651, 279.9846, 285.0679, 282.6013), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 2], +c(285.9373, 283.6340, 280.6685, 279.6016, 279.5081), +tolerance = 0.0001 +) + +}) + +#> drop(arr2$data_array)[,] +# [,1] [,2] +#[1,] 286.3033 285.9373 +#[2,] 283.8651 283.6340 +#[3,] 279.9846 280.6685 +#[4,] 285.0679 279.6016 +#[5,] 282.6013 279.5081 + +#-------------------------------------------------- + +test_that("5. across meridian", { + +lons.min <- 170 +lons.max <- 190 +lats.min <- 20 +lats.max <- 40 + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(-180, 180), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +expect_equal( +dim(drop(exp)), +c(latitude = 5, longitude = 5) +) +expect_equal( +drop(exp)[, 1], +c(295.9371, 294.0865, 291.8104, 289.0014, 284.9630), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 2], +c(296.2407, 294.4130, 291.8895, 289.5334, 285.7766), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 3], +c(296.2065, 294.4305, 291.9352, 289.5931, 286.0924), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 4], +c(295.7689, 293.6672, 291.2874, 288.4160, 284.6429), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 5], +c(295.8033, 293.8032, 291.5909, 288.5543, 284.6293), +tolerance = 0.0001 +) + +}) + +# cdo result is [170:190]; Start() is [-180:-170; 170:180]. +# So drop(exp)[, 1:3] == drop(arr2$data_array)[, 3:5]. +#drop(arr2$data_array)[,] +# [,1] [,2] [,3] [,4] [,5] +#[1,] 295.7689 295.8034 295.9371 296.2407 296.2065 +#[2,] 293.6672 293.8032 294.0865 294.4130 294.4306 +#[3,] 291.2874 291.5910 291.8104 291.8895 291.9352 +#[4,] 288.4159 288.5543 289.0014 289.5334 289.5931 +#[5,] 284.6429 284.6293 284.9630 285.7766 286.0924 + + +test_that("6. normal case; [-180, 180]", { +# The lon range is too close to border for the original longitude [0, 360], but +# is normal case for [-180, 180]. In zzz.R, it is counted as normal case, and the +# result is the same as 3. +lons.min <- 0 +lons.max <- 20 +lats.min <- 20 +lats.max <- 40 + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(-180, 180), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +expect_equal( +dim(drop(exp)), +c(latitude = 5, longitude = 6) +) +expect_equal( +drop(exp)[, 1], +c(286.4231, 283.8847, 280.4234, 277.7688, 284.3575), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 2], +c(288.0916, 283.9386, 280.9974, 278.4432, 284.8728), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 3], +c(285.7436, 283.1867, 281.7465, 280.2615, 284.6408), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 4], +c(284.9907, 282.9883, 281.2574, 284.1387, 285.6547), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 5], +c(285.4820, 282.9362, 282.6088, 287.3716, 285.0194), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 6], +c(286.1208, 284.3523, 285.9198, 287.7389, 286.1099), +tolerance = 0.0001 +) + +}) + +#---------------------------------------------------------- + +test_that("6. left side too close to border; [-180, 180]", { +# The lon range is too close to border for the original longitude [0, 360], but +# is normal case for [-180, 180]. In zzz.R, it is counted as normal case, and the +# result is the same as 3. +lons.min <- -179 +lons.max <- -170 +lats.min <- 20 +lats.max <- 40 + +suppressWarnings( +exp <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(lats.min, lats.max)), + latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = CircularSort(-180, 180), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(lons.min, lons.max, lats.min, lats.max)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve = T) +) +expect_equal( +dim(drop(exp)), +c(latitude = 5, longitude = 2) +) +expect_equal( +drop(exp)[, 1], +c(296.2407, 294.4130, 291.8895, 289.5334, 285.7766), +tolerance = 0.0001 +) +expect_equal( +drop(exp)[, 2], +c(296.2065, 294.4305, 291.9352, 289.5931, 286.0924), +tolerance = 0.0001 +) + +}) + +# cdo result is [181:190] +#> drop(arr2$data_array)[,] +# [,1] [,2] +#[1,] 296.2407 296.2065 +#[2,] 294.4130 294.4306 +#[3,] 291.8895 291.9352 +#[4,] 289.5334 289.5931 +#[5,] 285.7766 286.0924 + diff --git a/tests/testthat/test-Start-transform-lat-Sort-all.R b/tests/testthat/test-Start-transform-lat-Sort-all.R new file mode 100644 index 0000000000000000000000000000000000000000..28393699cd5c4e15194963f60135ee57a9e99b57 --- /dev/null +++ b/tests/testthat/test-Start-transform-lat-Sort-all.R @@ -0,0 +1,125 @@ +# This unit test uses 'all' to do the transformation and tests "lat_reorder". +# The results should be identical and consistent with cdo result (with precision difference). +# "lon_reorder = CircularSort(0, 360)" are used in all the tests. +# The test contains three calls: lat_reorder = Sort(), no lat_reorder, and lat_reorder = Sort(decreasing = T). +# Note that the original latitude is descending [90:-90]. cdo result is ascending [-90:90]. + +context("Transform and lat_reorder test: 'all'") + +#--------------------------------------------------------------- +# cdo is used to verify the data values +library(easyNCDF) +path <- "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc" +file <- NcOpen(path) +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)$data_array +#--------------------------------------------------------------- + +path <- '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc' + +test_that("1. 'all'", { + +# lat should be ascending +suppressWarnings( +res1 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending or ascending? Because Sort() is not specified and 'all' does not +# say the order either, it could follow the transformed order (if so, ascending). +suppressWarnings( +res2 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', +# latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, #'dat', + longitude = NULL, #'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res3 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(decreasing = T), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +expect_equal( +as.vector(res1), +as.vector(arr2), +tolerance = 0.0001 +) +expect_equal( +as.vector(res2), +as.vector(arr2), +tolerance = 0.0001 +) + +expect_equal( +as.vector(drop(res3)[50:1, ]), +as.vector(arr2), +tolerance = 0.0001 +) + +}) + diff --git a/tests/testthat/test-Start-transform-lat-Sort-indices.R b/tests/testthat/test-Start-transform-lat-Sort-indices.R new file mode 100644 index 0000000000000000000000000000000000000000..1a1f1ee924f0264c1371cafc33689ed0a29135eb --- /dev/null +++ b/tests/testthat/test-Start-transform-lat-Sort-indices.R @@ -0,0 +1,230 @@ +# This unit test uses indices() to do the transformation and tests "lat_reorder". +# The results should be identical and consistent with cdo result (with precision difference). +# The lat/lon range is all the grids here. +# "lon_reorder = CircularSort(0, 360)" are used in all the tests. +# Test 1 uses indices(1:640), and test 2 uses indices(640:1). +# Each of them contains three calls: lat_reorder = Sort(), no lat_reorder, and lat_reorder = Sort(decreasing = T). +# Note that the original latitude is descending [90:-90]. cdo result is ascending [-90:90]. + +#!!!!!!!!!!!!!!!!!!!!!PROBLEM in test 2, indices(640:1)!!!!!!!!!!!!!!!!!!!! +#TODO: Add regional test + +context("Transform and lat_reorder test: indices") + +#--------------------------------------------------------------- +# cdo is used to verify the data values +library(easyNCDF) +path <- "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc" +file <- NcOpen(path) +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)$data_array +#--------------------------------------------------------------- + +path <- '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc' + +test_that("1. indices(1:640)", { + +# lat should be ascending +suppressWarnings( +res1 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(1:640), + latitude_reorder = Sort(), + longitude = indices(1:1296), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res2 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(1:640), +# latitude_reorder = Sort(), + longitude = indices(1:1296), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res3 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(1:640), + latitude_reorder = Sort(decreasing = T), + longitude = indices(1:1296), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +expect_equal( +as.vector(res1), +as.vector(arr2), +tolerance = 0.0001 +) +expect_equal( +as.vector(drop(res2)[50:1, ]), +as.vector(arr2), +tolerance = 0.0001 +) +expect_equal( +as.vector(drop(res3)[50:1, ]), +as.vector(arr2), +tolerance = 0.0001 +) + +}) + +################################################### +################################################### +################################################### + +test_that("2. indices(640:1)", { + +# lat should be ascending +suppressWarnings( +res1 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(640:1), + latitude_reorder = Sort(), + longitude = indices(1:1296), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be ascending +suppressWarnings( +res2 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(640:1), +# latitude_reorder = Sort(), + longitude = indices(1:1296), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res3 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(640:1), + latitude_reorder = Sort(decreasing = T), + longitude = indices(1:1296), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +#WRONG!!!!!!!!!! it is descending now +#expect_equal( +#as.vector(res1), +#as.vector(arr2), +#tolerance = 0.0001 +#) + +expect_equal( +as.vector(res2), +as.vector(arr2), +tolerance = 0.0001 +) +#WRONG!!!!!!!!!! it is ascending now +#expect_equal( +#as.vector(drop(res3)[50:1, ]), +#as.vector(arr2), +#tolerance = 0.0001 +#) + +}) + diff --git a/tests/testthat/test-Start-transform-lat-Sort-values.R b/tests/testthat/test-Start-transform-lat-Sort-values.R new file mode 100644 index 0000000000000000000000000000000000000000..af00f73095106a416186ea4efa781aac0055f292 --- /dev/null +++ b/tests/testthat/test-Start-transform-lat-Sort-values.R @@ -0,0 +1,430 @@ +# This unit test uses values() to do the transformation and tests "lat_reorder". +# The results should be identical and consistent with cdo result (with precision difference). +# The lon range is all the grids here. +# "lon_reorder = CircularSort(0, 360)" are used in all the tests. +# Test 1 & 2 are global: test 1 uses values(list(-90, 90)) and test 2 uses values(list(90, -90)). +# Test 3 & 4 are regional: test 3 uses values(list(-90, -80)) and test 4 uses values(list(-80, -90)). +# Each of them contains three calls: lat_reorder = Sort(), no lat_reorder, and lat_reorder = Sort(decreasing = T). +# Note that the original latitude is descending [90:-90]. cdo result is ascending [-90:90]. + +context("Transform and lat_reorder test: values") + +#--------------------------------------------------------------- +# cdo is used to verify the data values +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)$data_array +#--------------------------------------------------------------- + +path <- '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc' + +test_that("1. values(list(-90, 90))", { + +# lat should be ascending +suppressWarnings( +res1 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-90, 90)), + latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be ascending +suppressWarnings( +res2 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-90, 90)), +# latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res3 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-90, 90)), + latitude_reorder = Sort(decreasing = T), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +expect_equal( +as.vector(res1), +as.vector(arr2), +tolerance = 0.0001 +) +expect_equal( +as.vector(res2), +as.vector(arr2), +tolerance = 0.0001 +) +expect_equal( +as.vector(drop(res3)[50:1, ]), +as.vector(arr2), +tolerance = 0.0001 +) + +}) + +############################################################ +############################################################ +############################################################ + +test_that("2. values(list(90, -90))", { + +# lat should be ascending +suppressWarnings( +res1 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(90, -90)), + latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res2 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(90, -90)), +# latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res3 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(90, -90)), + latitude_reorder = Sort(decreasing = T), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + + +expect_equal( +as.vector(res1), +as.vector(arr2), +tolerance = 0.0001 +) +expect_equal( +as.vector(drop(res2)[50:1, ]), +as.vector(arr2), +tolerance = 0.0001 +) +expect_equal( +as.vector(drop(res3)[50:1, ]), +as.vector(arr2), +tolerance = 0.0001 +) + +}) + + +############################################################ +############################################################ +############################################################ + +#NOTE: The numbers at lat = 3 are different with cdo if transform_extra_cells = 2. + +test_that("3. values(list(-90, -80))", { + +# lat should be ascending +suppressWarnings( +res1 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-90, -80)), + latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = T), # note that crop = T here + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be ascending +suppressWarnings( +res2 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-90, -80)), +# latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = TRUE), #FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res3 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-90, -80)), + latitude_reorder = Sort(decreasing = T), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = TRUE), #FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +expect_equal( +as.vector(res1), +as.vector(drop(arr2)[1:3, ]), +tolerance = 0.0001 +) +expect_equal( +as.vector(res2), +as.vector(drop(arr2)[1:3, ]), +tolerance = 0.0001 +) +expect_equal( +as.vector(res3), +as.vector(drop(arr2)[3:1, ]), +tolerance = 0.0001 +) + +}) + +############################################################ +############################################################ +############################################################ + + +#NOTE: The numbers at lat = 3 are different with cdo if transform_extra_cells = 2. + +test_that("4. values(list(-80, -90))", { + +# lat should be ascending +suppressWarnings( +res1 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-80, -90)), + latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = T), # note that crop = T here + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res2 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-80, -90)), +# latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = TRUE), #FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +# lat should be descending +suppressWarnings( +res3 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-80, -90)), + latitude_reorder = Sort(decreasing = T), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = TRUE), #FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +expect_equal( +as.vector(res1), +as.vector(drop(arr2)[1:3, ]), +tolerance = 0.0001 +) +expect_equal( +as.vector(res2), +as.vector(drop(arr2)[3:1, ]), +tolerance = 0.0001 +) +expect_equal( +as.vector(res3), +as.vector(drop(arr2)[3:1, ]), +tolerance = 0.0001 +) + +}) + diff --git a/tests/testthat/test-Start-transform-three-selectors.R b/tests/testthat/test-Start-transform-three-selectors.R new file mode 100644 index 0000000000000000000000000000000000000000..3eb00409ba6f0a65fc91943c026af9f57cacf592 --- /dev/null +++ b/tests/testthat/test-Start-transform-three-selectors.R @@ -0,0 +1,194 @@ +# This unit test uses three different selector forms: indices(), values(), and 'all', to do +# the transformation. "lat_reorder" is also tested. +# Their results should be all identical and consistent with cdo result (with precision difference). +# The selected lat/lon range is all the grids here. +# "lon_reorder = CircularSort(0, 360)" and "lat = Sort()" are used in all the tests. +# To see different lat_reorder options, go to "test-Start-transform-lat-Sort-*". +# If values, the lat selector is [-90, 90] or [90, -90]; if indices, c(1:640) or c(640:1). + +# Note that the original latitude is descending [90:-90]. + +context("Transform: three selector forms") + +#--------------------------------------------------------------- +# cdo is used to verify the data values +library(easyNCDF) +path <- "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc" +file <- NcOpen(path) +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)$data_array +#--------------------------------------------------------------- + +path <- '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc' + +test_that("1. indices", { + +suppressWarnings( +res1 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(1:640), + latitude_reorder = Sort(), + longitude = indices(1:1296), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + + +suppressWarnings( +res2 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = indices(640:1), + latitude_reorder = Sort(), + longitude = indices(1:1296), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +expect_equal( +as.vector(res1), +as.vector(arr2), +tolerance = 0.0001 +) +#WRONG!!!!! res2 lat is descending now +#expect_equal( +#as.vector(res1), +#as.vector(res2) +#) + +}) + + +test_that("2. values", { + +suppressWarnings( +res1 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(-90, 90)), + latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) +suppressWarnings( +res2 <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(90, -90)), + latitude_reorder = Sort(), + longitude = values(list(0, 359.9)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +expect_equal( +as.vector(res1), +as.vector(arr2), +tolerance = 0.0001 +) +expect_equal( +as.vector(res2), +as.vector(arr2), +tolerance = 0.0001 +) + + +}) + + +test_that("3. all", { + +suppressWarnings( +res <- Start(dat = path, + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = 'all', + latitude_reorder = Sort(), + longitude = 'all', + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = FALSE), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 2, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = 'dat', + longitude = 'dat', + time = 'sdate'), + retrieve = T) +) + +expect_equal( +as.vector(res), +as.vector(arr2), +tolerance = 0.0001 +) + +}) diff --git a/tests/testthat/test-Start-values_list_vector.R b/tests/testthat/test-Start-values_list_vector.R new file mode 100644 index 0000000000000000000000000000000000000000..d93b2473526a0237e59b6174902ab6761513deab --- /dev/null +++ b/tests/testthat/test-Start-values_list_vector.R @@ -0,0 +1,245 @@ +# This unit test tests the consistence between list of values and vector of values. +# 1. transform +# 2. no transform +# 3. transform, indices reversed +# 4. no transform, indices reversed + +context("List of values and vector of values") + +#----------------------------------------------------------------- +# To get lat and lon vectors +library(easyNCDF) +pathh <- "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20000101.nc" +file <- NcOpen(pathh) +lats <- NcToArray(file, + dim_indices = list(latitude = 1:35), vars_to_read = 'latitude') +lons <- NcToArray(file, + dim_indices = list(longitude = 1:33), vars_to_read = 'longitude') +NcClose(file) +#------------------------------------------------------------------ + +test_that("1. transform", { + +# lat and lon are lists of values +suppressWarnings( +exp1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(81, 88.2)), + latitude_reorder = Sort(), + longitude = values(list(0, 7.2)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(0, 9, 80, 90)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +# lat and lon are vectors of values. This one is a weird usage though... +suppressWarnings( +exp2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(c(81, 84.6, 88.2)), + latitude_reorder = Sort(), + longitude = values(c(0, 3.6, 7.2)), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(0, 9, 80, 90)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +expect_equal( +as.vector(exp1), +as.vector(exp2) +) + +}) + +############################################################# +############################################################# +############################################################# + +test_that("2. no transform", { + +# lat and lon are lists of indices +suppressWarnings( +exp1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(80, 90)), + latitude_reorder = Sort(), + longitude = values(list(0, 9)), + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +# lat and lon are vectors of indices +suppressWarnings( +exp2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(as.vector(lats)), + latitude_reorder = Sort(), + longitude = values(as.vector(lons)), + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +expect_equal( +as.vector(exp1), +as.vector(exp2) +) + +}) + + +############################################################# +############################################################# +############################################################# + +test_that("3. transform, vector reverse", { + +# lat and lon are lists of values +suppressWarnings( +exp1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(88.2, 81)), + latitude_reorder = Sort(), + longitude = values(list(0, 7.2)), # It can't be reversed; different meanings + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(0, 9, 80, 90)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +#WRONG!!!!!!!!!! +# lat and lon are vectors of values +suppressWarnings( +exp2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(rev(c(81, 84.6, 88.2))), + latitude_reorder = Sort(), + longitude = values(rev(c(0, 3.6, 7.2))), + longitude_reorder = CircularSort(0, 360), + transform = CDORemapper, + transform_params = list(grid = 'r100x50', + method = 'con', + crop = c(0, 9, 80, 90)), + transform_vars = c('latitude', 'longitude'), + transform_extra_cells = 8, + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +expect_equal( +as.vector(exp1), +as.vector(exp2) +) + +}) + +################################################################ +################################################################ +################################################################ + +test_that("4. no transform, vector reverse", { + +# lat and lon are lists of values +suppressWarnings( +exp1 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(list(90, 80)), + latitude_reorder = Sort(), + longitude = values(list(0, 9)), # it can't be reversed; different meanings + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +# lat and lon are vectors of values +suppressWarnings( +exp2 <- Start(dat = '/esarchive/exp/ecmwf/system5_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc', + var = 'tas', + sdate = '20000101', + ensemble = indices(1), + time = indices(1), + latitude = values(rev(as.vector(lats))), + latitude_reorder = Sort(), + longitude = values(rev(as.vector(lons))), + longitude_reorder = CircularSort(0, 360), + synonims = list(latitude = c('lat', 'latitude'), + longitude = c('longitude', 'lon')), + return_vars = list(latitude = NULL, + longitude = NULL, + time = 'sdate'), + retrieve= T) +) + +expect_equal( +as.vector(exp1), +as.vector(exp2) +) + +})