diff --git a/R/CDORemap.R b/R/CDORemap.R index efd9bf216c27e02795f6864938360c5f23e7d4f4..88c2cd42187cf84d6b1d99bd58ffd6b10554ddc3 100644 --- a/R/CDORemap.R +++ b/R/CDORemap.R @@ -398,6 +398,7 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method, cdo_version <- as.numeric_version( strsplit(suppressWarnings(system2("cdo", args = '-V', stderr = TRUE))[[1]], ' ')[[1]][5] ) + warning("CDORemap: Using CDO version ", cdo_version, ".") if ((cdo_version >= as.numeric_version('1.7.0')) && (method == 'con')) { method <- 'ycon' } diff --git a/R/Utils.R b/R/Utils.R index 31f9ae071fbb78f4d398cf447cdf907f9aaaa645..d8ffa744c673527c7085450e96528e1ef7a2c290 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -733,16 +733,16 @@ if (fnc$var[[namevar]]$hasScaleFact || fnc$var[[namevar]]$hasAddOffset) { tmp <- (tmp - add_offset) / scale_factor } - nc_close(fnc) - fnc <- nc_create(filein2, list(ncdf_var)) - ncvar_put(fnc, ncdf_var, tmp) + #nc_close(fnc) + fnc2 <- nc_create(filein2, list(ncdf_var)) + ncvar_put(fnc2, ncdf_var, tmp) if (add_offset != 0) { - ncatt_put(fnc, ncdf_var, 'add_offset', add_offset) + ncatt_put(fnc2, ncdf_var, 'add_offset', add_offset) } if (scale_factor != 1) { - ncatt_put(fnc, ncdf_var, 'scale_factor', scale_factor) + ncatt_put(fnc2, ncdf_var, 'scale_factor', scale_factor) } - nc_close(fnc) + nc_close(fnc2) system(paste0("cdo -s -sellonlatbox,", if (lonmin > lonmax) { "0,360," } else { @@ -751,27 +751,27 @@ " -remap", work_piece[['remap']], ",", common_grid_name, " ", filein2, " ", filein, " 2>/dev/null", sep = "")) file.remove(filein2) - fnc <- nc_open(filein) - lon <- ncvar_get(fnc, 'lon') - lat <- ncvar_get(fnc, 'lat') + fnc2 <- nc_open(filein) + sub_lon <- ncvar_get(fnc2, 'lon') + sub_lat <- ncvar_get(fnc2, 'lat') ## We read the longitudes and latitudes from the file. ## In principle cdo should put in order the longitudes ## and slice them properly unless data is across greenwich - lon[which(lon < 0)] <- lon[which(lon < 0)] + 360 - lon_indices <- 1:length(lon) + sub_lon[which(sub_lon < 0)] <- sub_lon[which(sub_lon < 0)] + 360 + sub_lon_indices <- 1:length(sub_lon) if (lonmax < lonmin) { - lon_indices <- lon_indices[which((lon <= lonmax) | (lon >= lonmin))] + sub_lon_indices <- sub_lon_indices[which((sub_lon <= lonmax) | (sub_lon >= lonmin))] } - lat_indices <- 1:length(lat) + sub_lat_indices <- 1:length(sub_lat) ## In principle cdo should put in order the latitudes - if (lat[1] < lat[length(lat)]) { - lat_indices <- length(lat):1 + if (sub_lat[1] < sub_lat[length(sub_lat)]) { + sub_lat_indices <- length(sub_lat):1 } - final_dims[c(1, 2)] <- c(length(lon_indices), length(lat_indices)) - subset_indices[[dim_matches[1]]] <- lon_indices - subset_indices[[dim_matches[2]]] <- lat_indices + final_dims[c(1, 2)] <- c(length(sub_lon_indices), length(sub_lat_indices)) + subset_indices[[dim_matches[1]]] <- sub_lon_indices + subset_indices[[dim_matches[2]]] <- sub_lat_indices - tmp <- take(ncvar_get(fnc, namevar, collapse_degen = FALSE), + tmp <- take(ncvar_get(fnc2, namevar, collapse_degen = FALSE), 1:length(subset_indices), subset_indices) if (!is.null(mask)) { @@ -802,13 +802,13 @@ mask_lon_indices <- which((mask_lons >= lonmin) | (mask_lons <= lonmax)) } mask_lat_indices <- which((mask_lats >= latmin) & (mask_lats <= latmax)) - if (lat[1] < lat[length(lat)]) { + if (sub_lat[1] < sub_lat[length(sub_lat)]) { mask_lat_indices <- mask_lat_indices[length(mask_lat_indices):1] } mask <- mask[mask_lon_indices, mask_lat_indices] } - lon <- lon[lon_indices] - lat <- lat[lat_indices] + sub_lon <- sub_lon[sub_lon_indices] + sub_lat <- sub_lat[sub_lat_indices] ### nc_close(fnc_mask) ### system(paste0("cdo -s -sellonlatbox,", if (lonmin > lonmax) { ### "0,360," @@ -834,6 +834,13 @@ dim(tmp) <- final_dims # If we are exploring the file we don't need to process and arrange # the retrieved data. We only need to keep the dimension sizes. + if (is_2d_var && lonlat_subsetting_requested && remap_needed) { + final_lons <- sub_lon + final_lats <- sub_lat + } else { + final_lons <- lon + final_lats <- lat + } if (explore_dims) { if (work_piece[['is_file_per_member']]) { ## TODO: When the exp_full_path contains asterisks and is file_per_member @@ -846,7 +853,7 @@ nmemb <- length(files) } } - dims <- list(member = nmemb, ftime = nltime, lon = lon, lat = lat) + dims <- list(member = nmemb, ftime = nltime, lon = final_lons, lat = final_lats) } else { # If we are not exploring, then we have to process the retrieved data if (is_2d_var) { @@ -863,13 +870,13 @@ } if (output == 'areave' || output == 'lon') { - weights <- InsertDim(cos(lat * pi / 180), 1, length(lon)) + weights <- InsertDim(cos(final_lats * pi / 180), 1, length(final_lons)) weights[which(is.na(x))] <- NA if (output == 'areave') { weights <- weights / mean(weights, na.rm = TRUE) mean(x * weights, na.rm = TRUE) } else { - weights <- weights / InsertDim(Mean1Dim(weights, 2, narm = TRUE), 2, length(lat)) + weights <- weights / InsertDim(Mean1Dim(weights, 2, narm = TRUE), 2, length(final_lats)) Mean1Dim(x * weights, 2, narm = TRUE) } } else if (output == 'lat') { @@ -1511,11 +1518,20 @@ # 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)] + attr_bk <- attributes(x) + if ('dim' %in% names(attr_bk)) { + attr_bk[['dim']] <- NULL + } + if (is.numeric(x)) { + x <- aperm(x, new_order) + } else { + y <- array(1:length(x), dim = dim(x)) + y <- aperm(y, new_order) + x <- x[as.vector(y)] + } dim(x) <- old_dims[new_order] + attributes(x) <- c(attributes(x), attr_bk) x }