From 86a2f54f69757e0b52d28ad5d174b15b6b1d1477 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 19 May 2021 11:34:48 +0200 Subject: [PATCH 1/2] Replace ArrayToNetCDF in CDORemap with easyNCDF::ArrayToNc --- DESCRIPTION | 3 ++- NAMESPACE | 1 + R/CDORemap.R | 9 +++++---- man/CDORemap.Rd | 2 +- 4 files changed, 9 insertions(+), 6 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index f8f5567..9a5a771 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -36,7 +36,8 @@ Imports: plyr, ncdf4, NbClust, - multiApply (>= 2.1.1) + multiApply (>= 2.1.1), + easyNCDF Suggests: easyVerification, testthat diff --git a/NAMESPACE b/NAMESPACE index 87937e2..6e7d545 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -69,6 +69,7 @@ importFrom(ClimProjDiags,Subset) importFrom(ClimProjDiags,WeightedMean) importFrom(abind,abind) importFrom(abind,adrop) +importFrom(easyNCDF,ArrayToNc) importFrom(grDevices,bmp) importFrom(grDevices,col2rgb) importFrom(grDevices,colorRampPalette) diff --git a/R/CDORemap.R b/R/CDORemap.R index ae3c988..0ae0be6 100644 --- a/R/CDORemap.R +++ b/R/CDORemap.R @@ -1,4 +1,4 @@ -#'Interpolates arrays with longitude and latitude dimensions using CDO +#'Interpolate arrays with longitude and latitude dimensions using CDO #' #'This function takes as inputs a multidimensional array (optional), a vector #'or matrix of longitudes, a vector or matrix of latitudes, a destination grid @@ -202,6 +202,7 @@ #'tas2 <- CDORemap(tas, lon, lat, 'external_file.nc', 'bil') #'} #'@import ncdf4 +#'@importFrom easyNCDF ArrayToNc #'@importFrom stats lm predict setNames #'@export CDORemap <- function(data_array = NULL, lons, lats, grid, method, @@ -684,7 +685,7 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, names(dim(data_array)) <- paste0('dim', 1:length(dim(data_array))) names(dim(data_array))[c(lon_pos, lat_pos)] <- c(lon_dim, lat_dim) if (!is.null(unlimited_dim)) { - # This will make ArrayToNetCDF create this dim as unlimited. + # This will make ArrayToNc create this dim as unlimited. names(dim(data_array))[unlimited_dim] <- 'time' } if (length(dim(lons)) == 1) { @@ -725,10 +726,10 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, 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) + easyNCDF::ArrayToNc(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) + easyNCDF::ArrayToNc(setNames(list(data_array, lons, lats), c('var', lon_var_name, lat_var_name)), tmp_file) } sellonlatbox <- '' if (crop) { diff --git a/man/CDORemap.Rd b/man/CDORemap.Rd index aabb434..c7a00a8 100644 --- a/man/CDORemap.Rd +++ b/man/CDORemap.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/CDORemap.R \name{CDORemap} \alias{CDORemap} -\title{Interpolates arrays with longitude and latitude dimensions using CDO} +\title{Interpolate arrays with longitude and latitude dimensions using CDO} \usage{ CDORemap( data_array = NULL, -- GitLab From 5d0412a8e6b0105d51aa20dac16b80994ffe3fc6 Mon Sep 17 00:00:00 2001 From: aho Date: Fri, 4 Jun 2021 14:29:56 +0200 Subject: [PATCH 2/2] Enable to regrid irregular curvilinear grid to regular grid --- R/CDORemap.R | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/R/CDORemap.R b/R/CDORemap.R index 0ae0be6..fc25b52 100644 --- a/R/CDORemap.R +++ b/R/CDORemap.R @@ -724,10 +724,45 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, 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') + # Fix issue 259, curvilinear grid, the order of the dimensions in slices and + # coordinates needs to match + if (is_irregular) { + pos_lon <- which(names(dim(subset)) == lon_dim) + pos_lat <- which(names(dim(subset)) == lat_dim) + pos_lon_dim_in_lons <- which(names(dim(lons)) == lon_dim) + pos_lat_dim_in_lons <- which(names(dim(lons)) == lat_dim) + if ((pos_lon > pos_lat && pos_lon_dim_in_lons < pos_lat_dim_in_lons) || + (pos_lon < pos_lat && pos_lon_dim_in_lons > pos_lat_dim_in_lons)) { + new_pos <- 1:length(dim(subset)) + new_pos[pos_lon] <- pos_lat + new_pos[pos_lat] <- pos_lon + subset <- .aperm2(subset, new_pos) + } + # The unlimited dimension should be placed in the last position + if ('time' %in% names(dim(subset)) && + which(names(dim(subset)) == 'time') != length(dim(subset))) { + new_pos <- 2:length(dim(subset)) + new_pos[length(dim(subset))] <- 1 + subset <- .aperm2(subset, new_pos) + } + } # dims_before_crop <- dim(subset) # Make sure subset goes along with metadata easyNCDF::ArrayToNc(setNames(list(subset, lons, lats), c('var', lon_var_name, lat_var_name)), tmp_file) } else { + if (is_irregular) { + pos_lon <- which(names(dim(data_array)) == lon_dim) + pos_lat <- which(names(dim(data_array)) == lat_dim) + pos_lon_dim_in_lons <- which(names(dim(lons)) == lon_dim) + pos_lat_dim_in_lons <- which(names(dim(lons)) == lat_dim) + if ((pos_lon > pos_lat && pos_lon_dim_in_lons < pos_lat_dim_in_lons) || + (pos_lon < pos_lat && pos_lon_dim_in_lons > pos_lat_dim_in_lons)) { + new_pos <- 1:length(dim(data_array)) + new_pos[pos_lon] <- pos_lat + new_pos[pos_lat] <- pos_lon + data_array <- .aperm2(data_array, new_pos) + } + } # dims_before_crop <- dim(data_array) easyNCDF::ArrayToNc(setNames(list(data_array, lons, lats), c('var', lon_var_name, lat_var_name)), tmp_file) } @@ -816,24 +851,63 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, if (return_array) { new_dims <- dim(data_array) new_dims[c(lon_dim, lat_dim)] <- c(found_lon_dim_size, found_lat_dim_size) + lon_pos <- which(names(new_dims) == lon_dim) + lat_pos <- which(names(new_dims) == lat_dim) + # Fix issue 259, expected order from CDO output is lon lat + # If is irregular, lat and lon position need to be checked: + if (is_irregular) { + if (lon_pos > lat_pos) { + new_pos <- 1:length(new_dims) + new_pos[lon_pos] <- lat_pos + new_pos[lat_pos] <- lon_pos + new_dims <- new_dims[new_pos] + } + } result_array <- array(dim = new_dims) store_indices <- as.list(rep(TRUE, length(dim(result_array)))) } } if (return_array) { store_indices[dims_to_iterate] <- as.list(slice_indices) + # If is irregular, the order of dimenesions in result_array and file may be different and need to be checked before reading the temporal file: + if (is_irregular) { + test_dims <- dim(ncvar_get(ncdf_remapped, 'var', + collapse_degen = FALSE)) + test_dims <- test_dims[which(test_dims > 1)] + pos_test_dims <- match(dim(result_array), test_dims) + if (is.unsorted(pos_test_dims, na.rm = TRUE)) { + # pos_new_dims is used later in the code. Don't overwrite + pos_new_dims <- 1:length(dim(result_array)) + pos_new_dims[which(!is.na(pos_test_dims))] <- + match(test_dims, dim(result_array)) + backup_result_array_dims <- dim(result_array) + dim(result_array) <- dim(result_array)[pos_new_dims] + } + } result_array <- do.call('[<-', c(list(x = result_array), store_indices, list(value = ncvar_get(ncdf_remapped, 'var', collapse_degen = FALSE)))) } } else { 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) dim(result_array) <- new_dims } nc_close(ncdf_remapped) file.remove(tmp_file2) } + # If is irregular, the order of dimension may need to be recovered after reading all the file: + if (is_irregular & (!is.null(dims_to_iterate))) { + if (exists('pos_new_dims')) { + pos_new_dims <- 1:length(dim(result_array)) + dims_to_change <- match(backup_result_array_dims, dim(result_array)) + pos_new_dims[which(dims_to_change != 1)] <- + dims_to_change[which(dims_to_change != 1)] + result_array <- .aperm2(result_array, pos_new_dims) + } + } + if (!is.null(permutation)) { dim_backup <- dim(result_array) result_array <- aperm(result_array, permutation_back) -- GitLab