diff --git a/DESCRIPTION b/DESCRIPTION index cf50cf03ac3c494a290ad886fc18c5d9fc6431cf..d36c898a2b20a453e41201f55896a3c82852da87 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: s2dverification Title: Set of Common Tools for Forecast Verification -Version: 2.8.1 +Version: 2.8.2 Authors@R: c( person("BSC-CNS", role = c("aut", "cph")), person("Virginie", "Guemas", , "virginie.guemas@bsc.es", role = "aut"), diff --git a/R/ArrayToNetCDF.R b/R/ArrayToNetCDF.R index 2b573f639be3db0d7daadd2f278efbe4d5b50ff7..47fe13ef20ee9003cc937210915f820277411dea 100644 --- a/R/ArrayToNetCDF.R +++ b/R/ArrayToNetCDF.R @@ -40,7 +40,7 @@ ArrayToNetCDF <- function(arrays, file_path) { } dim_names <- names(dim(arrays[[i]])) if (!is.null(dim_names)) { - if (any(is.na(dim_names) || (sapply(dim_names, nchar) == 0))) { + if (any(is.na(dim_names) | (sapply(dim_names, nchar) == 0))) { stop("The provided arrays must have all named dimensions or ", "all unnamed dimensions.") } @@ -74,7 +74,7 @@ ArrayToNetCDF <- function(arrays, file_path) { if (!is.numeric(dim_info[['len']])) { stop("The provided 'len' for the ", k, "th dimension in the ", i, "th array must be a numeric value.") } - dim_info[['len']] <- round(dim_info[['len']][1]) + dim_info[['len']] <- as.integer(round(dim_info[['len']][1])) if (dim_info[['len']] != dim(arrays[[i]])[k]) { stop("The provided 'len' for the ", k, "th dimension in the ", i, "th array does not match the actual length of the provided array.") } @@ -98,7 +98,7 @@ ArrayToNetCDF <- function(arrays, file_path) { if (!('vals' %in% names(dim_info))) { dim_info[['vals']] <- 1:dim_info[['len']] } else { - if (!is.numeric(dim_info[['vals']])) { + if (!(is.numeric(dim_info[['vals']]))) { stop("The provided 'vals' for the ", k, "th dimension in the ", i, "th array must be a numeric vector.") } if (dim_info[['units']] == '') { @@ -292,10 +292,25 @@ ArrayToNetCDF <- function(arrays, file_path) { if (!is.character(var_info[['coordinates']])) { stop("The attribute 'coordinates' must be a character string.") } - if (!(all(strsplit(var_info[['coordinates']], ' ')[[1]] %in% sapply(defined_vars, '[[', 'name')))) { - stop("All the dimensions appearing in 'coordinates' must point to defined variables.") + coords <- strsplit(var_info[['coordinates']], ' ')[[1]] + if (!(all(coords %in% sapply(defined_vars, '[[', 'name') | + coords %in% sapply(defined_dims[which(sapply(defined_dims, '[[', 'create_dimvar'))], '[[', 'name')))) { + coords <- coords[which(coords %in% sapply(defined_vars, '[[', 'name') | + coords %in% sapply(defined_dims[which(sapply(defined_dims, '[[', 'create_dimvar'))], '[[', 'name'))] + .warning("Some of the dimensions appearing in 'coordinates' have been removed because they point to undefined variables.") + } + ncatt_put(ncdf_object, defined_vars[[var_counter]]$name, 'coordinates', paste(coords, collapse = ' ')) + } + attrs_to_skip <- which(names(var_info) %in% c('addOffset', 'scaleFact', 'coordinates', 'dim')) + attrs_to_add <- names(var_info) + if (length(attrs_to_skip) > 0) { + attrs_to_add <- attrs_to_add[-attrs_to_skip] + } + for (attribute_name in attrs_to_add) { + if (is.numeric(var_info[[attribute_name]]) || + is.character(var_info[[attribute_name]])) { + ncatt_put(ncdf_object, defined_vars[[var_counter]]$name, attribute_name, var_info[[attribute_name]]) } - ncatt_put(ncdf_object, defined_vars[[var_counter]]$name, 'coordinates', var_info[['coordinates']]) } var_counter <- var_counter + 1 } diff --git a/R/CDORemap.R b/R/CDORemap.R index 24809c3f94f4ae37190068360af75d288f92188f..49bd9dbe1d280f64248442f9a66c37b3c741682f 100644 --- a/R/CDORemap.R +++ b/R/CDORemap.R @@ -42,11 +42,24 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, return_array <- FALSE if (length(dim(lons)) == 1) { array_dims <- c(length(lats), length(lons)) - names(array_dims) <- c('lat', 'lon') + new_lon_dim_name <- 'lon' + new_lat_dim_name <- 'lat' } else { array_dims <- dim(lons) - names(array_dims) <- c('j', 'i') + new_lon_dim_name <- 'i' + new_lat_dim_name <- 'j' } + if (!is.null(names(dim(lons)))) { + if (any(known_lon_names %in% names(dim(lons)))) { + new_lon_dim_name <- known_lon_names[which(known_lon_names %in% names(dim(lons)))[1]] + } + } + if (!is.null(names(dim(lats)))) { + if (any(known_lat_names %in% names(dim(lats)))) { + new_lat_dim_name <- known_lat_names[which(known_lat_names %in% names(dim(lats)))[1]] + } + } + names(array_dims) <- c(new_lat_dim_name, new_lon_dim_name) data_array <- array(NA, array_dims) } if (!(is.logical(data_array) || is.numeric(data_array)) || !is.array(data_array)) { @@ -86,7 +99,9 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, } } if (is.null(names(dim(lats)))) { - if (length(dim(lats)) > 1) { + if (length(dim(lats)) == 1) { + names(dim(lats)) <- lat_dim + } else { stop("Parameter 'lats' must be provided with dimension names.") } } else { @@ -203,7 +218,8 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, tmp_lon <- Subset(lons, lat_dim, min_pos[which(names(dim(lons)) == lat_dim)], drop = 'selected') } i <- 1:length(tmp_lon) - lon_model <- lm(tmp_lon ~ poly(i, 3)) + degree <- min(3, length(i) - 1) + lon_model <- lm(tmp_lon ~ poly(i, degree)) lon_extremes <- c(NA, NA) left_is_min <- FALSE right_is_max <- FALSE @@ -211,6 +227,7 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, left_is_min <- TRUE prev_lon <- predict(lon_model, data.frame(i = 0)) first_lon_cell_width <- (tmp_lon[1] - prev_lon) + # The signif is needed because cdo sellonlatbox crashes with too many digits lon_extremes[1] <- tmp_lon[1] - first_lon_cell_width / 2 } else { lon_extremes[1] <- min(tmp_lon) @@ -261,7 +278,8 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, tmp_lat <- Subset(lats, lon_dim, min_pos[which(names(dim(lats)) == lon_dim)], drop = 'selected') } i <- 1:length(tmp_lat) - lat_model <- lm(tmp_lat ~ poly(i, 3)) + degree <- min(3, length(i) - 1) + lat_model <- lm(tmp_lat ~ poly(i, degree)) lat_extremes <- c(NA, NA) if (which.min(tmp_lat) == 1) { prev_lat <- predict(lat_model, data.frame(i = 0)) @@ -474,24 +492,28 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, } names(attr(lons, 'variables')) <- lon_var_name names(attr(lats, 'variables')) <- lat_var_name + lons_lats_taken <- FALSE for (i in 1:total_slices) { tmp_file <- tempfile('R_CDORemap_', write_dir, fileext = '.nc') tmp_file2 <- tempfile('R_CDORemap_', write_dir, fileext = '.nc') if (!is.null(dims_to_iterate)) { slice_indices <- which(slices_to_iterate == i, arr.ind = TRUE) subset <- Subset(data_array, dims_to_iterate, as.list(slice_indices), drop = 'selected') +# dims_before_crop <- dim(subset) # Make sure subset goes along with metadata ArrayToNetCDF(setNames(list(subset, lons, lats), c('var', lon_var_name, lat_var_name)), tmp_file) } else { +# dims_before_crop <- dim(data_array) ArrayToNetCDF(setNames(list(data_array, lons, lats), c('var', lon_var_name, lat_var_name)), tmp_file) } sellonlatbox <- '' if (crop) { - sellonlatbox <- paste0('sellonlatbox,', lon_extremes[1], ',', lon_extremes[2], - ',', lat_extremes[1], ',', lat_extremes[2], ' -') + sellonlatbox <- paste0('sellonlatbox,', format(lon_extremes[1], scientific = FALSE), + ',', format(lon_extremes[2], scientific = FALSE), + ',', format(lat_extremes[1], scientific = FALSE), + ',', format(lat_extremes[2], scientific = FALSE), ' -') } err <- try({ -## TODO: Here add sellonlatbox. Also check constantin's issue, may contain hint. Also search if possible to crop without system(paste0("cdo -s ", sellonlatbox, "remap", method, ",", grid, " ", tmp_file, " ", tmp_file2)) }) file.remove(tmp_file) @@ -499,30 +521,45 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, stop("CDO remap failed.") } ncdf_remapped <- nc_open(tmp_file2) - found_dim_names <- sapply(ncdf_remapped$var$var$dim, '[[', 'name') - found_lon_dim <- found_dim_names[which(found_dim_names %in% .KnownLonNames())[1]] - found_lat_dim <- found_dim_names[which(found_dim_names %in% .KnownLatNames())[1]] - found_lon_dim_size <- length(ncdf_remapped$dim[[found_lon_dim]]$vals) - found_lat_dim_size <- length(ncdf_remapped$dim[[found_lat_dim]]$vals) - found_lons <- ncvar_get(ncdf_remapped, 'lon', collapse_degen = FALSE) - found_lats <- ncvar_get(ncdf_remapped, 'lat', collapse_degen = FALSE) - if (length(dim(found_lons)) > 1) { - if (found_lon_dim < found_lat_dim) { - names(dim(found_lons)) <- c(found_lon_dim, found_lat_dim) + if (!lons_lats_taken) { + found_dim_names <- sapply(ncdf_remapped$var$var$dim, '[[', 'name') + found_lon_dim <- found_dim_names[which(found_dim_names %in% .KnownLonNames())[1]] + found_lat_dim <- found_dim_names[which(found_dim_names %in% .KnownLatNames())[1]] + found_lon_dim_size <- length(ncdf_remapped$dim[[found_lon_dim]]$vals) + found_lat_dim_size <- length(ncdf_remapped$dim[[found_lat_dim]]$vals) + found_lons <- ncvar_get(ncdf_remapped, 'lon', collapse_degen = FALSE) + found_lats <- ncvar_get(ncdf_remapped, 'lat', collapse_degen = FALSE) + if (length(dim(lons)) == length(dim(found_lons))) { + new_lon_name <- lon_dim } else { - names(dim(found_lons)) <- c(found_lat_dim, found_lon_dim) + new_lon_name <- found_lon_dim } - } else { - names(dim(found_lons)) <- found_lon_dim - } - if (length(dim(found_lats)) > 1) { - if (found_lon_dim < found_lat_dim) { - names(dim(found_lats)) <- c(found_lon_dim, found_lat_dim) + if (length(dim(lats)) == length(dim(found_lats))) { + new_lat_name <- lat_dim } else { - names(dim(found_lats)) <- c(found_lat_dim, found_lon_dim) + new_lat_name <- found_lat_dim } - } else { - names(dim(found_lats)) <- found_lat_dim + if (length(dim(found_lons)) > 1) { + if (which(sapply(ncdf_remapped$var$lon$dim, '[[', 'name') == found_lon_dim) < + which(sapply(ncdf_remapped$var$lon$dim, '[[', 'name') == found_lat_dim)) { + names(dim(found_lons)) <- c(new_lon_name, new_lat_name) + } else { + names(dim(found_lons)) <- c(new_lat_name, new_lon_name) + } + } else { + names(dim(found_lons)) <- new_lon_name + } + if (length(dim(found_lats)) > 1) { + if (which(sapply(ncdf_remapped$var$lat$dim, '[[', 'name') == found_lon_dim) < + which(sapply(ncdf_remapped$var$lat$dim, '[[', 'name') == found_lat_dim)) { + names(dim(found_lats)) <- c(new_lon_name, new_lat_name) + } else { + names(dim(found_lats)) <- c(new_lat_name, new_lon_name) + } + } else { + names(dim(found_lats)) <- new_lat_name + } + lons_lats_taken <- TRUE } if (!is.null(dims_to_iterate)) { if (is.null(result_array)) { @@ -542,7 +579,7 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, new_dims <- dim(data_array) new_dims[c(lon_dim, lat_dim)] <- c(found_lon_dim_size, found_lat_dim_size) result_array <- ncvar_get(ncdf_remapped, 'var', collapse_degen = FALSE) - names(dim(result_array)) <- names(new_dims) + dim(result_array) <- new_dims } nc_close(ncdf_remapped) file.remove(tmp_file2) @@ -559,8 +596,6 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, } attribute_backup[['dim']][which(names(dim(result_array)) == lon_dim)] <- dim(result_array)[lon_dim] attribute_backup[['dim']][which(names(dim(result_array)) == lat_dim)] <- dim(result_array)[lat_dim] - new_lon_name <- names(dim(found_lons))[which(names(dim(found_lons)) %in% .KnownLonNames())] - new_lat_name <- names(dim(found_lats))[which(names(dim(found_lats)) %in% .KnownLatNames())] names(attribute_backup[['dim']])[which(names(dim(result_array)) == lon_dim)] <- new_lon_name names(attribute_backup[['dim']])[which(names(dim(result_array)) == lat_dim)] <- new_lat_name if (!is.null(attribute_backup[['variables']]) && (length(attribute_backup[['variables']]) > 0)) { diff --git a/R/Corr.R b/R/Corr.R index 10858c29494a38681c96381946144e228fd46101..4654db2ced466e1ff53be85824cf9f0c4b5725bc 100644 --- a/R/Corr.R +++ b/R/Corr.R @@ -149,8 +149,8 @@ Corr <- function(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL, } } if (pval && (method == "pearson")) { - t <- CORR * sqrt((eno - 2) / (1 - (CORR ^ 2))) - p_val <- 1 - pt(t, eno - 2) + t <-sqrt(CORR * CORR * (eno - 2) / (1 - (CORR ^ 2))) + p_val <- pt(t, eno - 2, lower.tail = FALSE) } if (conf && (method == "pearson")) { conf_low <- (1 - siglev) / 2 diff --git a/R/Load.R b/R/Load.R index 96b696a8d28b1830922e41c818eb8624a671aea5..efcc74c13ea9d37dbd282da71901f54ec2c84946 100644 --- a/R/Load.R +++ b/R/Load.R @@ -268,10 +268,14 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, # grid if (!is.null(grid)) { if (is.character(grid)) { - supported_grids <- list('r[0-9]{1,}x[0-9]{1,}', 't[0-9]{1,}grid') - grid_matches <- unlist(lapply(lapply(supported_grids, regexpr, grid), .IsFullMatch, grid)) - if (sum(grid_matches) < 1) { - stop("The specified grid in the parameter 'grid' is incorrect. Must be one of rx or tgrid.") + if (grid == 'none') { + grid <- NULL + } else { + supported_grids <- list('r[0-9]{1,}x[0-9]{1,}', 't[0-9]{1,}grid') + grid_matches <- unlist(lapply(lapply(supported_grids, regexpr, grid), .IsFullMatch, grid)) + if (sum(grid_matches) < 1) { + stop("The specified grid in the parameter 'grid' is incorrect. Must be one of rx or tgrid.") + } } } else { stop("Error: parameter 'grid' should be a character string, if specified.") diff --git a/R/PlotEquiMap.R b/R/PlotEquiMap.R index b8b7bd60fc89d6bd8d341f2162f35f280c6729f1..51599dc70f90f5440f61f0fab8bae9a5e3cf5663 100644 --- a/R/PlotEquiMap.R +++ b/R/PlotEquiMap.R @@ -83,16 +83,28 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL, dims <- dim(var) # Transpose the input matrices because the base plot functions work directly # with dimensions c(lon, lat). + transpose <- FALSE + if (!is.null(names(dims))) { + if (any(names(dims) %in% .KnownLonNames()) && + any(names(dims) %in% .KnownLatNames())) { + if (which(names(dims) %in% .KnownLonNames()) != 1) { + transpose <- TRUE + } + } + } if (dims[1] != length(lon) || dims[2] != length(lat)) { if (dims[1] == length(lat) && dims[2] == length(lon)) { - var <- t(var) - if (!is.null(varu)) varu <- t(varu) - if (!is.null(varv)) varv <- t(varv) - if (!is.null(contours)) contours <- t(contours) - if (!is.null(dots)) dots <- aperm(dots, c(1, 3, 2)) - dims <- dim(var) + transpose <- TRUE } } + if (transpose) { + var <- t(var) + if (!is.null(varu)) varu <- t(varu) + if (!is.null(varv)) varv <- t(varv) + if (!is.null(contours)) contours <- t(contours) + if (!is.null(dots)) dots <- aperm(dots, c(1, 3, 2)) + dims <- dim(var) + } # Check lon if (length(lon) != dims[1]) { diff --git a/R/PlotLayout.R b/R/PlotLayout.R index 36719d6d8505c16e359f22705dea0cbc6ebd0d9e..7d5eee1efc2d184e50bb535ba4088fcf9b736c65 100644 --- a/R/PlotLayout.R +++ b/R/PlotLayout.R @@ -127,21 +127,19 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL, } # Check the rest of parameters (unless the user simply wants to build an empty layout) - if (!(drawleg == FALSE)) { - var_limits <- NULL - if (!all(sapply(var, is_single_na))) { - var_limits <- c(min(unlist(var), na.rm = TRUE), max(unlist(var), na.rm = TRUE)) - if ((any(is.infinite(var_limits)) || var_limits[1] == var_limits[2])) { - stop("Arrays in parameter 'var' must contain at least 2 different values.") - } + var_limits <- NULL + if (!all(sapply(var, is_single_na))) { + var_limits <- c(min(unlist(var), na.rm = TRUE), max(unlist(var), na.rm = TRUE)) + if ((any(is.infinite(var_limits)) || var_limits[1] == var_limits[2])) { + stop("Arrays in parameter 'var' must contain at least 2 different values.") } - colorbar <- ColorBar(brks, cols, FALSE, subsampleg, bar_limits, - var_limits, triangle_ends, col_inf, col_sup, color_fun, - plot = FALSE, draw_bar_ticks, - draw_separators, triangle_ends_scale, bar_extra_labels, - units, units_scale, bar_label_scale, bar_tick_scale, - bar_extra_margin, bar_label_digits) } + colorbar <- ColorBar(brks, cols, FALSE, subsampleg, bar_limits, + var_limits, triangle_ends, col_inf, col_sup, color_fun, + plot = FALSE, draw_bar_ticks, + draw_separators, triangle_ends_scale, bar_extra_labels, + units, units_scale, bar_label_scale, bar_tick_scale, + bar_extra_margin, bar_label_digits) # Check bar_scale if (!is.numeric(bar_scale)) { @@ -227,8 +225,8 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL, dim_ids <- plot_dims[[plot_array_i]] if (is.character(dim_ids)) { dimnames <- NULL - if (!is.null(names(plot_array))) { - dimnames <- names(plot_array) + if (!is.null(names(dim(plot_array)))) { + dimnames <- names(dim(plot_array)) } else if (!is.null(attr(plot_array, 'dimensions'))) { dimnames <- attr(plot_array, 'dimensions') } @@ -237,7 +235,7 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL, stop("All arrays provided in parameter 'var' must have all the dimensions in 'plot_dims'.") } dim_ids <- sapply(dim_ids, function(x) which(dimnames == x)[1]) - var[[plot_array_i]] <- aperm(var[[plot_array_i]], c((1:length(dim(plot_array)))[-dim_ids], dim_ids)) + var[[plot_array_i]] <- .aperm2(var[[plot_array_i]], c((1:length(dim(plot_array)))[-dim_ids], dim_ids)) } else { .warning(paste0("Assuming the ", plot_array_i, "th array provided in 'var' has 'plot_dims' as last dimensions (right-most).")) dims <- tail(c(rep(1, length(dim_ids)), dim(plot_array)), length(dim_ids)) @@ -434,8 +432,13 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL, } plot_number <<- plot_number + 1 } else { + if (is.character(plot_dims[[array_number]])) { + plot_dim_indices <- which(names(dim(x)) %in% plot_dims[[array_number]]) + } else { + plot_dim_indices <- plot_dims[[array_number]] + } # For each of the arrays provided in that array - apply(x, (1:length(dim(x)))[1:(length(dim(x)) - length(plot_dims[[array_number]]))], + apply(x, (1:length(dim(x)))[-plot_dim_indices], function(y) { # Do the plot fun_args <- c(list(y, toptitle = titles[plot_number]), list(...), diff --git a/R/Utils.R b/R/Utils.R index 806ad92d448b709655a06d34aae3d28f3f7f4310..c748a50876375c615c7b8d7c8a387475bc727614 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -220,6 +220,9 @@ # Here we read the grid type and its number of longitudes and latitudes file_info <- system(paste('cdo -s griddes', filein, '2> /dev/null'), intern = TRUE) grids_positions <- grep('# gridID', file_info) + if (length(grids_positions) < 1) { + stop("The grid should be defined in the files.") + } grids_first_lines <- grids_positions + 2 grids_last_lines <- c((grids_positions - 2)[-1], length(file_info)) grids_info <- as.list(1:length(grids_positions)) @@ -230,21 +233,23 @@ grids_types <- unlist(lapply(grids_info, function (x) x[grep('gridtype', x) + 1])) grids_matches <- unlist(lapply(grids_info, function (x) { nlons <- if (length(grep('xsize', x)) > 0) { - as.integer(x[grep('xsize', x) + 1]) + as.numeric(x[grep('xsize', x) + 1]) } else { NA } nlats <- if (length(grep('ysize', x)) > 0) { - as.integer(x[grep('ysize', x) + 1]) + as.numeric(x[grep('ysize', x) + 1]) } else { NA } - if (identical(nlons, length(lon)) && - identical(nlats, length(lat))) { - TRUE - } else { - FALSE + result <- FALSE + if (!any(is.na(c(nlons, nlats)))) { + if ((nlons == length(lon)) && + (nlats == length(lat))) { + result <- TRUE + } } + result })) grids_matches <- grids_matches[which(grids_types %in% c('gaussian', 'lonlat'))] grids_info <- grids_info[which(grids_types %in% c('gaussian', 'lonlat'))] @@ -261,8 +266,10 @@ } else { stop("Error: Load() can't disambiguate: More than one lonlat/gaussian grids with the same size as the requested variable defined in ", filename) } - } else { + } else if (sum(grids_matches) == 1) { grid_type <- grids_types[which(grids_matches)] + } else { + stop("Unexpected error.") } grid_lons <- length(lon) grid_lats <- length(lat) @@ -279,14 +286,14 @@ # Now we calculate the common grid type and its lons and lats if (length(grep('^t\\d{1,+}grid$', work_piece[['grid']])) > 0) { common_grid_type <- 'gaussian' - common_grid_res <- as.integer(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][2]) + common_grid_res <- as.numeric(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][2]) nlonlat <- .t2nlatlon(common_grid_res) common_grid_lats <- nlonlat[1] common_grid_lons <- nlonlat[2] } else if (length(grep('^r\\d{1,+}x\\d{1,+}$', work_piece[['grid']])) > 0) { common_grid_type <- 'lonlat' - common_grid_lons <- as.integer(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][2]) - common_grid_lats <- as.integer(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][3]) + common_grid_lons <- as.numeric(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][2]) + common_grid_lats <- as.numeric(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][3]) } else { stop("Error: Only supported grid types in parameter 'grid' are tgrid and rx") } @@ -341,12 +348,14 @@ if (!is.null(work_piece[['progress_amount']])) { cat("\n") } - cat(paste0("! Warning: the dataset with index ", tail(work_piece[['indices']], 1), - " in '", work_piece[['dataset_type']], "' is originally on ", - "a grid coarser than the common grid and it has been ", - "extrapolated. Check the results carefully. It is ", - "recommended to specify as common grid the coarsest grid ", - "among all requested datasets via the parameter 'grid'.\n")) + if (!explore_dims) { + cat(paste0("! Warning: the dataset with index ", tail(work_piece[['indices']], 1), + " in '", work_piece[['dataset_type']], "' is originally on ", + "a grid coarser than the common grid and it has been ", + "extrapolated. Check the results carefully. It is ", + "recommended to specify as common grid the coarsest grid ", + "among all requested datasets via the parameter 'grid'.\n")) + } } # Now calculate if the user requests for a lonlat subset or for the # entire field @@ -610,15 +619,15 @@ ' 2>/dev/null'), intern = TRUE), split = ' ') years <- strsplit(system(paste('cdo showyear ', filein, ' 2>/dev/null'), intern = TRUE), split = ' ') - mons <- as.integer(mons[[1]][which(mons[[1]] != "")]) - years <- as.integer(years[[1]][which(years[[1]] != "")]) + mons <- as.numeric(mons[[1]][which(mons[[1]] != "")]) + years <- as.numeric(years[[1]][which(years[[1]] != "")]) time_indices <- ts(time_indices, start = c(years[1], mons[1]), end = c(years[length(years)], mons[length(mons)]), frequency = 12) ltimes_list <- list() for (sdate in work_piece[['startdates']]) { - selected_time_indices <- window(time_indices, start = c(as.integer( - substr(sdate, 1, 4)), as.integer(substr(sdate, 5, 6))), + selected_time_indices <- window(time_indices, start = c(as.numeric( + substr(sdate, 1, 4)), as.numeric(substr(sdate, 5, 6))), end = c(3000, 12), frequency = 12, extend = TRUE) selected_time_indices <- selected_time_indices[work_piece[['leadtimes']]] ltimes_list <- c(ltimes_list, list(selected_time_indices)) @@ -1156,7 +1165,7 @@ if (give_warning) { .warning(paste0("Too complex path pattern specified for ", dataset_name, - ". Double check carefully the 'source_files' fetched for this dataset or specify a simpler path pattern.")) + ". Double check carefully the '$Files' fetched for this dataset or specify a simpler path pattern.")) } if (permissive) { @@ -1166,6 +1175,60 @@ } } +.FindTagValue <- function(path_with_globs_and_tag, actual_path, tag) { + tag <- paste0('\\$', tag, '\\$') + path_with_globs_and_tag <- paste0('^', path_with_globs_and_tag, '$') + parts <- strsplit(path_with_globs_and_tag, '*', fixed = TRUE)[[1]] + parts <- as.list(parts[grep(tag, parts)]) + longest_couples <- c() + pos_longest_couples <- c() + found_value <- NULL + for (i in 1:length(parts)) { + parts[[i]] <- strsplit(parts[[i]], tag)[[1]] + if (length(parts[[i]]) == 1) { + parts[[i]] <- c(parts[[i]], '') + } + len_parts <- sapply(parts[[i]], nchar) + len_couples <- len_parts[-length(len_parts)] + len_parts[2:length(len_parts)] + pos_longest_couples <- c(pos_longest_couples, which.max(len_couples)) + longest_couples <- c(longest_couples, max(len_couples)) + } + chosen_part <- which.max(longest_couples) + parts[[chosen_part]] <- parts[[chosen_part]][pos_longest_couples[chosen_part]:(pos_longest_couples[chosen_part] + 1)] + if (nchar(parts[[chosen_part]][1]) >= nchar(parts[[chosen_part]][2])) { + if (nchar(parts[[chosen_part]][1]) > 0) { + matches <- gregexpr(parts[[chosen_part]][1], actual_path)[[1]] + if (length(matches) == 1) { + match_left <- matches + actual_path <- substr(actual_path, match_left + attr(match_left, 'match.length'), nchar(actual_path)) + } + } + if (nchar(parts[[chosen_part]][2]) > 0) { + matches <- gregexpr(parts[[chosen_part]][2], actual_path)[[1]] + if (length(matches) == 1) { + match_right <- matches + found_value <- substr(actual_path, 0, match_right - 1) + } + } + } else { + if (nchar(parts[[chosen_part]][2]) > 0) { + matches <- gregexpr(parts[[chosen_part]][2], actual_path)[[1]] + if (length(matches) == 1) { + match_right <- matches + actual_path <- substr(actual_path, 0, match_right - 1) + } + } + if (nchar(parts[[chosen_part]][1]) > 0) { + matches <- gregexpr(parts[[chosen_part]][1], actual_path)[[1]] + if (length(matches) == 1) { + match_left <- matches + found_value <- substr(actual_path, match_left + attr(match_left, 'match.length'), nchar(actual_path)) + } + } + } + found_value +} + .FilterUserGraphicArgs <- function(excludedArgs, ...) { # This function filter the extra graphical parameters passed by the user in # a plot function, excluding the ones that the plot function uses by default. @@ -1438,3 +1501,90 @@ plot(0, type = 'n', axes = FALSE, ann = FALSE) par(mfg = next_attempt) } + +# Function to permute arrays of non-atomic elements (e.g. POSIXct) +.aperm2 <- function(x, new_order) { + y <- array(1:length(x), dim = dim(x)) + y <- aperm(y, new_order) + old_dims <- dim(x) + x <- x[as.vector(y)] + dim(x) <- old_dims[new_order] + x +} + +# This function is a helper for the function .MergeArrays. +# It expects as inputs two named numeric vectors, and it extends them +# with dimensions of length 1 until an ordered common dimension +# format is reached. +.MergeArrayDims <- function(dims1, dims2) { + new_dims1 <- c() + new_dims2 <- c() + while (length(dims1) > 0) { + if (names(dims1)[1] %in% names(dims2)) { + pos <- which(names(dims2) == names(dims1)[1]) + dims_to_add <- rep(1, pos - 1) + if (length(dims_to_add) > 0) { + names(dims_to_add) <- names(dims2[1:(pos - 1)]) + } + new_dims1 <- c(new_dims1, dims_to_add, dims1[1]) + new_dims2 <- c(new_dims2, dims2[1:pos]) + dims1 <- dims1[-1] + dims2 <- dims2[-c(1:pos)] + } else { + new_dims1 <- c(new_dims1, dims1[1]) + new_dims2 <- c(new_dims2, 1) + names(new_dims2)[length(new_dims2)] <- names(dims1)[1] + dims1 <- dims1[-1] + } + } + if (length(dims2) > 0) { + dims_to_add <- rep(1, length(dims2)) + names(dims_to_add) <- names(dims2) + new_dims1 <- c(new_dims1, dims_to_add) + new_dims2 <- c(new_dims2, dims2) + } + list(new_dims1, new_dims2) +} + +# This function takes two named arrays and merges them, filling with +# NA where needed. +# dim(array1) +# 'b' 'c' 'e' 'f' +# 1 3 7 9 +# dim(array2) +# 'a' 'b' 'd' 'f' 'g' +# 2 3 5 9 11 +# dim(.MergeArrays(array1, array2, 'b')) +# 'a' 'b' 'c' 'e' 'd' 'f' 'g' +# 2 4 3 7 5 9 11 +.MergeArrays <- function(array1, array2, along) { + if (!(identical(names(dim(array1)), names(dim(array2))) && + identical(dim(array1)[-which(names(dim(array1)) == along)], + dim(array2)[-which(names(dim(array2)) == along)]))) { + new_dims <- .MergeArrayDims(dim(array1), dim(array2)) + dim(array1) <- new_dims[[1]] + dim(array2) <- new_dims[[2]] + for (j in 1:length(dim(array1))) { + if (names(dim(array1))[j] != along) { + if (dim(array1)[j] != dim(array2)[j]) { + if (which.max(c(dim(array1)[j], dim(array2)[j])) == 1) { + na_array_dims <- dim(array2) + na_array_dims[j] <- dim(array1)[j] - dim(array2)[j] + na_array <- array(dim = na_array_dims) + array2 <- abind(array2, na_array, along = j) + names(dim(array2)) <- names(na_array_dims) + } else { + na_array_dims <- dim(array1) + na_array_dims[j] <- dim(array2)[j] - dim(array1)[j] + na_array <- array(dim = na_array_dims) + array1 <- abind(array1, na_array, along = j) + names(dim(array1)) <- names(na_array_dims) + } + } + } + } + } + array1 <- abind(array1, array2, along = which(names(dim(array1)) == along)) + names(dim(array1)) <- names(dim(array2)) + array1 +} diff --git a/man/ArrayToNetCDF.Rd b/man/ArrayToNetCDF.Rd index b857eaf2396b05aa9f0dcd26eea5dbc249bf520e..fcfb77005610e6eb13fa044e3b0c49e2a158c876 100644 --- a/man/ArrayToNetCDF.Rd +++ b/man/ArrayToNetCDF.Rd @@ -13,8 +13,8 @@ ArrayToNetCDF(arrays = list(temperature = array(1:9, c(3, 3))), \item{Via the dimension names of each provided array:}{The dimension names of each of the provided arrays will be interpreted as names for the dimensions of the NetCDF files. Read further for special dimension names that will trigger special behaviours, such as 'time' and 'var'.\cr E.g:\cr \code{ -temperature <- array(rnorm(10 * 50 * 100), dim = c(10, 50, 100)) -names(dim(temperature)) <- c('time', 'latitude', 'longitude') +temperature <- array(rnorm(100 * 50 * 10), dim = c(100, 50, 10)) +names(dim(temperature)) <- c('longitude', 'latitude', 'time') ArrayToNetCDF(list(temperature = temperature), file_path = 'example.nc') } } diff --git a/man/CDORemap.Rd b/man/CDORemap.Rd index 7960d7bbf5395f06d921f9465f55c9e7aeb59365..fa539dc6569f2b1540ba765b174b7e0edf854ce0 100644 --- a/man/CDORemap.Rd +++ b/man/CDORemap.Rd @@ -14,7 +14,7 @@ CDORemap(data_array = NULL, lons, lats, grid, method, \item{data_array}{Multidimensional numeric array to be interpolated. If provided, it must have at least a longitude and a latitude dimensions, identified by the array dimension names. The names for these dimensions must be one of the recognized by s2dverification (can be checked with \code{s2dverification:::.KnownLonNames()} and \code{s2dverification:::.KnownLatNames()}).} \item{lons}{Numeric vector or array of longitudes of the centers of the grid cells. Its size must match the size of the longitude/latitude dimensions of the input array.} \item{lats}{Numeric vector or array of latitudes of the centers of the grid cells. Its size must match the size of the longitude/latitude dimensions of the input array.} - \item{grid}{Character string specifying either a name of a grid (recognized by CDO; e.g.: 'r256x128', 't106grid') or a path to another NetCDF file which to read the grid from (a single grid must be defined in such file).} + \item{grid}{Character string specifying either a name of a target grid (recognized by CDO; e.g.: 'r256x128', 't106grid') or a path to another NetCDF file which to read the target grid from (a single grid must be defined in such file).} \item{method}{Character string specifying an interpolation method (recognized by CDO; e.g.: 'con', 'bil', 'bic', 'dis'). The following long names are also supported: 'conservative', 'bilinear', 'bicubic' and 'distance-weighted'.} \item{avoid_writes}{The step of permutation is needed when the input array has more than 3 dimensions and none of the longitude or latitude dimensions in the right-most position (CDO would not accept it without permuting previously). This step, executed by default when needed, can be avoided for the price of writing more intermediate files (whis usually is unconvenient) by setting the parameter \code{avoid_writes = TRUE}.} \item{crop}{Whether to crop the data after interpolation with 'cdo sellonlatbox' (TRUE) or to extend interpolated data to the whole world as CDO does by default (FALSE). If \code{crop = TRUE} then the longitude and latitude borders which to crop at are taken as the limits of the cells at the borders ('lons' and 'lats' are perceived as cell centers), i.e. the resulting array will contain data that covers the same area as the input array. This is equivalent to specifying \code{crop = 'preserve'}, i.e. preserving area. If \code{crop = 'tight'} then the borders which to crop at are taken as the minimum and maximum cell centers in 'lons' and 'lats', i.e. the area covered by the resulting array may be smaller if interpolating from a coarse grid to a fine grid. The parameter 'crop' also accepts a numeric vector of custom borders which to crop at: c(western border, eastern border, southern border, northern border). } diff --git a/man/s2dverification.Rd b/man/s2dverification.Rd index af6beaaea1857176dad23d6f3b4fd700cb55baea..ee3fca8036ed945a86a8253a531e73bf8124c685 100644 --- a/man/s2dverification.Rd +++ b/man/s2dverification.Rd @@ -10,8 +10,8 @@ Set of tools to verify forecasts through the computation of typical prediction s \tabular{ll}{ Package: \tab s2dverification\cr Type: \tab Package\cr -Version: \tab 2.8.1\cr -Date: \tab 2017-02-15\cr +Version: \tab 2.8.2\cr +Date: \tab 2017-05-15\cr License: \tab LGPLv3\cr } Check an overview of the package functionalities and its modules at \url{https://earth.bsc.es/gitlab/es/s2dverification/wikis/home}. diff --git a/s2dverification-manual.pdf b/s2dverification-manual.pdf index 01c10e3721942fc4dfc4179606bda389e0b0a1ba..8c8bdc7fc49c2a8c052aefcaa8a4d462386c393c 100644 Binary files a/s2dverification-manual.pdf and b/s2dverification-manual.pdf differ