diff --git a/NAMESPACE b/NAMESPACE index 696a3e89078ea52b1a70f32e6cb00738166e76d1..342847f35714e9473b71200e7c58f7f19515e5ac 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -17,6 +17,7 @@ export(CST_MultiMetric) export(CST_MultivarRMSE) export(CST_QuantileMapping) export(CST_RFSlope) +export(CST_RFTemp) export(CST_RFWeights) export(CST_RainFARM) export(CST_RegimesAssign) @@ -33,6 +34,7 @@ export(PlotMostLikelyQuantileMap) export(PlotPDFsOLE) export(PlotTriangles4Categories) export(RFSlope) +export(RFTemp) export(RainFARM) export(RegimesAssign) export(SplitDim) diff --git a/NEWS.md b/NEWS.md index 208fec3928a3834ae63edfc05f78da7354981a6f..19caa068dc1f95e7f3a710ae5d05e45292f78317 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,7 @@ + PlotTriangles4Categories new plotting function to convert any 3-d numerical array to a grid of coloured triangles. + CST_WeatherRegimes/WeatherRegimes and CST_RegimeAssign/RegimeAssign + PlotPDFsOLE plots two probability density gaussian functions and the optimal linear estimation + + CST_RFTemp/RF_Temp functions available for downscaling temperature - Fixes + CST_Anomaly handles exp, obs or both diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R new file mode 100644 index 0000000000000000000000000000000000000000..cebac85aa555252b19ba156d7d2d4e2a60f946a3 --- /dev/null +++ b/R/CST_RFTemp.R @@ -0,0 +1,563 @@ +#' @rdname CST_RFTemp +#' @title Temperature downscaling of a CSTools object using lapse rate +#' correction or a reference field +#' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +#' @description This function implements a simple lapse rate correction of a +#' temperature field (an object of class 's2dv_cube' as provided by +#' `CST_Load`) as input. +#' The input lon grid must be increasing (but can be modulo 360). +#' The input lat grid can be irregularly spaced (e.g. a Gaussian grid) +#' The output grid can be irregularly spaced in lon and/or lat. +#' @references Method described in ERA4CS MEDSCOPE milestone M3.2: +#' High-quality climate prediction data available to WP4 +#' [https://www.medscope-project.eu/the-project/deliverables-reports/]([https://www.medscope-project.eu/the-project/deliverables-reports/) +#' and in H2020 ECOPOTENTIAL Deliverable No. 8.1: +#' High resolution (1-10 km) climate, land use and ocean change scenarios +#' [https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf](https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf) +#' @param data An object of the class 's2dv_cube' as returned by `CST_Load`, +#' containing the temperature fields to downscale. +#' The data object is expected to have an element named \code{$data} +#' with at least two spatial dimensions named "lon" and "lat". +#' (these default names can be changed with the \code{lon_dim} and +#' \code{lat_dim} parameters) +#' @param oro An object of the class 's2dv_cube' as returned by `CST_Load`, +#' containing fine scale orography (in meters). +#' The destination downscaling area must be contained in the orography field. +#' @param xlim vector with longitude bounds for downscaling; +#' the full input field is downscaled if `xlim` and `ylim` are not specified. +#' @param ylim vector with latitude bounds for downscaling +#' @param lapse float with environmental lapse rate +#' @param lon_dim string with name of longitude dimension +#' @param lat_dim string with name of latitude dimension +#' @param time_dim a vector of character string indicating the name of temporal dimension. By default, it is set to NULL and it considers "ftime", "sdate" and "time" as temporal dimensions. +#' @param verbose logical if to print diagnostic output +#' @param nolapse logical, if true `oro` is interpreted as a fine-scale +#' climatology and used directly for bias correction +#' @param compute_delta logical if true returns only a delta to be used for +#' out-of-sample forecasts. Returns an object of the class 's2dv_cube', +#' containing a delta. Activates `nolapse = TRUE`. +#' @param delta An object of the class 's2dv_cube', containing a delta +#' to be applied to the downscaled input data. Activates `nolapse = TRUE`. +#' The grid of this object must coincide with that of the required output. +#' @param method string indicating the method used for interpolation: +#' "nearest" (nearest neighbours followed by smoothing with a circular +#' uniform weights kernel), "bilinear" (bilinear interpolation) +#' The two methods provide similar results, but nearest is slightly better +#' provided that the fine-scale grid is correctly centered as a subdivision +#' of the large-scale grid +#' @return CST_RFTemp() returns a downscaled CSTools object +#' (i.e., of the class 's2dv_cube'). +#' @export +#' @import multiApply +#' @examples +#' # Generate simple synthetic data and downscale by factor 4 +#' t <- rnorm(7 * 6 * 2 * 3 * 4)*10 + 273.15 + 10 +#' dim(t) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 6, lon = 7) +#' lon <- seq(3, 9, 1) +#' lat <- seq(42, 47, 1) +#' exp <- list(data = t, lat = lat, lon = lon) +#' attr(exp, 'class') <- 's2dv_cube' +#' o <- runif(29*29)*3000 +#' dim(o) <- c(lat = 29, lon = 29) +#' lon <- seq(3, 10, 0.25) +#' lat <- seq(41, 48, 0.25) +#' oro <- list(data = o, lat = lat, lon = lon) +#' attr(oro, 'class') <- 's2dv_cube' +#' res <- CST_RFTemp(exp, oro, xlim=c(4,8), ylim=c(43, 46), lapse=6.5) + +CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse = 6.5, + lon_dim = "lon", lat_dim = "lat", time_dim = NULL, + nolapse = FALSE, verbose = FALSE, compute_delta = FALSE, + method = "bilinear", delta = NULL) { + + if (!inherits(data, "s2dv_cube")) { + stop("Parameter 'data' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (!inherits(oro, "s2dv_cube")) { + stop("Parameter 'oro' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (!is.null(delta)) { + if (!inherits(delta, "s2dv_cube")) { + stop("Parameter 'delta' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + } + + res <- RFTemp(data$data, data$lon, data$lat, + oro$data, oro$lon, oro$lat, xlim, ylim, lapse, + lon_dim = lon_dim, lat_dim = lat_dim, time_dim = time_dim, + nolapse = nolapse, verbose = verbose, method = method, + compute_delta = compute_delta, delta = delta$data) + + data$data <- res$data + data$lon <- res$lon + data$lat <- res$lat + + return(data) +} + +#' @rdname RFTemp +#' @title Temperature downscaling of a CSTools object using lapse rate +#' correction (reduced version) +#' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +#' @description This function implements a simple lapse rate correction of a +#' temperature field (a multidimensional array) as input. +#' The input lon grid must be increasing (but can be modulo 360). +#' The input lat grid can be irregularly spaced (e.g. a Gaussian grid) +#' The output grid can be irregularly spaced in lon and/or lat. +#' @references Method described in ERA4CS MEDSCOPE milestone M3.2: +#' High-quality climate prediction data available to WP4 +#' [https://www.medscope-project.eu/the-project/deliverables-reports/]([https://www.medscope-project.eu/the-project/deliverables-reports/) +#' and in H2020 ECOPOTENTIAL Deliverable No. 8.1: +#' High resolution (1-10 km) climate, land use and ocean change scenarios +#' [https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf](https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf) +#' @param data Temperature array to downscale. +#' The input array is expected to have at least two dimensions named +#' "lon" and "lat" by default +#' (these default names can be changed with the \code{lon_dim} and +#' \code{lat_dim} parameters) +#' @param lon Vector or array of longitudes. +#' @param lat Vector or array of latitudes. +#' @param lonoro Vector or array of longitudes corresponding to the fine orography. +#' @param latoro Vector or array of latitudes corresponding to the fine orography. +#' @param oro Array containing fine-scale orography (in m) +#' The destination downscaling area must be contained in the orography field. +#' @param xlim vector with longitude bounds for downscaling; +#' the full input field is downscaled if `xlim` and `ylim` are not specified. +#' @param ylim vector with latitude bounds for downscaling +#' @param lapse float with environmental lapse rate +#' @param lon_dim string with name of longitude dimension +#' @param lat_dim string with name of latitude dimension +#' @param time_dim a vector of character string indicating the name of temporal dimension. By default, it is set to NULL and it considers "ftime", "sdate" and "time" as temporal dimensions. +#' @param verbose logical if to print diagnostic output +#' @param nolapse logical, if true `oro` is interpreted as a +#' fine-scale climatology and used directly for bias correction +#' @param compute_delta logical if true returns only a delta to be used for +#' out-of-sample forecasts. +#' @param delta matrix containing a delta to be applied to the downscaled +#' input data. The grid of this matrix is supposed to be same as that of +#' the required output field +#' @param method string indicating the method used for interpolation: +#' "nearest" (nearest neighbours followed by smoothing with a circular +#' uniform weights kernel), "bilinear" (bilinear interpolation) +#' The two methods provide similar results, but nearest is slightly better +#' provided that the fine-scale grid is correctly centered as a subdivision +#' of the large-scale grid +#' @return CST_RFTemp() returns a downscaled CSTools object +#' @return RFTemp() returns a list containing the fine-scale +#' longitudes, latitudes and the downscaled fields. +#' @export +#' @import multiApply +#' @examples +#' # Generate simple synthetic data and downscale by factor 4 +#' t <- rnorm(7 * 6 * 4 * 3) * 10 + 273.15 + 10 +#' dim(t) <- c(sdate = 3, ftime = 4, lat = 6, lon = 7) +#' lon <- seq(3, 9, 1) +#' lat <- seq(42, 47, 1) +#' o <- runif(29 * 29) * 3000 +#' dim(o) <- c(lat = 29, lon = 29) +#' lono <- seq(3, 10, 0.25) +#' lato <- seq(41, 48, 0.25) +#' res <- RFTemp(t, lon, lat, o, lono, lato, xlim = c(4, 8), ylim = c(43, 46), +#' lapse = 6.5) + +RFTemp <- function(data, lon, lat, oro, lonoro, latoro, + xlim = NULL, ylim = NULL, lapse = 6.5, + lon_dim = "lon", lat_dim = "lat", time_dim = NULL, + nolapse = FALSE, verbose = FALSE, compute_delta = FALSE, + method = "bilinear", delta = NULL) { + + # Check/detect time_dim + if (is.null(time_dim)) { + time_dim_names <- c("ftime", "sdate", "time") + time_dim_num <- which(time_dim_names %in% names(dim(data))) + if (length(time_dim_num) > 0) { + # Find time dimension with length > 0 + ilong <- which(dim(data)[time_dim_names[time_dim_num]] > 0) + if (length(ilong) > 0) { + time_dim <- time_dim_names[time_dim_num[ilong]] + } else { + stop("No time dimension longer than zero found.") + } + } else { + stop("Could not automatically detect a target time dimension ", + "in the provided data in 'data'.") + } + warning(paste("Selected time dim:", time_dim, "\n")) + } + + # Repeatedly apply .downscalet + if (is.null(delta)) { + result <- Apply(data, target_dims = c(lon_dim, lat_dim, time_dim), + fun = .downscalet, lon, lat, oro, lonoro, latoro, + xlim = xlim, ylim = ylim, lapse = lapse, + nolapse = nolapse, verbose = verbose, method = method, + compute_delta = compute_delta) + } else { + result <- Apply(list(data, delta), + target_dims = list(c(lon_dim, lat_dim, time_dim), + c(lon_dim, lat_dim)), + fun = .downscalet_delta, lon, lat, oro, lonoro, latoro, + xlim = xlim, ylim = ylim, lapse = lapse, + nolapse = nolapse, verbose = verbose, method = method, + compute_delta = compute_delta) + } + + names(dim(result$data))[1] <- names(dim(data))[names(dim(data)) == lon_dim] + names(dim(result$data))[2] <- names(dim(data))[names(dim(data)) == lat_dim] + result$lon <- array(result$lon[1:dim(result$lon)[1]]) + result$lat <- array(result$lat[1:dim(result$lat)[1]]) + names(dim(result$lon)) <- lon_dim + names(dim(result$lat)) <- lat_dim + return(result) +} + +#' Lapse-rate temperature correction downscaling +#' +#' @description Downscales a temperature field using a lapse-rate +#' correction based on a reference orography. Time-averaging is done on all +#' dimensions after the first two. +#' @author Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it} +#' @param lon vector of input longitudes +#' @param lat vector of input latitudes +#' @param t matrix of input temperature data +#' @param lono vector of orography longitudes +#' @param lato vector of orography latitudes +#' @param oro matrix of topographical elevations (in meters) +#' The destination downscaling area must be contained in the orography field. +#' @param xlim vector of longitude bounds; the full input field is downscaled if `xlim` and `ylim` are not specified. +#' @param ylim vector of latitude bounds +#' @param radius smoothing radius expressed in longitude units +#' (default is half a large-scale pixel) +#' @param lapse environmental lapse rate (in K/Km) +#' @param nolapse logical, if true `oro` is interpreted as a fine-scale +#' climatology and used directly for bias correction +#' @param compute_delta logical if true returns only a delta to be used for +#' out-of-sample forecasts. +#' @param delta matrix containing a delta to be applied to the input data. +#' The grid of this matrix is supposed to be same as +#' that of the required output field +#' @param verbose logical if to print diagnostic output +#' @return A downscaled temperature matrix +#' @examples +#' lon = 5:20 +#' lat = 35:40 +#' t = runif(16 * 6); dim(t) = c(16, 6) +#' lono = seq(5, 20, 0.1) +#' lato = seq(35, 40, 0.1) +#' o = runif(151 * 51) * 2000; dim(o) = c(151, 51) +#' td = .downscalet(t, lon, lat, o, lono, lato, c(8, 12), c(36, 38)) +#' @noRd +.downscalet <- function(t, lon, lat, oro, lono, lato, + xlim = NULL, ylim = NULL, + radius = 0, lapse = 6.5, nolapse = FALSE, + verbose = FALSE, compute_delta = FALSE, + method = "bilinear", delta = NULL) { + + if (!is.null(delta) & compute_delta) { + stop("Cannot `compute_delta` and provide `delta` at the same time.") + } + + tdim <- dim(t) + ldim <- FALSE + if (length(tdim) == 3) { + nt <- tdim[3] + } else if (length(tdim) == 2) { + nt <- 1 + dim(t) <- c(tdim, 1) + } else if (length(tdim) < 2) { + stop("Input array must have at least two dimensions") + } else { + ldim <- TRUE + dim(t) <- c(tdim[1:2], time = prod(tdim[-(1:2)])) + nt <- dim(t)[3] + } + if (lon[2] <= lon[1]) { + stop("Longitudes must be monotone increasing.") + } + # Regularize lon coords to monotone increasing + lon[lon >= lon[1]] <- lon[lon >= lon[1]] - 360 + if (lon[length(lon)] < 0) { lon <- lon + 360 } + lono[lono >= lono[1]] <- lono[lono >= lono[1]] - 360 + if (lono[length(lono)] < 0) { lono <- lono + 360 } + + dxl <- (lon[2] - lon[1]) / 2 + dyl <- (lat[2] - lat[1]) / 2 + if (radius == 0) { + radius <- dxl + } + if (is.null(xlim)) { + xlim <- c(lon[1] + dxl, lon[length(lon)] - dxl) + } + if (is.null(ylim)) { + ylim <- c(lat[1] + dyl, lat[length(lat)] - dyl) + if (ylim[1] > ylim[2]) { + ylim <- ylim[2:1] + } + } + #Add buffer + lon1 <- xlim[1] - radius + lon2 <- xlim[2] + radius + lat1 <- ylim[1] - radius + lat2 <- ylim[2] + radius + + orocut <- oro[(lono <= lon2) & (lono >= lon1), + (lato <= lat2) & (lato >= lat1)] + lonocut <- lono[(lono <= lon2) & (lono >= lon1)] + latocut <- lato[(lato <= lat2) & (lato >= lat1)] + + dxol <- lonocut[2] - lonocut[1] + nrad <- as.integer(radius / abs(dxol) + 0.5) + + if (length(lonocut) == 0 | length(latocut) == 0) { + stop("Orography not available for selected area") + } + + if (is.null(delta) & compute_delta & nolapse) { + if (verbose) { + print("Time averaging") + } + # If we just need the delta we can work with time averages + t <- apply(t, c(1, 2), mean) + dim(t) <- c(dim(t), time = 1) + } + + if (!(is.null(delta) & compute_delta & !nolapse)) { + # This calculation is not needed if we just need the delta + # and lapse rate is used + if (verbose) { + print(paste("Interpolation using", method, "method")) + } + tout <- .interp2d(t, lon, lat, lonocut, latocut, method = method) + # Fine-scale smooth interpolated input field + if (method == "nearest") { + if (verbose) { + print(paste("Smoothing interpolated field")) + } + tout <- .smooth(tout, nrad) + } + } + if (is.null(delta)) { + # Compute delta + if (nolapse) { + # oro is a reference fine-scale field: use that directly + # for bias correcting the interpolated input field's temporal mean + t1 <- orocut - apply(tout, c(1, 2), mean) + } else { + # Lapse-rate correction + orocuts <- .smooth(orocut, nrad) + t1 <- -(orocut - orocuts) * lapse / 1000. + } + if (compute_delta) { + # Return delta + tout <- t1[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1])] + } else { + # Apply delta + if (verbose) { + print("Adding computed delta") + } + for (i in seq_len(nt)) { + tout[, , i] <- tout[, , i] + t1 + } + tout <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1]), ] + if (ldim) { + dim(tout) <- c(dim(tout)[1:2], tdim[-(1:2)]) + } + } + } else { + # Just apply the provided delta + tout <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1]), ] + if (any(dim(delta)[1:2] != dim(tout)[1:2])) { + stop("Provided delta has not the same dimensions as required output") + } + if (verbose) { + print("Adding provided delta") + } + for (i in seq_len(nt)) { + tout[, , i] <- tout[, , i] + delta + } + if (ldim) { + dim(tout) <- c(dim(tout)[1:2], tdim[-(1:2)]) + } + } + lonocut <- lonocut[(lonocut <= xlim[2]) & (lonocut >= xlim[1])] + latocut <- latocut[(latocut <= ylim[2]) & (latocut >= ylim[1])] + + return(list(data = tout, lon = lonocut, lat = latocut)) +} + +# Special driver version of .downscalet to apply delta +.downscalet_delta <- function(t, delta, lon, lat, oro, lono, lato, + xlim = NULL, ylim = NULL, radius = 0, + lapse = 6.5, nolapse = FALSE, verbose = FALSE, + compute_delta = FALSE, method = "bilinear") { + res <- .downscalet(t, lon, lat, oro, lono, lato, xlim = xlim, + ylim = ylim, radius = radius, lapse = lapse, + nolapse = nolapse, verbose = verbose, + compute_delta = compute_delta, delta = delta, + method = method) +} + +#' Nearest neighbour interpolation +#' +#' @description The input field is interpolated onto the output +#' coordinate grid using nearest neighbours or bilinear interpolation. +#' The input lon grid must be monotone increasing. +#' The input lat grid can be irregularly spaced (e.g. a Gaussian grid) +#' The output grid can be irregularly spaced in lon and/or lat. +#' @author Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it} +#' @param z matrix with the input field to interpolate (assumed to +#' include also a third time dimension) +#' @param lon vector of input longitudes +#' @param lat vector of input latitudes +#' @param lonp vector of output longitudes +#' @param latp vector of output latitudes +#' @param method string indicating the interpolation method +#' ("nearest" or "bilinear" (default)) +#' @return The interpolated field. +#' @examples +#' lon = 5:11 +#' lat = 35:40 +#' z = runif(7 * 6 * 2); dim(z) = c(7, 6, 2) +#' lonp = seq(5, 10, 0.2) +#' latp = seq(35, 40, 0.2) +#' zo <- .interp2d(z, lon, lat, lonp, latp, method = "nearest") +#' @noRd +.interp2d <- function(z, lon, lat, lonp, latp, method="bilinear") { + + nx <- length(lonp) + ny <- length(latp) + nt <- dim(z)[3] + # Interpolate nn assuming a regular grid + zo <- array(0., c(nx, ny, nt)) + dy <- lat[2] - lat[1] + dx <- lon[2] - lon[1] + + if (method == "nearest") { + jj <- 1:length(latp) + for (j in 1:length(latp)) { + jj[j] <- which.min(abs(latp[j]-lat)) + } + ii <- ((lonp - (lon[1] - dx / 2)) %/% dx + 1) + # If lon are global and periodic attempt to fix issues + if ((lon[1] - lon[length(lon)]) %% 360 == dx) { + ii[ii <= 0] <- ii[ii <= 0] + length(lon) + ii[ii > length(lon)] <- ii[ii > length(lon)] - length(lon) + } + if ((ii[1] <= 0) | (ii[length(ii)] > length(lon)) | + (jj[1] <= 0) | (jj[length(jj)] > length(lat))) { + stop("Downscaling area not contained in input data") + } + zo[, , ] <- z[ii, jj, ] + } else { + jj <- 1:length(latp) + jj2 <- jj + for (j in 1:length(latp)) { + jmin <- which.min(abs(latp[j]-lat)) + if ( (((latp[j]-lat[jmin]) >= 0) & ( dy >= 0)) | + (((latp[j]-lat[jmin]) < 0) & ( dy <= 0))) { + jj[j] <- jmin + jj2[j] <- jmin + 1 + } else { + jj2[j] <- jmin + jj[j] <- jmin - 1 + } + } + ii <- ((lonp - lon[1]) %/% dx + 1) + ii2 <- ii + 1 + # If lon are global and periodic attempt to fix issues + if ((lon[1] - lon[length(lon)]) %% 360 == dx) { + ii[ii <= 0] <- ii[ii <= 0] + length(lon) + ii[ii > length(lon)] <- ii[ii > length(lon)] - length(lon) + ii2[ii2 <= 0] <- ii2[ii2 <= 0] + length(lon) + ii2[ii2 > length(lon)] <- ii2[ii2 > length(lon)] - length(lon) + } + if ((ii[1] <= 0) | (ii[length(ii)] > length(lon)) | + (jj[1] <= 0) | (jj[length(jj)] > length(lat)) | + (ii2[1] <= 0) | (ii2[length(ii2)] > length(lon)) | + (jj2[1] <= 0) | (jj2[length(jj2)] > length(lat))) { + stop("Downscaling area not contained in input data") + } + xx <- ((lonp - lon[ii]) / dx) %% 360 + yy <- (latp - lat[jj]) / (lat[jj2] - lat[jj]) + xx <- xx[row(zo[, , 1])] + yy <- yy[col(zo[, , 1])] + dim(xx) <- c(nx, ny) + dim(yy) <- c(nx, ny) + w00 <- (1 - xx) * (1 - yy) + w10 <- xx * (1 - yy) + w01 <- (1 - xx) * yy + w11 <- xx * yy + for (k in seq_len(nt)) { + zo[, , k] <- z[ii, jj, k] * w00 + z[ii2, jj, k] * w10 + + z[ii, jj2, k] * w01 + z[ii2, jj2, k] * w11 + } + } + names(dim(zo)) <- names(dim(z)) + return(zo) +} + +#' Smoothening using convolution with a circular kernel +#' +#' @description The input field is convolved with a circular kernel with equal +#' weights. Takes into account missing values. +#' @author Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it} +#' @param z matrix with the input field to smoothen, with dimensions `c(ns, ns)` +#' @param sdim the smoothing kernel radius in pixel +#' @return The smoothened field. +#' @examples +#' z <- rnorm(64 * 64) +#' dim(z) <- c(64, 64) +#' zs <- smooth(z, 8) +#' sd(zs) +#' # [1] 0.1334648 +#' @noRd +.smooth <- function(z, sdim) { + nsx <- dim(z)[1] + nsy <- dim(z)[2] + zdim <- dim(z) + if (length(dim(z)) == 2) { + dim(z) <- c(dim(z), time = 1) + nt <- 1 + } else { + nt <- dim(z)[3] + } + imask <- !is.finite(z) + z[imask] <- 0. + mask <- matrix(0., nsx, nsy) + for (i in 1:nsx) { + for (j in 1:nsy) { + kx <- i - 1 + ky <- j - 1 + if (i > nsx / 2 + 1) { + kx <- i - nsx - 1 + } + if (j > nsy / 2 + 1) { + ky <- j - nsy - 1 + } + r2 <- kx * kx + ky * ky + mask[i, j] <- (r2 <= (sdim * sdim)) + } + } + fm <- fft(mask) + zf <- array(0., c(nsx, nsy, nt)) + for (k in seq_len(nt)) { + zf[, , k] <- Re(fft(fm * fft(z[, , k]), inverse = TRUE) + ) / sum(mask) / length(fm) + if (sum(imask[, , k]) > 0) { + zz <- z[, , k] + zz[!imask[, , k]] <- 1.0 + zz <- zf[, , k] / (Re(fft(fm * fft(zz), inverse = TRUE)) / + sum(mask) / length(fm)) + zz[imask[, , k]] <- NA + zf[, , k] <- zz + } + } + dim(zf) <- zdim + return(zf) +} diff --git a/man/CST_RFTemp.Rd b/man/CST_RFTemp.Rd new file mode 100644 index 0000000000000000000000000000000000000000..8ab5b6f329b63174cab568199ee4e5245b8927e5 --- /dev/null +++ b/man/CST_RFTemp.Rd @@ -0,0 +1,107 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_RFTemp.R +\name{CST_RFTemp} +\alias{CST_RFTemp} +\title{Temperature downscaling of a CSTools object using lapse rate +correction or a reference field} +\usage{ +CST_RFTemp( + data, + oro, + xlim = NULL, + ylim = NULL, + lapse = 6.5, + lon_dim = "lon", + lat_dim = "lat", + time_dim = NULL, + nolapse = FALSE, + verbose = FALSE, + compute_delta = FALSE, + method = "bilinear", + delta = NULL +) +} +\arguments{ +\item{data}{An object of the class 's2dv_cube' as returned by `CST_Load`, +containing the temperature fields to downscale. +The data object is expected to have an element named \code{$data} +with at least two spatial dimensions named "lon" and "lat". +(these default names can be changed with the \code{lon_dim} and +\code{lat_dim} parameters)} + +\item{oro}{An object of the class 's2dv_cube' as returned by `CST_Load`, +containing fine scale orography (in meters). +The destination downscaling area must be contained in the orography field.} + +\item{xlim}{vector with longitude bounds for downscaling; +the full input field is downscaled if `xlim` and `ylim` are not specified.} + +\item{ylim}{vector with latitude bounds for downscaling} + +\item{lapse}{float with environmental lapse rate} + +\item{lon_dim}{string with name of longitude dimension} + +\item{lat_dim}{string with name of latitude dimension} + +\item{time_dim}{a vector of character string indicating the name of temporal dimension. By default, it is set to NULL and it considers "ftime", "sdate" and "time" as temporal dimensions.} + +\item{nolapse}{logical, if true `oro` is interpreted as a fine-scale +climatology and used directly for bias correction} + +\item{verbose}{logical if to print diagnostic output} + +\item{compute_delta}{logical if true returns only a delta to be used for +out-of-sample forecasts. Returns an object of the class 's2dv_cube', +containing a delta. Activates `nolapse = TRUE`.} + +\item{method}{string indicating the method used for interpolation: +"nearest" (nearest neighbours followed by smoothing with a circular +uniform weights kernel), "bilinear" (bilinear interpolation) +The two methods provide similar results, but nearest is slightly better +provided that the fine-scale grid is correctly centered as a subdivision +of the large-scale grid} + +\item{delta}{An object of the class 's2dv_cube', containing a delta +to be applied to the downscaled input data. Activates `nolapse = TRUE`. +The grid of this object must coincide with that of the required output.} +} +\value{ +CST_RFTemp() returns a downscaled CSTools object +(i.e., of the class 's2dv_cube'). +} +\description{ +This function implements a simple lapse rate correction of a +temperature field (an object of class 's2dv_cube' as provided by +`CST_Load`) as input. +The input lon grid must be increasing (but can be modulo 360). +The input lat grid can be irregularly spaced (e.g. a Gaussian grid) +The output grid can be irregularly spaced in lon and/or lat. +} +\examples{ +# Generate simple synthetic data and downscale by factor 4 +t <- rnorm(7 * 6 * 2 * 3 * 4)*10 + 273.15 + 10 +dim(t) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 6, lon = 7) +lon <- seq(3, 9, 1) +lat <- seq(42, 47, 1) +exp <- list(data = t, lat = lat, lon = lon) +attr(exp, 'class') <- 's2dv_cube' +o <- runif(29*29)*3000 +dim(o) <- c(lat = 29, lon = 29) +lon <- seq(3, 10, 0.25) +lat <- seq(41, 48, 0.25) +oro <- list(data = o, lat = lat, lon = lon) +attr(oro, 'class') <- 's2dv_cube' +res <- CST_RFTemp(exp, oro, xlim=c(4,8), ylim=c(43, 46), lapse=6.5) +} +\references{ +Method described in ERA4CS MEDSCOPE milestone M3.2: +High-quality climate prediction data available to WP4 +[https://www.medscope-project.eu/the-project/deliverables-reports/]([https://www.medscope-project.eu/the-project/deliverables-reports/) +and in H2020 ECOPOTENTIAL Deliverable No. 8.1: +High resolution (1-10 km) climate, land use and ocean change scenarios +[https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf](https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf) +} +\author{ +Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/man/RFTemp.Rd b/man/RFTemp.Rd new file mode 100644 index 0000000000000000000000000000000000000000..106ae6e218f2037ac5987d16e49641ac90c044d8 --- /dev/null +++ b/man/RFTemp.Rd @@ -0,0 +1,114 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_RFTemp.R +\name{RFTemp} +\alias{RFTemp} +\title{Temperature downscaling of a CSTools object using lapse rate +correction (reduced version)} +\usage{ +RFTemp( + data, + lon, + lat, + oro, + lonoro, + latoro, + xlim = NULL, + ylim = NULL, + lapse = 6.5, + lon_dim = "lon", + lat_dim = "lat", + time_dim = NULL, + nolapse = FALSE, + verbose = FALSE, + compute_delta = FALSE, + method = "bilinear", + delta = NULL +) +} +\arguments{ +\item{data}{Temperature array to downscale. +The input array is expected to have at least two dimensions named +"lon" and "lat" by default +(these default names can be changed with the \code{lon_dim} and +\code{lat_dim} parameters)} + +\item{lon}{Vector or array of longitudes.} + +\item{lat}{Vector or array of latitudes.} + +\item{oro}{Array containing fine-scale orography (in m) +The destination downscaling area must be contained in the orography field.} + +\item{lonoro}{Vector or array of longitudes corresponding to the fine orography.} + +\item{latoro}{Vector or array of latitudes corresponding to the fine orography.} + +\item{xlim}{vector with longitude bounds for downscaling; +the full input field is downscaled if `xlim` and `ylim` are not specified.} + +\item{ylim}{vector with latitude bounds for downscaling} + +\item{lapse}{float with environmental lapse rate} + +\item{lon_dim}{string with name of longitude dimension} + +\item{lat_dim}{string with name of latitude dimension} + +\item{time_dim}{a vector of character string indicating the name of temporal dimension. By default, it is set to NULL and it considers "ftime", "sdate" and "time" as temporal dimensions.} + +\item{nolapse}{logical, if true `oro` is interpreted as a +fine-scale climatology and used directly for bias correction} + +\item{verbose}{logical if to print diagnostic output} + +\item{compute_delta}{logical if true returns only a delta to be used for +out-of-sample forecasts.} + +\item{method}{string indicating the method used for interpolation: +"nearest" (nearest neighbours followed by smoothing with a circular +uniform weights kernel), "bilinear" (bilinear interpolation) +The two methods provide similar results, but nearest is slightly better +provided that the fine-scale grid is correctly centered as a subdivision +of the large-scale grid} + +\item{delta}{matrix containing a delta to be applied to the downscaled +input data. The grid of this matrix is supposed to be same as that of +the required output field} +} +\value{ +CST_RFTemp() returns a downscaled CSTools object + +RFTemp() returns a list containing the fine-scale +longitudes, latitudes and the downscaled fields. +} +\description{ +This function implements a simple lapse rate correction of a +temperature field (a multidimensional array) as input. +The input lon grid must be increasing (but can be modulo 360). +The input lat grid can be irregularly spaced (e.g. a Gaussian grid) +The output grid can be irregularly spaced in lon and/or lat. +} +\examples{ +# Generate simple synthetic data and downscale by factor 4 +t <- rnorm(7 * 6 * 4 * 3) * 10 + 273.15 + 10 +dim(t) <- c(sdate = 3, ftime = 4, lat = 6, lon = 7) +lon <- seq(3, 9, 1) +lat <- seq(42, 47, 1) +o <- runif(29 * 29) * 3000 +dim(o) <- c(lat = 29, lon = 29) +lono <- seq(3, 10, 0.25) +lato <- seq(41, 48, 0.25) +res <- RFTemp(t, lon, lat, o, lono, lato, xlim = c(4, 8), ylim = c(43, 46), + lapse = 6.5) +} +\references{ +Method described in ERA4CS MEDSCOPE milestone M3.2: +High-quality climate prediction data available to WP4 +[https://www.medscope-project.eu/the-project/deliverables-reports/]([https://www.medscope-project.eu/the-project/deliverables-reports/) +and in H2020 ECOPOTENTIAL Deliverable No. 8.1: +High resolution (1-10 km) climate, land use and ocean change scenarios +[https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf](https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf) +} +\author{ +Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/tests/testthat/test-CST_RFTemp.R b/tests/testthat/test-CST_RFTemp.R new file mode 100644 index 0000000000000000000000000000000000000000..536e44f0e0092d0266cf253271c5df6138c7666e --- /dev/null +++ b/tests/testthat/test-CST_RFTemp.R @@ -0,0 +1,71 @@ +context("Generic tests") +test_that("Sanity checks and simple use cases", { + # Generate simple synthetic data + t <- rnorm(2 * 6 * 6 * 2 * 3 * 4) * 10 + 273.15 + 10 + dim(t) <- c(dataset = 2, member = 2, sdate = 3, ftime = 4, lat = 6, lon = 6) + lon <- seq(4, 9, 1) + lat <- seq(42, 47, 1) + exp <- list(data = t, lat = lat, lon = lon) + o <- runif(29 * 29) * 3000 + dim(o) <- c(lat = 29, lon = 29) + lon <- seq(3.125, 10.125, 0.25) - 100 + lat <- seq(41.125, 48.125, 0.25) - 60 + oro <- list(data = o, lat = lat, lon = lon) + attr(oro, "class") <- "s2dv_cube" + + expect_error( + res <- CST_RFTemp(exp, oro, xlim = c(1, 3), ylim = c(1, 3)), + paste("Parameter 'data' must be of the class", + "'s2dv_cube', as output by CSTools::CST_Load.")) + attr(exp, "class") <- "s2dv_cube" + + expect_error( + res <- CST_RFTemp(exp, oro, xlim = c(1, 3), ylim = c(1, 3)), + "Orography not available for selected area" + ) + + oro$lon <- oro$lon + 100 + oro$lat <- oro$lat + 60 + + expect_error( + res <- CST_RFTemp(exp, oro, xlim = c(3, 8), ylim = c(43, 46)), + "Downscaling area not contained in input data" + ) + + expect_warning( + resl <- CST_RFTemp(exp, oro, lapse = 6.5), + "Selected time dim: ftime" + ) + + dimexp <- dim(exp$data) + expect_equal(dim(resl$data), c(lon = 16, lat = 16, dimexp["ftime"], + dimexp["sdate"], dimexp["dataset"], + dimexp["member"])) + expect_warning( + resd <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), + nolapse = TRUE), NA) + dim(resd$data) <- c(16, 16, 4 * 3 * 2 * 2) + resm <- apply(resd$data, c(1, 2), mean) + expect_equal(mean(oro$data[7:22, 7:22]), mean(resm), tolerance = 1e-10) + + # Test precomputation of delta based with nolapse + delta <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), + nolapse = TRUE, compute_delta = TRUE) + resdnl <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), + nolapse = TRUE, delta = delta) + expect_equal(mean(resdnl$data), mean(resd$data), tolerance = 1e-10) + + # Test precomputation of delta based with lapse rate + delta <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), + lapse = 6.5, compute_delta = TRUE) + resdl <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), + delta = delta) + expect_equal(mean(resdl$data), mean(resl$data), tolerance = 1e-10) + + expect_warning( + resd <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), + nolapse = TRUE, method = "nearest"), NA) + dim(resd$data) <- c(16, 16, 4 * 3 * 2 * 2) + resm <- apply(resd$data, c(1, 2), mean) + expect_equal(mean(oro$data[7:22, 7:22]), mean(resm), tolerance = 1e-10) +}) diff --git a/tests/testthat/test-CST_WeatherRegimes.R b/tests/testthat/test-CST_WeatherRegimes.R index 33e75e24768c839b4c3458fe06f742208f86330e..e2bc6feb989f2c11aa7ecae275950b288e2aaefa 100644 --- a/tests/testthat/test-CST_WeatherRegimes.R +++ b/tests/testthat/test-CST_WeatherRegimes.R @@ -13,8 +13,8 @@ test_that("Sanity checks", { paste0("Parameter 'data' must be an array with named dimensions.")) data1 <- 1 : 20 - dim(data1) <- c(lat = 5, lon=4) - data1 <- list(data = data1 , lat=1:5) + dim(data1) <- c(lat = 5, lon = 4) + data1 <- list(data = data1 , lat = 1:5) class(data1) <- 's2dv_cube' expect_error( CST_WeatherRegimes(data = data1), @@ -29,38 +29,38 @@ test_that("Sanity checks", { paste0("Parameter 'lat' must be specified.")) data1 <- 1 : 400 - dim(data1) <- c(time = 20, lat = 5, lon=4) - data1 <- list(data = data1, lat=1:5) + dim(data1) <- c(time = 20, lat = 5, lon = 4) + data1 <- list(data = data1, lat = 1:5) class(data1) <- 's2dv_cube' expect_error( CST_WeatherRegimes(data = data1), paste0("Parameter 'ncenters' must be specified.")) expect_error( - CST_WeatherRegimes(data = data1, ncenters=3), + CST_WeatherRegimes(data = data1, ncenters = 3), paste0("Parameter 'lon' must be specified.")) expect_equal( - names(dim(CST_WeatherRegimes(data = data1, ncenters=3, EOFS= FALSE)$data)), + names(dim(CST_WeatherRegimes(data = data1, ncenters = 3, EOFS = FALSE)$data)), c('lat', 'lon', 'cluster')) data1 <- 1 : 400 - dim(data1) <- c(sdate = 2, ftime = 10, lat = 5, lon=4) - data1 <- list(data = data1, lat=1:5) + dim(data1) <- c(sdate = 2, ftime = 10, lat = 5, lon = 4) + data1 <- list(data = data1, lat = 1:5) class(data1) <- 's2dv_cube' nclusters <- 3 expect_equal( dim(CST_WeatherRegimes(data = data1 , ncenters = nclusters, - EOFS = FALSE)$statistics$frequency),c(2, nclusters)) + EOFS = FALSE)$statistics$frequency), c(2, nclusters)) expect_equal( - names(dim(CST_WeatherRegimes(data = data1, nclusters, EOFS= FALSE)$data)), + names(dim(CST_WeatherRegimes(data = data1, nclusters, EOFS = FALSE)$data)), c('lat', 'lon', 'cluster')) data1 <- 1 : 400 - dim(data1) <- c(sdate = 2, ftime = 10, lat = 5, lon=4) - data1 <- list(data = data1, lat=1:5 ,lon=1:4) + dim(data1) <- c(sdate = 2, ftime = 10, lat = 5, lon = 4) + data1 <- list(data = data1, lat = 1:5 ,lon = 1:4) class(data1) <- 's2dv_cube' expect_equal( @@ -68,32 +68,32 @@ test_that("Sanity checks", { c('pvalue', 'cluster', 'frequency', 'persistence')) expect_equal( - names(CST_WeatherRegimes(data = data1 , ncenters = 4, method='ward.D')$statistics), + names(CST_WeatherRegimes(data = data1 , ncenters = 4, method = 'ward.D')$statistics), c('pvalue', 'cluster')) expect_equal( - names(dim(CST_WeatherRegimes(data = data1, ncenters=4)$data)), + names(dim(CST_WeatherRegimes(data = data1, ncenters = 4)$data)), c('lat', 'lon', 'cluster')) data1 <- 1 : 400 - dim(data1) <- c(time = 20, lat = 5, lon=4) + dim(data1) <- c(time = 20, lat = 5, lon = 4) data1[4,,] <- NA - data1 <- list(data = data1, lat=1:5 ,lon=1:4) + data1 <- list(data = data1, lat = 1:5 ,lon = 1:4) class(data1) <- 's2dv_cube' expect_error( - CST_WeatherRegimes(data = data1, ncenters=3, EOFS = FALSE), + CST_WeatherRegimes(data = data1, ncenters = 3, EOFS = FALSE), paste0("Parameter 'data' contains NAs in the 'time' dimensions.")) data1 <- 1 : 400 - dim(data1) <- c(time = 20, lat = 5, lon=4) + dim(data1) <- c(time = 20, lat = 5, lon = 4) data1[,2,3] <- NA - data1 <- list(data = data1, lat=1:5 ,lon=1:4) + data1 <- list(data = data1, lat = 1:5 ,lon = 1:4) class(data1) <- 's2dv_cube' expect_equal( - any(is.na(CST_WeatherRegimes(data = data1, ncenters=3, EOFS = FALSE)$data)), + any(is.na(CST_WeatherRegimes(data = data1, ncenters =3, EOFS = FALSE)$data)), TRUE) expect_equal( - names(dim(CST_WeatherRegimes(data = data1, ncenters=3, EOFS = FALSE)$data)), + names(dim(CST_WeatherRegimes(data = data1, ncenters =3, EOFS = FALSE)$data)), c('lat', 'lon', 'cluster')) }) diff --git a/vignettes/ENSclustering_vignette.Rmd b/vignettes/ENSclustering_vignette.Rmd index 52707d73a6dd63b26c2e6ec10728fb268607b569..af8d5c6fd2fac1820509ff1426bef452bf861379 100644 --- a/vignettes/ENSclustering_vignette.Rmd +++ b/vignettes/ENSclustering_vignette.Rmd @@ -14,7 +14,7 @@ Ensemble clustering ### Introduction Ensemble forecasts provide insight on average weather conditions from sub-seasonal to seasonal timescales. With large ensembles, it is useful to group members according to similar characteristics and to select the most representative member of each cluster. This allows to characterize forecast scenarios in a multi-model (or single model) ensemble prediction. In particular, at the regional level, this function can be used to identify the subset of ensemble members that best represent the full range of possible outcomes that may be fed in turn to downscaling applications. The choice of the ensemble members can be done by selecting different options within the function in order to meet the requirements of specific climate information products that can be tailored to different regions and user needs. -For a detailed description of the tool please refer to the CSTools guide : +For a detailed description of the tool please refer to the CSTools guide : ### Steps of the vignette