From 53d8920494c72896fedec891173d2feef69c2893 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Wed, 26 Feb 2020 00:19:43 +0100 Subject: [PATCH 01/41] Working version - no dependences on new rainfarmr --- R/CST_RFTemp.R | 212 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 212 insertions(+) create mode 100644 R/CST_RFTemp.R diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R new file mode 100644 index 00000000..6ebfb381 --- /dev/null +++ b/R/CST_RFTemp.R @@ -0,0 +1,212 @@ +#' @rdname CST_RFTemp +#' @title Temperature downscaling of a CSTools object using lapse rate correction +#' @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. +#' @references Terzago, S. et al. (2018). NHESS 18(11), 2825-2840 http://doi.org/10.5194/nhess-18-2825-2018 +#' @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 +#' The number of longitudes and latitudes in the input data is expected to be even and the same. If not +#' spatial dimensions named "lon" and "lat". +#' @param oro An object of the class 's2dv_cube' as returned by `CST_Load`, +#' containing fine scale orography (in meters). +#' @param xlim vector with longitude bounds for downscaling +#' @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 +#' @return CST_RFTemp() returns a downscaled CSTools object (i.e., of the class 's2dv_cube'). +#' @export +#' @import multiApply +#' @import rainfarmr +#' @examples +#' + +CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, lon_dim = "lon", lat_dim = "lat") { + + 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) + + 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 (an object of class 's2dv_cube' as provided by `CST_Load`) as input. +#' @references Terzago, S. et al. (2018). NHESS 18(11), 2825-2840 http://doi.org/10.5194/nhess-18-2825-2018 +#' @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 oro Array containing fine-scale orography (in m) +#' @param xlim vector with longitude bounds for downscaling +#' @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 +#' @return RFTemp() returns a list containing the fine-scale longitudes, latitudes and the downscaled fields. +#' @export +#' @import multiApply +#' @import rainfarmr +#' @examples +#' + +RFTemp <- function(data, lon, lat, oro, lonoro, latoro, + xlim, ylim, lapse = 6.5, lon_dim = "lon", lat_dim = "lat") { + + # Repeatedly apply .downscalet + result <- Apply(data, target_dims=c(lon_dim, lat_dim), fun=.downscalet, + lon, lat, oro, lonoro, latoro, xlim, ylim, lapse=lapse)$output1 + names(dim(result))[1] <- names(dim(tas))[names(dim(tas))==lon_dim] + names(dim(result))[2] <- names(dim(tas))[names(dim(tas))==lat_dim] + + return(list(data = result, lon = lonoro, lat = latoro)) +} + + +#' Lapse-rate temperature correction downscaling +#' +#' @description Downscales a temperature field using a lapse-rate +#' correction based on a reference orography +#' @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) +#' @param xlim vector of longitude bounds +#' @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) +#' @return A downscaled temperature matrix +#' @export +#' @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(lon, lat,t,lono,lato,o,c(8,12),c(36,38)) + +.downscalet <- function(t, lon, lat, oro, lono, lato, xlim, ylim, radius = 0, lapse = 6.5) { + + if (length(dim(t))==3) { + nt = dim(t)[3] + } else { + nt = 1 + dim(t) <- c(dim(t), 1) + } + + dxl <- lon[2] - lon[1] + if(radius==0) { + radius <- dxl/2 + } + #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)] + + if (length(dim(tcut0))<3) { + dim(tcut0) <- c(dim(tcut0), 1) + } + + nx <- length(lonocut) + ny <- length(latocut) + + # Interpolate nn assuming a regular grid + tcut <- array(0., c(nx, ny, nt)) + dy <- lat[2]-lat[1] + dx <- lon[2]-lon[1] + jj <- ((latocut-(lat[1]-dy/2)) %/% dy + 1) + ii <- ((lonocut-(lon[1]-dx/2)) %/% dx + 1) + + 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") + } + for (k in seq_len(nt)) { + for (j in seq_along(jj)) { + tcut[,j,k] <- t[ii,jj[j],k] + } + } + + dxol <- lonocut[2] - lonocut[1] + nrad <- as.integer(radius / abs(dxol) + 0.5) + + ns = length(lonocut) + orocuts=.smooth(orocut, nrad) + t1 <- -(orocut-orocuts)*lapse/1000. + + tout <- array(0., c(nx, ny, nt)) + for (i in seq_len(nt)) { + tout[,, i] <- .smooth(tcut[,,i], nrad) + t1 + } + downscalet <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1]),] +} + +#' 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. +#' @export +#' @examples +#' z <- rnorm(64 * 64) +#' dim(z) <- c(64, 64) +#' zs <- smooth(z, 8) +#' sd(zs) +#' # [1] 0.1334648 + +.smooth <- function(z, sdim) { + imask <- !is.finite(z) + z[imask] <- 0. + nsx <- dim(z)[1] + nsy <- dim(z)[2] + + 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] <- exp( - (r2 / (sdim * sdim)) / 2) + mask[i, j] <- (r2 <= (sdim * sdim)) + } + } + fm <- fft(mask) + zf <- Re(fft(fm * fft(z), inverse = TRUE)) / sum(mask) / length(fm) + if (sum(imask) > 0) { + z[!imask] <- 1.0 + zf <- zf / (Re(fft(fm * fft(z), inverse = TRUE)) / + sum(mask) / length(fm)) + } + zf[imask] <- NA + return(zf) +} -- GitLab From fcfc92ce5744635a2539785bbe47488fa04a8cc3 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Wed, 26 Feb 2020 12:07:23 +0100 Subject: [PATCH 02/41] added reference --- R/CST_RFTemp.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 6ebfb381..3615df05 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -3,7 +3,7 @@ #' @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. -#' @references Terzago, S. et al. (2018). NHESS 18(11), 2825-2840 http://doi.org/10.5194/nhess-18-2825-2018 +#' @references Method described 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) #' @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 @@ -41,7 +41,7 @@ CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, lon_dim = "lon", lat_di #' @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. -#' @references Terzago, S. et al. (2018). NHESS 18(11), 2825-2840 http://doi.org/10.5194/nhess-18-2825-2018 +#' @references Method described 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) #' @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) -- GitLab From b48db9b3348d38b8e53166987bb308831ca056fd Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Wed, 26 Feb 2020 15:50:33 +0100 Subject: [PATCH 03/41] minor doc improvements --- R/CST_RFTemp.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 3615df05..2a42a483 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -3,7 +3,7 @@ #' @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. -#' @references Method described 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) +#' @references Method described in H2020 ECOPOTENTIAL Deliverable No. 8.1: High resolution (1-10 km) climate, land use and ocean change scenarios \link[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 @@ -40,8 +40,8 @@ CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, lon_dim = "lon", lat_di #' @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 (an object of class 's2dv_cube' as provided by `CST_Load`) as input. -#' @references Method described 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) +#' temperature field (a multidimensional array) as input +#' @references Method described in H2020 ECOPOTENTIAL Deliverable No. 8.1: High resolution (1-10 km) climate, land use and ocean change scenarios \link[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) -- GitLab From 5f11e80197d77f037b728ee05498a11d6c034e22 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Wed, 26 Feb 2020 17:41:11 +0100 Subject: [PATCH 04/41] added test --- tests/testthat/test-CST_RFTemp.R | 40 ++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 tests/testthat/test-CST_RFTemp.R diff --git a/tests/testthat/test-CST_RFTemp.R b/tests/testthat/test-CST_RFTemp.R new file mode 100644 index 00000000..9d577d64 --- /dev/null +++ b/tests/testthat/test-CST_RFTemp.R @@ -0,0 +1,40 @@ +context("Generic tests") +test_that("Sanity checks and simple use cases", { + # Generate simple synthetic data + t <- rnorm(6 * 6 * 2 * 3 * 4)*10 + 273.15 + 10 + dim(t) <- c(dataset = 1, 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, 10, 0.25) - 100 + lat <- seq(41, 48, 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)), + "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" + ) + + res <- CST_RFTemp(exp, oro, xlim=c(4,8), ylim=c(43, 46), lapse=6.5) + + dimexp=dim(exp$data) + expect_equal(dim(res$data), c(lon=17,lat=13, + dimexp["dataset"], dimexp["member"], + dimexp["sdate"], dimexp["ftime"])) +}) -- GitLab From 062c6ef34fdde4b8ff6b90f7cca9171611e9f829 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Wed, 26 Feb 2020 17:41:48 +0100 Subject: [PATCH 05/41] undefined variables --- R/CST_RFTemp.R | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 2a42a483..0bc7d631 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -25,6 +25,11 @@ CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, lon_dim = "lon", lat_dim = "lat") { + if (!inherits(data, 's2dv_cube')) { + stop("Parameter 'data' 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) @@ -66,8 +71,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, # Repeatedly apply .downscalet result <- Apply(data, target_dims=c(lon_dim, lat_dim), fun=.downscalet, lon, lat, oro, lonoro, latoro, xlim, ylim, lapse=lapse)$output1 - names(dim(result))[1] <- names(dim(tas))[names(dim(tas))==lon_dim] - names(dim(result))[2] <- names(dim(tas))[names(dim(tas))==lat_dim] + names(dim(result))[1] <- names(dim(data))[names(dim(data))==lon_dim] + names(dim(result))[2] <- names(dim(data))[names(dim(data))==lat_dim] return(list(data = result, lon = lonoro, lat = latoro)) } @@ -103,9 +108,11 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, if (length(dim(t))==3) { nt = dim(t)[3] - } else { + } else if (length(dim(t))==2) { nt = 1 dim(t) <- c(dim(t), 1) + } else { + stop("Input array must be 2D or 3D") } dxl <- lon[2] - lon[1] @@ -123,8 +130,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, lonocut <- lono[(lono <= lon2) & (lono >= lon1)] latocut <- lato[(lato <= lat2) & (lato >= lat1)] - if (length(dim(tcut0))<3) { - dim(tcut0) <- c(dim(tcut0), 1) + if(length(lonocut)==0 | length(latocut)==0) { + stop("Orography not available for selected area") } nx <- length(lonocut) @@ -137,8 +144,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, jj <- ((latocut-(lat[1]-dy/2)) %/% dy + 1) ii <- ((lonocut-(lon[1]-dx/2)) %/% dx + 1) - if((ii[1]<0)|(ii[length(ii)]>length(lon))| - (jj[1]<0)|(jj[length(jj)]>length(lat))) { + 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") } for (k in seq_len(nt)) { -- GitLab From fc78f4c705ceed4cc48ef5a834505f9d5debdb37 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Wed, 26 Feb 2020 17:48:41 +0100 Subject: [PATCH 06/41] added examples --- R/CST_RFTemp.R | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 0bc7d631..95a46b67 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -21,7 +21,20 @@ #' @import multiApply #' @import rainfarmr #' @examples -#' +#' # Generate simple synthetic data and downscale by factor 4 +#' t <- rnorm(6 * 6 * 2 * 3 * 4)*10 + 273.15 + 10 +#' dim(t) <- c(dataset = 1, 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) +#' 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, ylim, lapse=6.5, lon_dim = "lon", lat_dim = "lat") { @@ -63,7 +76,16 @@ CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, lon_dim = "lon", lat_di #' @import multiApply #' @import rainfarmr #' @examples -#' +#' # Generate simple synthetic data and downscale by factor 4 +#' t <- rnorm(6 * 6 * 4)*10 + 273.15 + 10 +#' dim(t) <- c(ftime = 4, lat = 6, lon = 6) +#' lon <- seq(4, 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, ylim, lapse = 6.5, lon_dim = "lon", lat_dim = "lat") { -- GitLab From d1855836b3fe7178454220a9bff3c9983ff66864 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Sun, 1 Mar 2020 18:20:36 +0100 Subject: [PATCH 07/41] added delta method --- R/CST_RFTemp.R | 170 +++++++++++++++++++++++++++++++++++++------------ 1 file changed, 130 insertions(+), 40 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 95a46b67..e393e88a 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -1,5 +1,5 @@ #' @rdname CST_RFTemp -#' @title Temperature downscaling of a CSTools object using lapse rate correction +#' @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. @@ -16,6 +16,8 @@ #' @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 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 #' @return CST_RFTemp() returns a downscaled CSTools object (i.e., of the class 's2dv_cube'). #' @export #' @import multiApply @@ -36,7 +38,9 @@ #' 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, ylim, lapse=6.5, lon_dim = "lon", lat_dim = "lat") { +CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, + lon_dim = "lon", lat_dim = "lat", time_dim = NULL, + nolapse = FALSE, verbose=FALSE) { if (!inherits(data, 's2dv_cube')) { stop("Parameter 'data' must be of the class 's2dv_cube', ", @@ -44,8 +48,9 @@ CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, lon_dim = "lon", lat_di } 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) + oro$data, oro$lon, oro$lat, xlim, ylim, lapse, + lon_dim = lon_dim, lat_dim = lat_dim, + nolapse = nolapse, verbose = verbose) data$data <- res$data data$lon <- res$lon @@ -71,14 +76,16 @@ CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, lon_dim = "lon", lat_di #' @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 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 #' @return RFTemp() returns a list containing the fine-scale longitudes, latitudes and the downscaled fields. #' @export #' @import multiApply #' @import rainfarmr #' @examples #' # Generate simple synthetic data and downscale by factor 4 -#' t <- rnorm(6 * 6 * 4)*10 + 273.15 + 10 -#' dim(t) <- c(ftime = 4, lat = 6, lon = 6) +#' t <- rnorm(6 * 6 * 4 * 3)*10 + 273.15 + 10 +#' dim(t) <- c(sdate = 3, ftime = 4, lat = 6, lon = 6) #' lon <- seq(4, 9, 1) #' lat <- seq(42, 47, 1) #' o <- runif(29*29)*3000 @@ -88,11 +95,34 @@ CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, lon_dim = "lon", lat_di #' 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, ylim, lapse = 6.5, lon_dim = "lon", lat_dim = "lat") { + xlim, ylim, lapse = 6.5, lon_dim = "lon", lat_dim = "lat", + time_dim = NULL, nolapse = FALSE, verbose=FALSE) { + + # 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[1]]] + 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 - result <- Apply(data, target_dims=c(lon_dim, lat_dim), fun=.downscalet, - lon, lat, oro, lonoro, latoro, xlim, ylim, lapse=lapse)$output1 + result <- Apply(data, target_dims=c(lon_dim, lat_dim, time_dim), + fun=.downscalet, lon, lat, oro, lonoro, latoro, + xlim, ylim, lapse=lapse, + nolapse=nolapse, verbose=verbose)$output1 names(dim(result))[1] <- names(dim(data))[names(dim(data))==lon_dim] names(dim(result))[2] <- names(dim(data))[names(dim(data))==lat_dim] @@ -103,7 +133,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' Lapse-rate temperature correction downscaling #' #' @description Downscales a temperature field using a lapse-rate -#' correction based on a reference orography +#' 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 @@ -115,6 +146,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' @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 verbose logical if to print diagnostic output #' @return A downscaled temperature matrix #' @export #' @examples @@ -124,17 +157,24 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' lono=seq(5,20,0.1) #' lato=seq(35,40,0.1) #' o=runif(151*51)*2000; dim(o)=c(151,51) -#' td=downscalet(lon, lat,t,lono,lato,o,c(8,12),c(36,38)) - -.downscalet <- function(t, lon, lat, oro, lono, lato, xlim, ylim, radius = 0, lapse = 6.5) { +#' td=.downscalet(t,lon, lat,o,lono,lato,c(8,12),c(36,38)) - if (length(dim(t))==3) { - nt = dim(t)[3] - } else if (length(dim(t))==2) { +.downscalet <- function(t, lon, lat, oro, lono, lato, xlim, ylim, + radius = 0, lapse = 6.5, nolapse = FALSE, + verbose = FALSE) { + tdim <- dim(t) + ldim <- FALSE + if (length(tdim)==3) { + nt = tdim[3] + } else if (length(tdim)==2) { nt = 1 - dim(t) <- c(dim(t), 1) + dim(t) <- c(tdim, 1) + } else if (length(tdim)<2) { + stop("Input array must have at least two dimensions") } else { - stop("Input array must be 2D or 3D") + ldim <- TRUE + dim(t) <- c(tdim[1:2], time = prod(tdim[-(1:2)])) + nt = dim(t)[3] } dxl <- lon[2] - lon[1] @@ -155,40 +195,91 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, if(length(lonocut)==0 | length(latocut)==0) { stop("Orography not available for selected area") } - + + if(verbose) { + print("Nearest neighbours interpolation") + } + tcut <- .interpnn(lon, lat, t, lonocut, latocut) + # Fine-scale smooth interpolated input field + dxol <- lonocut[2] - lonocut[1] + nrad <- as.integer(radius / abs(dxol) + 0.5) nx <- length(lonocut) ny <- length(latocut) - + tout <- array(0., c(nx, ny, nt)) + for (i in seq_len(nt)) { + if(verbose) { + print(paste("Smooth ",i)) + } + tout[,, i] <- .smooth(tcut[,,i], nrad) + } + + 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 + ns <- length(lonocut) + orocuts <- .smooth(orocut, nrad) + t1 <- -(orocut-orocuts)*lapse/1000. + } + if(verbose) { + print("Adding") + } + for (i in seq_len(nt)) { + tout[,, i] <- tout[,, i] + t1 + } + + tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1]),] + if(ldim) { + dim(tdown) <- c(dim(tdown)[1:2], tdim[-(1:2)]) +   } + return(tdown) +} + +#' Nearest neighbour interpolation +#' +#' @description The input field is interpolated onto the output coordinate grid using nearest neighbours +#' @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 +#' @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 <- .interpnn(lon, lat, z, lonp, latp) + +.interpnn <- function(lon, lat, z, lonp, latp) { + + nx <- length(lonp) + ny <- length(latp) + nt <- dim(z)[3] + # Interpolate nn assuming a regular grid - tcut <- array(0., c(nx, ny, nt)) + zo <- array(0., c(nx, ny, nt)) dy <- lat[2]-lat[1] dx <- lon[2]-lon[1] - jj <- ((latocut-(lat[1]-dy/2)) %/% dy + 1) - ii <- ((lonocut-(lon[1]-dx/2)) %/% dx + 1) + jj <- ((latp-(lat[1]-dy/2)) %/% dy + 1) + ii <- ((lonp-(lon[1]-dx/2)) %/% dx + 1) if((ii[1]<=0)|(ii[length(ii)]>length(lon))| - (jj[1]<=0)|(jj[length(jj)]>length(lat))) { + (jj[1]<=0)|(jj[length(jj)]>length(lat))) { stop("Downscaling area not contained in input data") } for (k in seq_len(nt)) { for (j in seq_along(jj)) { - tcut[,j,k] <- t[ii,jj[j],k] + zo[,j,k] <- z[ii,jj[j],k] } } - - dxol <- lonocut[2] - lonocut[1] - nrad <- as.integer(radius / abs(dxol) + 0.5) - - ns = length(lonocut) - orocuts=.smooth(orocut, nrad) - t1 <- -(orocut-orocuts)*lapse/1000. - - tout <- array(0., c(nx, ny, nt)) - for (i in seq_len(nt)) { - tout[,, i] <- .smooth(tcut[,,i], nrad) + t1 - } - downscalet <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), - (latocut <= ylim[2]) & (latocut >= ylim[1]),] + names(dim(zo)) <- names(dim(z)) + return(zo) } #' Smoothening using convolution with a circular kernel @@ -212,7 +303,6 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, z[imask] <- 0. nsx <- dim(z)[1] nsy <- dim(z)[2] - mask <- matrix(0., nsx, nsy) for (i in 1:nsx) { for (j in 1:nsy) { -- GitLab From e51dc6dcf2c4a710641118287c0666302368c446 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Sun, 1 Mar 2020 23:16:26 +0100 Subject: [PATCH 08/41] output correct lon and lat and dimnames --- R/CST_RFTemp.R | 50 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 35 insertions(+), 15 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index e393e88a..918ff5f4 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -38,7 +38,7 @@ #' 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, ylim, 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) { @@ -95,8 +95,9 @@ CST_RFTemp <- function(data, oro, xlim, ylim, lapse=6.5, #' 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, ylim, lapse = 6.5, lon_dim = "lon", lat_dim = "lat", - time_dim = NULL, nolapse = FALSE, verbose=FALSE) { + xlim = NULL, ylim = NULL, lapse = 6.5, + lon_dim = "lon", lat_dim = "lat", time_dim = NULL, + nolapse = FALSE, verbose=FALSE) { # Check/detect time_dim if (is.null(time_dim)) { @@ -121,12 +122,16 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, # Repeatedly apply .downscalet result <- Apply(data, target_dims=c(lon_dim, lat_dim, time_dim), fun=.downscalet, lon, lat, oro, lonoro, latoro, - xlim, ylim, lapse=lapse, - nolapse=nolapse, verbose=verbose)$output1 - names(dim(result))[1] <- names(dim(data))[names(dim(data))==lon_dim] - names(dim(result))[2] <- names(dim(data))[names(dim(data))==lat_dim] + xlim=xlim, ylim=ylim, lapse=lapse, + nolapse=nolapse, verbose=verbose) - return(list(data = result, lon = lonoro, lat = latoro)) + 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) } @@ -159,7 +164,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' o=runif(151*51)*2000; dim(o)=c(151,51) #' td=.downscalet(t,lon, lat,o,lono,lato,c(8,12),c(36,38)) -.downscalet <- function(t, lon, lat, oro, lono, lato, xlim, ylim, +.downscalet <- function(t, lon, lat, oro, lono, lato, + xlim = NULL, ylim = NULL, radius = 0, lapse = 6.5, nolapse = FALSE, verbose = FALSE) { tdim <- dim(t) @@ -177,16 +183,27 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, nt = dim(t)[3] } - dxl <- lon[2] - lon[1] + dxl <- (lon[2] - lon[1]) / 2 + dyl <- (lat[2] - lat[1]) / 2 if(radius==0) { - radius <- dxl/2 + 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)] @@ -229,13 +246,16 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, for (i in seq_len(nt)) { tout[,, i] <- tout[,, i] + t1 } - tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), - (latocut <= ylim[2]) & (latocut >= ylim[1]),] + (latocut <= ylim[2]) & (latocut >= ylim[1]),] if(ldim) { dim(tdown) <- c(dim(tdown)[1:2], tdim[-(1:2)])   } - return(tdown) + lonocut <- lonocut[(lonocut <= xlim[2]) & (lonocut >= xlim[1])] + latocut <- latocut[(latocut <= ylim[2]) & (latocut >= ylim[1])] + names(dim(lonocut))="lon" + names(dim(latocut))="lat" + return(list(data=tdown, lon=lonocut, lat=latocut)) } #' Nearest neighbour interpolation -- GitLab From f9c94a611eef4f3ba97356b6bc4ed64057ce0beb Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Mon, 2 Mar 2020 00:39:08 +0100 Subject: [PATCH 09/41] minor fixes --- R/CST_RFTemp.R | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 918ff5f4..a1fe30ec 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -49,7 +49,7 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, 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, + lon_dim = lon_dim, lat_dim = lat_dim, time_dim = time_dim, nolapse = nolapse, verbose = verbose) data$data <- res$data @@ -197,7 +197,6 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, ylim <- ylim[2:1] } } - #Add buffer lon1 <- xlim[1] - radius lon2 <- xlim[2] + radius @@ -253,8 +252,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro,   } lonocut <- lonocut[(lonocut <= xlim[2]) & (lonocut >= xlim[1])] latocut <- latocut[(latocut <= ylim[2]) & (latocut >= ylim[1])] - names(dim(lonocut))="lon" - names(dim(latocut))="lat" + return(list(data=tdown, lon=lonocut, lat=latocut)) } -- GitLab From fbe4599a114b91910850a024f4335609db1a4b98 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Mon, 2 Mar 2020 00:39:34 +0100 Subject: [PATCH 10/41] additional tests --- tests/testthat/test-CST_RFTemp.R | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/tests/testthat/test-CST_RFTemp.R b/tests/testthat/test-CST_RFTemp.R index 9d577d64..346b66e2 100644 --- a/tests/testthat/test-CST_RFTemp.R +++ b/tests/testthat/test-CST_RFTemp.R @@ -8,8 +8,8 @@ test_that("Sanity checks and simple use cases", { exp <- list(data = t, lat = lat, lon = lon) o <- runif(29*29)*3000 dim(o) <- c(lat = 29, lon = 29) - lon <- seq(3, 10, 0.25) - 100 - lat <- seq(41, 48, 0.25) -60 + 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' @@ -23,6 +23,7 @@ test_that("Sanity checks and simple use cases", { 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 @@ -30,11 +31,23 @@ test_that("Sanity checks and simple use cases", { res <- CST_RFTemp(exp, oro, xlim=c(3,8), ylim=c(43, 46)), "Downscaling area not contained in input data" ) - - res <- CST_RFTemp(exp, oro, xlim=c(4,8), ylim=c(43, 46), lapse=6.5) + + expect_warning( + res <- CST_RFTemp(exp, oro, xlim=c(4,8), ylim=c(43, 46), lapse=6.5), + "Selected time dim: ftime" + ) dimexp=dim(exp$data) - expect_equal(dim(res$data), c(lon=17,lat=13, - dimexp["dataset"], dimexp["member"], - dimexp["sdate"], dimexp["ftime"])) + expect_equal(dim(res$data), c(lon=16,lat=12, dimexp["ftime"], + dimexp["sdate"], + dimexp["dataset"], dimexp["member"])) + expect_warning( + res <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), nolapse=TRUE), + "Selected time dim: sdate" + ) + dim(res$data) <- c(16,16,4*3*2) + resm <- apply(res$data, c(1,2), mean) + + expect_equal(mean(oro$data[7:22,7:22]), mean(resm), tolerance = 1e-10) + }) -- GitLab From c70cafcaa89a80990fb676a8f71af206282d95f7 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Mon, 2 Mar 2020 01:38:28 +0100 Subject: [PATCH 11/41] faster version --- R/CST_RFTemp.R | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index a1fe30ec..b723a0f0 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -221,13 +221,11 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, nrad <- as.integer(radius / abs(dxol) + 0.5) nx <- length(lonocut) ny <- length(latocut) - tout <- array(0., c(nx, ny, nt)) - for (i in seq_len(nt)) { - if(verbose) { - print(paste("Smooth ",i)) - } - tout[,, i] <- .smooth(tcut[,,i], nrad) + + if(verbose) { + print(paste("Smooth")) } + tout <- .smooth(tcut, nrad) if (nolapse) { # oro is a reference fine-scale field: use that directly @@ -321,6 +319,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, z[imask] <- 0. nsx <- dim(z)[1] nsy <- dim(z)[2] + nt <- dim(z)[3] mask <- matrix(0., nsx, nsy) for (i in 1:nsx) { for (j in 1:nsy) { @@ -333,17 +332,21 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, ky <- j - nsy - 1 } r2 <- kx * kx + ky * ky - #mask[i, j] <- exp( - (r2 / (sdim * sdim)) / 2) mask[i, j] <- (r2 <= (sdim * sdim)) } } fm <- fft(mask) - zf <- Re(fft(fm * fft(z), inverse = TRUE)) / sum(mask) / length(fm) - if (sum(imask) > 0) { - z[!imask] <- 1.0 - zf <- zf / (Re(fft(fm * fft(z), inverse = TRUE)) / + 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 + } } - zf[imask] <- NA return(zf) } -- GitLab From 170f2d1e1980714affa8ea277303e93f5f52ba1b Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Mon, 2 Mar 2020 02:02:04 +0100 Subject: [PATCH 12/41] hidden typo --- R/CST_RFTemp.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index b723a0f0..3bcf3f84 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -247,7 +247,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, (latocut <= ylim[2]) & (latocut >= ylim[1]),] if(ldim) { dim(tdown) <- c(dim(tdown)[1:2], tdim[-(1:2)]) -   } + } lonocut <- lonocut[(lonocut <= xlim[2]) & (lonocut >= xlim[1])] latocut <- latocut[(latocut <= ylim[2]) & (latocut >= ylim[1])] -- GitLab From 52d20121aff9055331a1e0f1fa250879fe625b59 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Mon, 2 Mar 2020 02:40:01 +0100 Subject: [PATCH 13/41] delinting --- R/CST_RFTemp.R | 178 +++++++++++++++++++++++++++---------------------- 1 file changed, 97 insertions(+), 81 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 3bcf3f84..c9306ba8 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -1,24 +1,32 @@ #' @rdname CST_RFTemp -#' @title Temperature downscaling of a CSTools object using lapse rate correction or a reference field +#' @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. -#' @references Method described in H2020 ECOPOTENTIAL Deliverable No. 8.1: High resolution (1-10 km) climate, land use and ocean change scenarios \link[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`, +#' temperature field (an object of class 's2dv_cube' as provided by +#' `CST_Load`) as input. +#' @references Method described in H2020 ECOPOTENTIAL Deliverable No. 8.1: +#' High resolution (1-10 km) climate, land use and ocean change scenarios +#' \link[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 -#' The number of longitudes and latitudes in the input data is expected to be even and the same. If not +#' The data object is expected to have an element named \code{$data} +#' with at least two +#' The number of longitudes and latitudes in the input data is expected to +#' be even and the same. If not #' spatial dimensions named "lon" and "lat". -#' @param oro An object of the class 's2dv_cube' as returned by `CST_Load`, +#' @param oro An object of the class 's2dv_cube' as returned by `CST_Load`, #' containing fine scale orography (in meters). #' @param xlim vector with longitude bounds for downscaling #' @param ylim vector with latitude bounds for downscaling -#' @param lapse float with environmental lapse rate +#' @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 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 -#' @return CST_RFTemp() returns a downscaled CSTools object (i.e., of the class 's2dv_cube'). +#' @param nolapse logical, if true `oro` is interpreted as a fine-scale +#' climatology and used directly for bias correction +#' @return CST_RFTemp() returns a downscaled CSTools object +#' (i.e., of the class 's2dv_cube'). #' @export #' @import multiApply #' @import rainfarmr @@ -42,7 +50,7 @@ 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) { - if (!inherits(data, 's2dv_cube')) { + if (!inherits(data, "s2dv_cube")) { stop("Parameter 'data' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") } @@ -60,25 +68,32 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, } #' @rdname RFTemp -#' @title Temperature downscaling of a CSTools object using lapse rate correction (reduced version) +#' @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 -#' @references Method described in H2020 ECOPOTENTIAL Deliverable No. 8.1: High resolution (1-10 km) climate, land use and ocean change scenarios \link[https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf] +#' @references Method described in H2020 ECOPOTENTIAL Deliverable No. 8.1: +#' High resolution (1-10 km) climate, land use and ocean change scenarios +#' \link[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) +#' 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 oro Array containing fine-scale orography (in m) #' @param xlim vector with longitude bounds for downscaling #' @param ylim vector with latitude bounds for downscaling -#' @param lapse float with environmental lapse rate +#' @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 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 -#' @return RFTemp() returns a list containing the fine-scale longitudes, latitudes and the downscaled fields. +#' @param nolapse logical, if true `oro` is interpreted as a +#' fine-scale climatology and used directly for bias correction +#' @return RFTemp() returns a list containing the fine-scale +#' longitudes, latitudes and the downscaled fields. #' @export #' @import multiApply #' @import rainfarmr @@ -92,7 +107,8 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' 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) +#' 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, @@ -107,7 +123,6 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, # 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[1]]] time_dim <- time_dim_names[time_dim_num[ilong]] } else { stop("No time dimension longer than zero found.") @@ -116,29 +131,28 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, stop("Could not automatically detect a target time dimension ", "in the provided data in 'data'.") } - warning(paste("Selected time dim:", time_dim,"\n")) + warning(paste("Selected time dim:", time_dim, "\n")) } # Repeatedly apply .downscalet - 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) + 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) - 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 + 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 +#' 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 @@ -146,12 +160,14 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' @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) +#' @param oro matrix of topographical elevations (in meters) #' @param xlim vector of longitude bounds #' @param ylim vector of latitude bounds -#' @param radius smoothing radius expressed in longitude units (default is half a large-scale pixel) +#' @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 nolapse logical, if true `oro` is interpreted as a fine-scale +#' climatology and used directly for bias correction #' @param verbose logical if to print diagnostic output #' @return A downscaled temperature matrix #' @export @@ -170,31 +186,31 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, verbose = FALSE) { tdim <- dim(t) ldim <- FALSE - if (length(tdim)==3) { - nt = tdim[3] - } else if (length(tdim)==2) { - nt = 1 + if (length(tdim) == 3) { + nt <- tdim[3] + } else if (length(tdim) == 2) { + nt <- 1 dim(t) <- c(tdim, 1) - } else if (length(tdim)<2) { + } 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] + dim(t) <- c(tdim[1:2], time = prod(tdim[-(1:2)])) + nt <- dim(t)[3] } dxl <- (lon[2] - lon[1]) / 2 dyl <- (lat[2] - lat[1]) / 2 - if(radius==0) { + if (radius == 0) { radius <- dxl } - if(is.null(xlim)) { + if (is.null(xlim)) { xlim <- c(lon[1] + dxl, lon[length(lon)] - dxl) } - if(is.null(ylim)) { + if (is.null(ylim)) { ylim <- c(lat[1] + dyl, lat[length(lat)] - dyl) if (ylim[1] > ylim[2]) { - ylim <- ylim[2:1] + ylim <- ylim[2:1] } } #Add buffer @@ -202,27 +218,25 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, 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)] - if(length(lonocut)==0 | length(latocut)==0) { + if (length(lonocut) == 0 | length(latocut) == 0) { stop("Orography not available for selected area") } - - if(verbose) { + + if (verbose) { print("Nearest neighbours interpolation") } tcut <- .interpnn(lon, lat, t, lonocut, latocut) # Fine-scale smooth interpolated input field dxol <- lonocut[2] - lonocut[1] nrad <- as.integer(radius / abs(dxol) + 0.5) - nx <- length(lonocut) - ny <- length(latocut) - if(verbose) { + if (verbose) { print(paste("Smooth")) } tout <- .smooth(tcut, nrad) @@ -230,35 +244,36 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, 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) + t1 <- orocut - apply(tout, c(1, 2), mean) } else { # Lapse-rate correction - ns <- length(lonocut) orocuts <- .smooth(orocut, nrad) - t1 <- -(orocut-orocuts)*lapse/1000. + t1 <- -(orocut - orocuts) * lapse / 1000. } - if(verbose) { + if (verbose) { print("Adding") } for (i in seq_len(nt)) { - tout[,, i] <- tout[,, i] + t1 + tout[, , i] <- tout[, , i] + t1 } tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), - (latocut <= ylim[2]) & (latocut >= ylim[1]),] - if(ldim) { + (latocut <= ylim[2]) & (latocut >= ylim[1]), ] + if (ldim) { dim(tdown) <- c(dim(tdown)[1:2], tdim[-(1:2)]) } lonocut <- lonocut[(lonocut <= xlim[2]) & (lonocut >= xlim[1])] latocut <- latocut[(latocut <= ylim[2]) & (latocut >= ylim[1])] - return(list(data=tdown, lon=lonocut, lat=latocut)) + return(list(data = tdown, lon = lonocut, lat = latocut)) } #' Nearest neighbour interpolation #' -#' @description The input field is interpolated onto the output coordinate grid using nearest neighbours +#' @description The input field is interpolated onto the output +#' coordinate grid using nearest neighbours #' @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 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 @@ -270,7 +285,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' z=runif(7*6*2); dim(z)=c(7,6,2) #' lonp=seq(5,10,0.2) #' latp=seq(35,40,0.2) -#' zo <- .interpnn(lon, lat, z, lonp, latp) +#' zo <- .interpnn(lon, lat, z, lonp, latp) .interpnn <- function(lon, lat, z, lonp, latp) { @@ -278,20 +293,20 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, ny <- length(latp) nt <- dim(z)[3] - # Interpolate nn assuming a regular grid + # Interpolate nn assuming a regular grid zo <- array(0., c(nx, ny, nt)) - dy <- lat[2]-lat[1] - dx <- lon[2]-lon[1] - jj <- ((latp-(lat[1]-dy/2)) %/% dy + 1) - ii <- ((lonp-(lon[1]-dx/2)) %/% dx + 1) + dy <- lat[2] - lat[1] + dx <- lon[2] - lon[1] + jj <- ((latp - (lat[1] - dy / 2)) %/% dy + 1) + ii <- ((lonp - (lon[1] - dx / 2)) %/% dx + 1) - if((ii[1]<=0)|(ii[length(ii)]>length(lon))| - (jj[1]<=0)|(jj[length(jj)]>length(lat))) { + 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") } for (k in seq_len(nt)) { for (j in seq_along(jj)) { - zo[,j,k] <- z[ii,jj[j],k] + zo[, j, k] <- z[ii, jj[j], k] } } names(dim(zo)) <- names(dim(z)) @@ -337,15 +352,16 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, } 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)) / + 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 + zz[imask[, , k]] <- NA + zf[, , k] <- zz } } return(zf) -- GitLab From becbc5cf21a35abd582cc282eac3126763b7d2b5 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Mon, 2 Mar 2020 11:53:43 +0100 Subject: [PATCH 14/41] smooth also 2d fields --- R/CST_RFTemp.R | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index c9306ba8..2bdbe192 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -330,11 +330,17 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' # [1] 0.1334648 .smooth <- function(z, sdim) { - imask <- !is.finite(z) - z[imask] <- 0. nsx <- dim(z)[1] nsy <- dim(z)[2] - nt <- dim(z)[3] + 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) { @@ -364,5 +370,6 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, zf[, , k] <- zz } } + dim(zf) <- zdim return(zf) } -- GitLab From bde12595b04ac3f01bcca673178df2f320a9c435 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Mon, 2 Mar 2020 12:33:11 +0100 Subject: [PATCH 15/41] fit test and delint --- tests/testthat/test-CST_RFTemp.R | 49 ++++++++++++++++---------------- 1 file changed, 24 insertions(+), 25 deletions(-) diff --git a/tests/testthat/test-CST_RFTemp.R b/tests/testthat/test-CST_RFTemp.R index 346b66e2..76ede1d0 100644 --- a/tests/testthat/test-CST_RFTemp.R +++ b/tests/testthat/test-CST_RFTemp.R @@ -1,53 +1,52 @@ context("Generic tests") test_that("Sanity checks and simple use cases", { # Generate simple synthetic data - t <- rnorm(6 * 6 * 2 * 3 * 4)*10 + 273.15 + 10 + t <- rnorm(6 * 6 * 2 * 3 * 4) * 10 + 273.15 + 10 dim(t) <- c(dataset = 1, 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 + 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 + lat <- seq(41.125, 48.125, 0.25) - 60 oro <- list(data = o, lat = lat, lon = lon) - attr(oro, 'class') <- 's2dv_cube' + attr(oro, "class") <- "s2dv_cube" expect_error( - res <- CST_RFTemp(exp, oro, xlim=c(1,3), ylim=c(1,3)), - "Parameter 'data' must be of the class 's2dv_cube', as output by CSTools::CST_Load." - ) - attr(exp, 'class') <- 's2dv_cube' + 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)), + 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 + + 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)), + res <- CST_RFTemp(exp, oro, xlim = c(3, 8), ylim = c(43, 46)), "Downscaling area not contained in input data" ) expect_warning( - res <- CST_RFTemp(exp, oro, xlim=c(4,8), ylim=c(43, 46), lapse=6.5), + res <- CST_RFTemp(exp, oro, xlim = c(4, 8), ylim = c(43, 46), lapse = 6.5), "Selected time dim: ftime" ) - dimexp=dim(exp$data) - expect_equal(dim(res$data), c(lon=16,lat=12, dimexp["ftime"], - dimexp["sdate"], - dimexp["dataset"], dimexp["member"])) + dimexp <- dim(exp$data) + expect_equal(dim(res$data), c(lon = 16, lat = 12, dimexp["ftime"], + dimexp["sdate"], dimexp["dataset"], + dimexp["member"])) expect_warning( - res <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), nolapse=TRUE), - "Selected time dim: sdate" - ) - dim(res$data) <- c(16,16,4*3*2) - resm <- apply(res$data, c(1,2), mean) + res <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), + nolapse = TRUE), NA) - expect_equal(mean(oro$data[7:22,7:22]), mean(resm), tolerance = 1e-10) + dim(res$data) <- c(16, 16, 4 * 3 * 2) + resm <- apply(res$data, c(1, 2), mean) -}) + expect_equal(mean(oro$data[7:22, 7:22]), mean(resm), tolerance = 1e-10) +}) -- GitLab From 39124dc9d92d6c0de52e82fc62582af821ce0d2e Mon Sep 17 00:00:00 2001 From: jhardenberg Date: Tue, 3 Mar 2020 00:29:06 +0100 Subject: [PATCH 16/41] regularize lon coords --- R/CST_RFTemp.R | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 2bdbe192..2b18ba3c 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -198,6 +198,11 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, dim(t) <- c(tdim[1:2], time = prod(tdim[-(1:2)])) nt <- dim(t)[3] } + # 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 -- GitLab From 942e3336fbee4a13bec32d9c7ef16d2dd9e6b496 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Tue, 3 Mar 2020 15:37:42 +0100 Subject: [PATCH 17/41] separate computation of delta --- R/CST_RFTemp.R | 133 +++++++++++++++++++++++++------ tests/testthat/test-CST_RFTemp.R | 28 +++++-- 2 files changed, 127 insertions(+), 34 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 2b18ba3c..853cbcf6 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -25,6 +25,11 @@ #' @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. #' @return CST_RFTemp() returns a downscaled CSTools object #' (i.e., of the class 's2dv_cube'). #' @export @@ -48,17 +53,29 @@ 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) { + nolapse = FALSE, verbose=FALSE, compute_delta = FALSE, + 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) + nolapse = nolapse, verbose = verbose, + compute_delta = compute_delta, delta = delta$data) data$data <- res$data data$lon <- res$lon @@ -92,6 +109,11 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' @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. Activates `nolapse = TRUE`. +#' @param delta matrix containing a delta to be applied to the downscaled input data. +#' Activates `nolapse = TRUE`. The grid of this matrix is supposed to be same as +#' that of the required output field #' @return RFTemp() returns a list containing the fine-scale #' longitudes, latitudes and the downscaled fields. #' @export @@ -113,7 +135,8 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, 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) { + nolapse = FALSE, verbose=FALSE, compute_delta = FALSE, + delta = NULL) { # Check/detect time_dim if (is.null(time_dim)) { @@ -135,10 +158,21 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, } # Repeatedly apply .downscalet - 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) + 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, + 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, + 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] @@ -168,6 +202,11 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' @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 #' @export @@ -183,7 +222,12 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, .downscalet <- function(t, lon, lat, oro, lono, lato, xlim = NULL, ylim = NULL, radius = 0, lapse = 6.5, nolapse = FALSE, - verbose = FALSE) { + verbose = FALSE, compute_delta = FALSE, 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) { @@ -246,25 +290,51 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, } tout <- .smooth(tcut, nrad) - 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) + 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 + tdown <- t1[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1])] + } else { + # Apply delta + if (verbose) { + print("Adding") + } + for (i in seq_len(nt)) { + tout[, , i] <- tout[, , i] + t1 + } + tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1]), ] + if (ldim) { + dim(tdown) <- c(dim(tdown)[1:2], tdim[-(1:2)]) + } + } } else { - # Lapse-rate correction - orocuts <- .smooth(orocut, nrad) - t1 <- -(orocut - orocuts) * lapse / 1000. - } - if (verbose) { - print("Adding") - } - for (i in seq_len(nt)) { - tout[, , i] <- tout[, , i] + t1 - } - tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), - (latocut <= ylim[2]) & (latocut >= ylim[1]), ] - if (ldim) { - dim(tdown) <- c(dim(tdown)[1:2], tdim[-(1:2)]) + # Just apply the provided delta + tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1]), ] + if (any(dim(delta)[1:2] != dim(tdown)[1:2])) { + stop("Provided delta has not the same dimensions as required output") + } + if (verbose) { + print("Adding") + } + for (i in seq_len(nt)) { + tdown[, , i] <- tdown[, , i] + delta + } + if (ldim) { + dim(tdown) <- c(dim(tdown)[1:2], tdim[-(1:2)]) + } } lonocut <- lonocut[(lonocut <= xlim[2]) & (lonocut >= xlim[1])] latocut <- latocut[(latocut <= ylim[2]) & (latocut >= ylim[1])] @@ -272,6 +342,17 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, return(list(data = tdown, 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) { + res <- .downscalet(t, lon, lat, oro, lono, lato, + xlim = NULL, ylim = NULL, radius = 0, + lapse = 6.5, nolapse = FALSE, verbose = FALSE, + compute_delta = FALSE, delta = delta) +} + #' Nearest neighbour interpolation #' #' @description The input field is interpolated onto the output diff --git a/tests/testthat/test-CST_RFTemp.R b/tests/testthat/test-CST_RFTemp.R index 76ede1d0..b9fbcf7f 100644 --- a/tests/testthat/test-CST_RFTemp.R +++ b/tests/testthat/test-CST_RFTemp.R @@ -1,8 +1,8 @@ context("Generic tests") test_that("Sanity checks and simple use cases", { # Generate simple synthetic data - t <- rnorm(6 * 6 * 2 * 3 * 4) * 10 + 273.15 + 10 - dim(t) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 6, lon = 6) + 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) @@ -33,20 +33,32 @@ test_that("Sanity checks and simple use cases", { ) expect_warning( - res <- CST_RFTemp(exp, oro, xlim = c(4, 8), ylim = c(43, 46), lapse = 6.5), + resl <- CST_RFTemp(exp, oro, lapse = 6.5), "Selected time dim: ftime" ) dimexp <- dim(exp$data) - expect_equal(dim(res$data), c(lon = 16, lat = 12, dimexp["ftime"], + expect_equal(dim(resl$data), c(lon = 16, lat = 16, dimexp["ftime"], dimexp["sdate"], dimexp["dataset"], dimexp["member"])) expect_warning( - res <- CST_RFTemp(exp, oro, time_dim = c("ftime", "sdate", "member"), + 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) - dim(res$data) <- c(16, 16, 4 * 3 * 2) - resm <- apply(res$data, c(1, 2), mean) + # 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) - expect_equal(mean(oro$data[7:22, 7:22]), mean(resm), 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) }) -- GitLab From 0d18b6b6d8a1c4c507374a94fe72e0ab34c09132 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Tue, 3 Mar 2020 15:59:04 +0100 Subject: [PATCH 18/41] delinting [skip ci] --- R/CST_RFTemp.R | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 853cbcf6..fc40ea6d 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -26,7 +26,8 @@ #' @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`. +#' 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. @@ -110,10 +111,10 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' @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. Activates `nolapse = TRUE`. -#' @param delta matrix containing a delta to be applied to the downscaled input data. -#' Activates `nolapse = TRUE`. The grid of this matrix is supposed to be same as -#' that of the required output field +#' 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 #' @return RFTemp() returns a list containing the fine-scale #' longitudes, latitudes and the downscaled fields. #' @export @@ -323,7 +324,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, # Just apply the provided delta tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), (latocut <= ylim[2]) & (latocut >= ylim[1]), ] - if (any(dim(delta)[1:2] != dim(tdown)[1:2])) { + if (any(dim(delta)[1:2] != dim(tdown)[1:2])) { stop("Provided delta has not the same dimensions as required output") } if (verbose) { @@ -420,7 +421,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, nsy <- dim(z)[2] zdim <- dim(z) if (length(dim(z)) == 2) { - dim(z) <- c(dim(z), time=1) + dim(z) <- c(dim(z), time = 1) nt <- 1 } else { nt <- dim(z)[3] -- GitLab From 73a2a9391912f36ebfdb051ba63c0e70b96935d2 Mon Sep 17 00:00:00 2001 From: jhardenberg Date: Wed, 4 Mar 2020 11:10:53 +0100 Subject: [PATCH 19/41] fixed applied delta case --- R/CST_RFTemp.R | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index fc40ea6d..0df18616 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -287,7 +287,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, nrad <- as.integer(radius / abs(dxol) + 0.5) if (verbose) { - print(paste("Smooth")) + print(paste("Smoothing interpolated field")) } tout <- .smooth(tcut, nrad) @@ -309,7 +309,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, } else { # Apply delta if (verbose) { - print("Adding") + print("Adding computed delta") } for (i in seq_len(nt)) { tout[, , i] <- tout[, , i] + t1 @@ -325,10 +325,12 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), (latocut <= ylim[2]) & (latocut >= ylim[1]), ] if (any(dim(delta)[1:2] != dim(tdown)[1:2])) { + print(dim(delta)) + print(dim(tdown)) stop("Provided delta has not the same dimensions as required output") } if (verbose) { - print("Adding") + print("Adding provided delta") } for (i in seq_len(nt)) { tdown[, , i] <- tdown[, , i] + delta @@ -349,9 +351,9 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, radius = 0, lapse = 6.5, nolapse = FALSE, verbose = FALSE, compute_delta = FALSE) { res <- .downscalet(t, lon, lat, oro, lono, lato, - xlim = NULL, ylim = NULL, radius = 0, - lapse = 6.5, nolapse = FALSE, verbose = FALSE, - compute_delta = FALSE, delta = delta) + xlim = xlim, ylim = ylim, radius = radius, + lapse = lapse, nolapse = nolapse, verbose = verbose, + compute_delta = compute_delta, delta = delta) } #' Nearest neighbour interpolation -- GitLab From 0c085d9eb75eed5000b8d5148eab298295d917fc Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Sat, 7 Mar 2020 11:22:24 +0100 Subject: [PATCH 20/41] faster interpnn --- R/CST_RFTemp.R | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 0df18616..37431b04 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -393,11 +393,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, (jj[1] <= 0) | (jj[length(jj)] > length(lat))) { stop("Downscaling area not contained in input data") } - for (k in seq_len(nt)) { - for (j in seq_along(jj)) { - zo[, j, k] <- z[ii, jj[j], k] - } - } + + zo[, , ] <- z[ii, jj, ] names(dim(zo)) <- names(dim(z)) return(zo) } -- GitLab From 522845f5c54480a522000c2a33b04ab904c70bb5 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Sat, 7 Mar 2020 18:17:53 +0100 Subject: [PATCH 21/41] added bilinear interp as default --- R/CST_RFTemp.R | 73 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 52 insertions(+), 21 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 37431b04..eb546188 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -31,6 +31,9 @@ #' @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) #' @return CST_RFTemp() returns a downscaled CSTools object #' (i.e., of the class 's2dv_cube'). #' @export @@ -55,7 +58,7 @@ 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, - delta = NULL) { + method = "bilinear", delta = NULL) { if (!inherits(data, "s2dv_cube")) { stop("Parameter 'data' must be of the class 's2dv_cube', ", @@ -75,7 +78,7 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, 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, + nolapse = nolapse, verbose = verbose, method = method, compute_delta = compute_delta, delta = delta$data) data$data <- res$data @@ -115,6 +118,9 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' @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) #' @return RFTemp() returns a list containing the fine-scale #' longitudes, latitudes and the downscaled fields. #' @export @@ -137,7 +143,7 @@ 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, - delta = NULL) { + method = "bilinear", delta = NULL) { # Check/detect time_dim if (is.null(time_dim)) { @@ -163,7 +169,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, 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, + nolapse = nolapse, verbose = verbose, method = method, compute_delta = compute_delta) } else { result <- Apply(list(data, delta), @@ -171,7 +177,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, c(lon_dim, lat_dim)), fun = .downscalet_delta, lon, lat, oro, lonoro, latoro, xlim = xlim, ylim = ylim, lapse = lapse, - nolapse = nolapse, verbose = verbose, + nolapse = nolapse, verbose = verbose, method = method, compute_delta = compute_delta) } @@ -223,7 +229,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, .downscalet <- function(t, lon, lat, oro, lono, lato, xlim = NULL, ylim = NULL, radius = 0, lapse = 6.5, nolapse = FALSE, - verbose = FALSE, compute_delta = FALSE, delta = NULL) { + 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.") @@ -279,16 +286,15 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, } if (verbose) { - print("Nearest neighbours interpolation") + print(paste("Interpolation using", method, "method")) } - tcut <- .interpnn(lon, lat, t, lonocut, latocut) + tcut <- .interp2d(lon, lat, t, lonocut, latocut, method = method) # Fine-scale smooth interpolated input field - dxol <- lonocut[2] - lonocut[1] - nrad <- as.integer(radius / abs(dxol) + 0.5) - if (verbose) { print(paste("Smoothing interpolated field")) } + dxol <- lonocut[2] - lonocut[1] + nrad <- as.integer(radius / abs(dxol) + 0.5) tout <- .smooth(tcut, nrad) if (is.null(delta)) { @@ -359,7 +365,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' Nearest neighbour interpolation #' #' @description The input field is interpolated onto the output -#' coordinate grid using nearest neighbours +#' coordinate grid using nearest neighbours or bilinear interpolation #' @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) @@ -367,6 +373,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' @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 @@ -374,9 +382,9 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' z=runif(7*6*2); dim(z)=c(7,6,2) #' lonp=seq(5,10,0.2) #' latp=seq(35,40,0.2) -#' zo <- .interpnn(lon, lat, z, lonp, latp) +#' zo <- .interp2d(lon, lat, z, lonp, latp, method="nearest") -.interpnn <- function(lon, lat, z, lonp, latp) { +.interp2d <- function(lon, lat, z, lonp, latp, method="bilinear") { nx <- length(lonp) ny <- length(latp) @@ -386,15 +394,38 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, zo <- array(0., c(nx, ny, nt)) dy <- lat[2] - lat[1] dx <- lon[2] - lon[1] - jj <- ((latp - (lat[1] - dy / 2)) %/% dy + 1) - ii <- ((lonp - (lon[1] - dx / 2)) %/% dx + 1) - 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") - } + if (method == "nearest") { + jj <- ((latp - (lat[1] - dy / 2)) %/% dy + 1) + ii <- ((lonp - (lon[1] - dx / 2)) %/% dx + 1) - zo[, , ] <- z[ii, jj, ] + 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 <- ((latp - lat[1]) %/% dy + 1) + ii <- ((lonp - lon[1]) %/% dx + 1) + + 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") + } + xx <- (lonp - lon[ii]) / dx + yy <- (latp - lat[jj]) / dy + xx <- xx[row(zo[, , 1])] + yy <- yy[col(zo[, , 1])] + dim(xx) <- c(nx, ny) + dim(yy) <- c(nx, ny) + + for (k in seq_len(nt)) { + zo[, , k] <- z[ii, jj, k] * (1 - xx) * (1 - yy) + + z[ii + 1, jj, k] * xx * (1 - yy) + + z[ii, jj + 1, k] * (1 - xx) * yy + + z[ii + 1, jj + 1, k] * xx * yy + } + } names(dim(zo)) <- names(dim(z)) return(zo) } -- GitLab From d0c52753c88e7057d2ca5e177eb7d192bdd2758a Mon Sep 17 00:00:00 2001 From: jhardenberg Date: Sat, 7 Mar 2020 23:33:45 +0100 Subject: [PATCH 22/41] small fix and added test --- R/CST_RFTemp.R | 36 +++++++++++++++++++------------- tests/testthat/test-CST_RFTemp.R | 7 +++++++ 2 files changed, 29 insertions(+), 14 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index eb546188..de421c12 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -34,11 +34,13 @@ #' @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 -#' @import rainfarmr #' @examples #' # Generate simple synthetic data and downscale by factor 4 #' t <- rnorm(6 * 6 * 2 * 3 * 4)*10 + 273.15 + 10 @@ -121,11 +123,14 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' @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 -#' @import rainfarmr #' @examples #' # Generate simple synthetic data and downscale by factor 4 #' t <- rnorm(6 * 6 * 4 * 3)*10 + 273.15 + 10 @@ -288,14 +293,16 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, if (verbose) { print(paste("Interpolation using", method, "method")) } - tcut <- .interp2d(lon, lat, t, lonocut, latocut, method = method) + tout <- .interp2d(lon, lat, t, lonocut, latocut, method = method) # Fine-scale smooth interpolated input field - if (verbose) { - print(paste("Smoothing interpolated field")) - } dxol <- lonocut[2] - lonocut[1] nrad <- as.integer(radius / abs(dxol) + 0.5) - tout <- .smooth(tcut, nrad) + if (method == "nearest") { + if (verbose) { + print(paste("Smoothing interpolated field")) + } + tout <- .smooth(tout, nrad) + } if (is.null(delta)) { # Compute delta @@ -353,13 +360,14 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, # 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) { - 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) + 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 diff --git a/tests/testthat/test-CST_RFTemp.R b/tests/testthat/test-CST_RFTemp.R index b9fbcf7f..536e44f0 100644 --- a/tests/testthat/test-CST_RFTemp.R +++ b/tests/testthat/test-CST_RFTemp.R @@ -61,4 +61,11 @@ test_that("Sanity checks and simple use cases", { 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) }) -- GitLab From 4f13f116c8ef0edc9b9bebc4f412497dc78c85bf Mon Sep 17 00:00:00 2001 From: jhardenberg Date: Tue, 17 Mar 2020 14:55:40 +0100 Subject: [PATCH 23/41] use less memory --- R/CST_RFTemp.R | 60 +++++++++++++++++++++++++++++--------------------- 1 file changed, 35 insertions(+), 25 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index de421c12..9162442b 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -286,24 +286,37 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, 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 (verbose) { - print(paste("Interpolation using", method, "method")) - } - tout <- .interp2d(lon, lat, t, lonocut, latocut, method = method) - # Fine-scale smooth interpolated input field - dxol <- lonocut[2] - lonocut[1] - nrad <- as.integer(radius / abs(dxol) + 0.5) - if (method == "nearest") { + if (is.null(delta) & compute_delta & nolapse) { if (verbose) { - print(paste("Smoothing interpolated field")) + print("Time averaging") } - tout <- .smooth(tout, nrad) + # 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) { @@ -317,8 +330,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, } if (compute_delta) { # Return delta - tdown <- t1[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), - (latocut <= ylim[2]) & (latocut >= ylim[1])] + tout <- t1[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1])] } else { # Apply delta if (verbose) { @@ -327,35 +340,33 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, for (i in seq_len(nt)) { tout[, , i] <- tout[, , i] + t1 } - tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), - (latocut <= ylim[2]) & (latocut >= ylim[1]), ] + tout <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), + (latocut <= ylim[2]) & (latocut >= ylim[1]), ] if (ldim) { - dim(tdown) <- c(dim(tdown)[1:2], tdim[-(1:2)]) + dim(tout) <- c(dim(tout)[1:2], tdim[-(1:2)]) } } } else { # Just apply the provided delta - tdown <- tout[(lonocut <= xlim[2]) & (lonocut >= xlim[1]), - (latocut <= ylim[2]) & (latocut >= ylim[1]), ] - if (any(dim(delta)[1:2] != dim(tdown)[1:2])) { - print(dim(delta)) - print(dim(tdown)) + 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)) { - tdown[, , i] <- tdown[, , i] + delta + tout[, , i] <- tout[, , i] + delta } if (ldim) { - dim(tdown) <- c(dim(tdown)[1:2], tdim[-(1:2)]) + 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 = tdown, lon = lonocut, lat = latocut)) + return(list(data = tout, lon = lonocut, lat = latocut)) } # Special driver version of .downscalet to apply delta @@ -392,12 +403,11 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' latp=seq(35,40,0.2) #' zo <- .interp2d(lon, lat, z, lonp, latp, method="nearest") -.interp2d <- function(lon, lat, z, lonp, latp, method="bilinear") { +.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] -- GitLab From 889605e83d4d184f4ac1c741f542d2b44373d191 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Sat, 13 Jun 2020 18:21:18 +0200 Subject: [PATCH 24/41] add reference to MEDSCOPE M3.2 --- R/CST_RFTemp.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 9162442b..6e865031 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -5,7 +5,10 @@ #' @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. -#' @references Method described in H2020 ECOPOTENTIAL Deliverable No. 8.1: +#' @references Method described in ERA4CS MEDSCOPE milestone M3.2: +#' High-quality climate prediction data available to WP4 +#' \link[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 #' \link[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`, @@ -96,7 +99,10 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' @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 -#' @references Method described in H2020 ECOPOTENTIAL Deliverable No. 8.1: +#' @references Method described in ERA4CS MEDSCOPE milestone M3.2: +#' High-quality climate prediction data available to WP4 +#' \link[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 #' \link[https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf] #' @param data Temperature array to downscale. -- GitLab From 7dcbc516da6dceebd1cd00e7bc4f4f1bc2a03124 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Sat, 13 Jun 2020 18:28:40 +0200 Subject: [PATCH 25/41] doc update --- R/CST_RFTemp.R | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 6e865031..d9cf041b 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -14,13 +14,14 @@ #' @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 -#' The number of longitudes and latitudes in the input data is expected to -#' be even and the same. If not -#' spatial dimensions named "lon" and "lat". +#' 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). -#' @param xlim vector with longitude bounds for downscaling +#' 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 @@ -113,7 +114,9 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' @param lon Vector or array of longitudes. #' @param lat Vector or array of latitudes. #' @param oro Array containing fine-scale orography (in m) -#' @param xlim vector with longitude bounds for downscaling +#' 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 @@ -213,7 +216,8 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' @param lono vector of orography longitudes #' @param lato vector of orography latitudes #' @param oro matrix of topographical elevations (in meters) -#' @param xlim vector of longitude bounds +#' 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) -- GitLab From 81132b4a0f8a161e8462a177b1379c75d88a7e67 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Sun, 14 Jun 2020 15:06:17 +0200 Subject: [PATCH 26/41] fix lon=0deg issue --- R/CST_RFTemp.R | 39 ++++++++++++++++++++++++++++++--------- 1 file changed, 30 insertions(+), 9 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index d9cf041b..066b2726 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -5,6 +5,7 @@ #' @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 and output grids must be regularly spaced. #' @references Method described in ERA4CS MEDSCOPE milestone M3.2: #' High-quality climate prediction data available to WP4 #' \link[https://www.medscope-project.eu/the-project/deliverables-reports/] @@ -99,7 +100,8 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' 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 +#' temperature field (a multidimensional array) as input. +#' The input and output grids must be regularly spaced. #' @references Method described in ERA4CS MEDSCOPE milestone M3.2: #' High-quality climate prediction data available to WP4 #' \link[https://www.medscope-project.eu/the-project/deliverables-reports/] @@ -265,6 +267,9 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, 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 } @@ -394,7 +399,9 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' Nearest neighbour interpolation #' #' @description The input field is interpolated onto the output -#' coordinate grid using nearest neighbours or bilinear interpolation +#' coordinate grid using nearest neighbours or bilinear interpolation. +#' The input and output grids must be regularly spaced and coordinates +#' must be monotone increasing. #' @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) @@ -426,7 +433,11 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, if (method == "nearest") { jj <- ((latp - (lat[1] - dy / 2)) %/% dy + 1) 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") @@ -435,12 +446,22 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, } else { jj <- ((latp - lat[1]) %/% dy + 1) ii <- ((lonp - lon[1]) %/% dx + 1) - + ii2 <- ii + 1 + jj2 <- jj + 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))) { + (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 + xx <- ((lonp - lon[ii]) / dx) %% 360 yy <- (latp - lat[jj]) / dy xx <- xx[row(zo[, , 1])] yy <- yy[col(zo[, , 1])] @@ -449,9 +470,9 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, for (k in seq_len(nt)) { zo[, , k] <- z[ii, jj, k] * (1 - xx) * (1 - yy) + - z[ii + 1, jj, k] * xx * (1 - yy) + - z[ii, jj + 1, k] * (1 - xx) * yy + - z[ii + 1, jj + 1, k] * xx * yy + z[ii2, jj, k] * xx * (1 - yy) + + z[ii, jj2, k] * (1 - xx) * yy + + z[ii2, jj2, k] * xx * yy } } names(dim(zo)) <- names(dim(z)) -- GitLab From c65a6e8ca21c0b56828ee532c8b2f6474a542ee6 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Sun, 14 Jun 2020 16:21:04 +0200 Subject: [PATCH 27/41] correct domain limits calculation --- R/CST_RFTemp.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 066b2726..fa281df0 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -455,10 +455,10 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, 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))) { + 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 -- GitLab From f096c4a3012a3cc1ee5f31fb6a9ee60c0349dda2 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Mon, 15 Jun 2020 20:00:12 +0200 Subject: [PATCH 28/41] fix example --- R/CST_RFTemp.R | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index fa281df0..f2b9112c 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -48,9 +48,9 @@ #' @import multiApply #' @examples #' # Generate simple synthetic data and downscale by factor 4 -#' t <- rnorm(6 * 6 * 2 * 3 * 4)*10 + 273.15 + 10 -#' dim(t) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 6, lon = 6) -#' lon <- seq(4, 9, 1) +#' 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' @@ -144,15 +144,15 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' @import multiApply #' @examples #' # Generate simple synthetic data and downscale by factor 4 -#' t <- rnorm(6 * 6 * 4 * 3)*10 + 273.15 + 10 -#' dim(t) <- c(sdate = 3, ftime = 4, lat = 6, lon = 6) -#' lon <- seq(4, 9, 1) +#' 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), +#' 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, -- GitLab From 4d5a7324bc84e237d5f9a03aba85e9f586f3a6f6 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Wed, 17 Jun 2020 13:08:44 +0200 Subject: [PATCH 29/41] allow irregular lat grids --- R/CST_RFTemp.R | 46 +++++++++++++++++++++++++++++++++------------- 1 file changed, 33 insertions(+), 13 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index f2b9112c..9a22f964 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -5,7 +5,9 @@ #' @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 and output grids must be regularly spaced. +#' 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 #' \link[https://www.medscope-project.eu/the-project/deliverables-reports/] @@ -101,7 +103,9 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' @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 and output grids must be regularly spaced. +#' 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 #' \link[https://www.medscope-project.eu/the-project/deliverables-reports/] @@ -400,8 +404,9 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' #' @description The input field is interpolated onto the output #' coordinate grid using nearest neighbours or bilinear interpolation. -#' The input and output grids must be regularly spaced and coordinates -#' must be monotone increasing. +#' 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) @@ -431,7 +436,10 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, dx <- lon[2] - lon[1] if (method == "nearest") { - jj <- ((latp - (lat[1] - dy / 2)) %/% dy + 1) + 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) { @@ -444,10 +452,21 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, } zo[, , ] <- z[ii, jj, ] } else { - jj <- ((latp - lat[1]) %/% dy + 1) + 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 - jj2 <- jj + 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) @@ -462,17 +481,18 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, stop("Downscaling area not contained in input data") } xx <- ((lonp - lon[ii]) / dx) %% 360 - yy <- (latp - lat[jj]) / dy + 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] * (1 - xx) * (1 - yy) + - z[ii2, jj, k] * xx * (1 - yy) + - z[ii, jj2, k] * (1 - xx) * yy + - z[ii2, jj2, k] * xx * yy + 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)) -- GitLab From e8e8baa879620ef74473f3bf8002c0f4cfe33361 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Thu, 18 Jun 2020 15:00:37 +0200 Subject: [PATCH 30/41] update interp2d example --- R/CST_RFTemp.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 9a22f964..e111e313 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -423,7 +423,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' 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(lon, lat, z, lonp, latp, method="nearest") +#' zo <- .interp2d(z, lon, lat, lonp, latp, method="nearest") .interp2d <- function(z, lon, lat, lonp, latp, method="bilinear") { -- GitLab From 4c86acd4df17797f3c1426786852fa35e86c6c67 Mon Sep 17 00:00:00 2001 From: Jost von Hardenberg Date: Thu, 18 Jun 2020 17:07:48 +0200 Subject: [PATCH 31/41] removed exports from .internal functions --- R/CST_RFTemp.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index e111e313..a974b65f 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -237,7 +237,6 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' that of the required output field #' @param verbose logical if to print diagnostic output #' @return A downscaled temperature matrix -#' @export #' @examples #' lon=5:20 #' lat=35:40 @@ -507,7 +506,6 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' @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. -#' @export #' @examples #' z <- rnorm(64 * 64) #' dim(z) <- c(64, 64) -- GitLab From 68677fbca96fb7e7943db9797405aa023eafe31b Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 17:42:50 +0200 Subject: [PATCH 32/41] documentation generated with devtools --- NAMESPACE | 2 + man/CST_RFTemp.Rd | 105 ++++++++++++++++++++++++++++++++++++++++ man/RFTemp.Rd | 108 ++++++++++++++++++++++++++++++++++++++++++ man/dot-downscalet.Rd | 79 ++++++++++++++++++++++++++++++ man/dot-interp2d.Rd | 44 +++++++++++++++++ man/dot-smooth.Rd | 30 ++++++++++++ 6 files changed, 368 insertions(+) create mode 100644 man/CST_RFTemp.Rd create mode 100644 man/RFTemp.Rd create mode 100644 man/dot-downscalet.Rd create mode 100644 man/dot-interp2d.Rd create mode 100644 man/dot-smooth.Rd diff --git a/NAMESPACE b/NAMESPACE index 696a3e89..342847f3 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/man/CST_RFTemp.Rd b/man/CST_RFTemp.Rd new file mode 100644 index 00000000..a7beced6 --- /dev/null +++ b/man/CST_RFTemp.Rd @@ -0,0 +1,105 @@ +% 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{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 +\link[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 +\link[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 00000000..3945f571 --- /dev/null +++ b/man/RFTemp.Rd @@ -0,0 +1,108 @@ +% 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{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{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 +\link[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 +\link[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/dot-downscalet.Rd b/man/dot-downscalet.Rd new file mode 100644 index 00000000..9b28e05e --- /dev/null +++ b/man/dot-downscalet.Rd @@ -0,0 +1,79 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_RFTemp.R +\name{.downscalet} +\alias{.downscalet} +\title{Lapse-rate temperature correction downscaling} +\usage{ +.downscalet( + 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 +) +} +\arguments{ +\item{t}{matrix of input temperature data} + +\item{lon}{vector of input longitudes} + +\item{lat}{vector of input latitudes} + +\item{oro}{matrix of topographical elevations (in meters) +The destination downscaling area must be contained in the orography field.} + +\item{lono}{vector of orography longitudes} + +\item{lato}{vector of orography latitudes} + +\item{xlim}{vector of longitude bounds; the full input field is downscaled if `xlim` and `ylim` are not specified.} + +\item{ylim}{vector of latitude bounds} + +\item{radius}{smoothing radius expressed in longitude units +(default is half a large-scale pixel)} + +\item{lapse}{environmental lapse rate (in K/Km)} + +\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{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} +} +\value{ +A downscaled temperature matrix +} +\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. +} +\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)) +} +\author{ +Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/man/dot-interp2d.Rd b/man/dot-interp2d.Rd new file mode 100644 index 00000000..b03c6eee --- /dev/null +++ b/man/dot-interp2d.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_RFTemp.R +\name{.interp2d} +\alias{.interp2d} +\title{Nearest neighbour interpolation} +\usage{ +.interp2d(z, lon, lat, lonp, latp, method = "bilinear") +} +\arguments{ +\item{z}{matrix with the input field to interpolate (assumed to +include also a third time dimension)} + +\item{lon}{vector of input longitudes} + +\item{lat}{vector of input latitudes} + +\item{lonp}{vector of output longitudes} + +\item{latp}{vector of output latitudes} + +\item{method}{string indicating the interpolation method +("nearest" or "bilinear" (default))} +} +\value{ +The interpolated field. +} +\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. +} +\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") +} +\author{ +Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/man/dot-smooth.Rd b/man/dot-smooth.Rd new file mode 100644 index 00000000..a992c541 --- /dev/null +++ b/man/dot-smooth.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_RFTemp.R +\name{.smooth} +\alias{.smooth} +\title{Smoothening using convolution with a circular kernel} +\usage{ +.smooth(z, sdim) +} +\arguments{ +\item{z}{matrix with the input field to smoothen, with dimensions `c(ns, ns)`} + +\item{sdim}{the smoothing kernel radius in pixel} +} +\value{ +The smoothened field. +} +\description{ +The input field is convolved with a circular kernel with equal +weights. Takes into account missing values. +} +\examples{ +z <- rnorm(64 * 64) +dim(z) <- c(64, 64) +zs <- smooth(z, 8) +sd(zs) +# [1] 0.1334648 +} +\author{ +Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it} +} -- GitLab From 903e69a3d18f168d41335fd2ba5394f59be51892 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 18:17:10 +0200 Subject: [PATCH 33/41] noRd for internal functions --- R/CST_RFTemp.R | 32 +++++++++--------- man/dot-downscalet.Rd | 79 ------------------------------------------- man/dot-interp2d.Rd | 44 ------------------------ man/dot-smooth.Rd | 30 ---------------- 4 files changed, 16 insertions(+), 169 deletions(-) delete mode 100644 man/dot-downscalet.Rd delete mode 100644 man/dot-interp2d.Rd delete mode 100644 man/dot-smooth.Rd diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index a974b65f..a8e43c74 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -238,14 +238,14 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' @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)) - +#' 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, @@ -417,13 +417,13 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' ("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") - +#' 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) @@ -512,7 +512,7 @@ RFTemp <- function(data, lon, lat, oro, lonoro, latoro, #' zs <- smooth(z, 8) #' sd(zs) #' # [1] 0.1334648 - +#' @noRd .smooth <- function(z, sdim) { nsx <- dim(z)[1] nsy <- dim(z)[2] diff --git a/man/dot-downscalet.Rd b/man/dot-downscalet.Rd deleted file mode 100644 index 9b28e05e..00000000 --- a/man/dot-downscalet.Rd +++ /dev/null @@ -1,79 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/CST_RFTemp.R -\name{.downscalet} -\alias{.downscalet} -\title{Lapse-rate temperature correction downscaling} -\usage{ -.downscalet( - 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 -) -} -\arguments{ -\item{t}{matrix of input temperature data} - -\item{lon}{vector of input longitudes} - -\item{lat}{vector of input latitudes} - -\item{oro}{matrix of topographical elevations (in meters) -The destination downscaling area must be contained in the orography field.} - -\item{lono}{vector of orography longitudes} - -\item{lato}{vector of orography latitudes} - -\item{xlim}{vector of longitude bounds; the full input field is downscaled if `xlim` and `ylim` are not specified.} - -\item{ylim}{vector of latitude bounds} - -\item{radius}{smoothing radius expressed in longitude units -(default is half a large-scale pixel)} - -\item{lapse}{environmental lapse rate (in K/Km)} - -\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{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} -} -\value{ -A downscaled temperature matrix -} -\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. -} -\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)) -} -\author{ -Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it} -} diff --git a/man/dot-interp2d.Rd b/man/dot-interp2d.Rd deleted file mode 100644 index b03c6eee..00000000 --- a/man/dot-interp2d.Rd +++ /dev/null @@ -1,44 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/CST_RFTemp.R -\name{.interp2d} -\alias{.interp2d} -\title{Nearest neighbour interpolation} -\usage{ -.interp2d(z, lon, lat, lonp, latp, method = "bilinear") -} -\arguments{ -\item{z}{matrix with the input field to interpolate (assumed to -include also a third time dimension)} - -\item{lon}{vector of input longitudes} - -\item{lat}{vector of input latitudes} - -\item{lonp}{vector of output longitudes} - -\item{latp}{vector of output latitudes} - -\item{method}{string indicating the interpolation method -("nearest" or "bilinear" (default))} -} -\value{ -The interpolated field. -} -\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. -} -\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") -} -\author{ -Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it} -} diff --git a/man/dot-smooth.Rd b/man/dot-smooth.Rd deleted file mode 100644 index a992c541..00000000 --- a/man/dot-smooth.Rd +++ /dev/null @@ -1,30 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/CST_RFTemp.R -\name{.smooth} -\alias{.smooth} -\title{Smoothening using convolution with a circular kernel} -\usage{ -.smooth(z, sdim) -} -\arguments{ -\item{z}{matrix with the input field to smoothen, with dimensions `c(ns, ns)`} - -\item{sdim}{the smoothing kernel radius in pixel} -} -\value{ -The smoothened field. -} -\description{ -The input field is convolved with a circular kernel with equal -weights. Takes into account missing values. -} -\examples{ -z <- rnorm(64 * 64) -dim(z) <- c(64, 64) -zs <- smooth(z, 8) -sd(zs) -# [1] 0.1334648 -} -\author{ -Jost von Hardenberg, \email{j.vonhardenberg@isac.cnr.it} -} -- GitLab From 00c52890db03e81e265feb9afed403e4af0b3b12 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 18:20:51 +0200 Subject: [PATCH 34/41] time_dim parameter defined in RFTemp --- R/CST_RFTemp.R | 6 ++++-- man/CST_RFTemp.Rd | 2 ++ man/RFTemp.Rd | 2 ++ 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index a8e43c74..e710a4d9 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -29,6 +29,7 @@ #' @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 @@ -64,9 +65,9 @@ #' 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, +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, + nolapse = FALSE, verbose = FALSE, compute_delta = FALSE, method = "bilinear", delta = NULL) { if (!inherits(data, "s2dv_cube")) { @@ -127,6 +128,7 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse=6.5, #' @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 diff --git a/man/CST_RFTemp.Rd b/man/CST_RFTemp.Rd index a7beced6..d3b0c6a2 100644 --- a/man/CST_RFTemp.Rd +++ b/man/CST_RFTemp.Rd @@ -44,6 +44,8 @@ the full input field is downscaled if `xlim` and `ylim` are not specified.} \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} diff --git a/man/RFTemp.Rd b/man/RFTemp.Rd index 3945f571..d09d351a 100644 --- a/man/RFTemp.Rd +++ b/man/RFTemp.Rd @@ -50,6 +50,8 @@ the full input field is downscaled if `xlim` and `ylim` are not specified.} \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} -- GitLab From 194e6febbcb8653edfbc84c4c4c0fafd10fa0ffa Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 18:25:54 +0200 Subject: [PATCH 35/41] remove link tag for hyperlinks --- R/CST_RFTemp.R | 8 ++++---- man/CST_RFTemp.Rd | 4 ++-- man/RFTemp.Rd | 4 ++-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index e710a4d9..14c41642 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -10,10 +10,10 @@ #' 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 -#' \link[https://www.medscope-project.eu/the-project/deliverables-reports/] +#' [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 -#' \link[https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf] +#' [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} @@ -109,10 +109,10 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse = 6.5, #' 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 -#' \link[https://www.medscope-project.eu/the-project/deliverables-reports/] +#' [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 -#' \link[https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf] +#' [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 diff --git a/man/CST_RFTemp.Rd b/man/CST_RFTemp.Rd index d3b0c6a2..8ab5b6f3 100644 --- a/man/CST_RFTemp.Rd +++ b/man/CST_RFTemp.Rd @@ -97,10 +97,10 @@ 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 -\link[https://www.medscope-project.eu/the-project/deliverables-reports/] +[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 -\link[https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf] +[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 index d09d351a..628faf2d 100644 --- a/man/RFTemp.Rd +++ b/man/RFTemp.Rd @@ -100,10 +100,10 @@ res <- RFTemp(t, lon, lat, o, lono, lato, xlim=c(4, 8), ylim=c(43, 46), \references{ Method described in ERA4CS MEDSCOPE milestone M3.2: High-quality climate prediction data available to WP4 -\link[https://www.medscope-project.eu/the-project/deliverables-reports/] +[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 -\link[https://www.ecopotential-project.eu/images/ecopotential/documents/D8.1.pdf] +[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} -- GitLab From 1373e0bdc229140637e32a4f7baeaec6c767d88d Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 18:39:42 +0200 Subject: [PATCH 36/41] definition of parameter latoro and lonoro --- R/CST_RFTemp.R | 12 +++++++----- man/RFTemp.Rd | 12 ++++++------ 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 14c41642..9d8fd10b 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -120,6 +120,8 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse = 6.5, #' \code{lat_dim} parameters) #' @param lon Vector or array of longitudes. #' @param lat Vector or array of latitudes. +#' @param lon Vector or array of longitudes corresponding to the fine orography. +#' @param lat 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; @@ -150,21 +152,21 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse = 6.5, #' @import multiApply #' @examples #' # Generate simple synthetic data and downscale by factor 4 -#' t <- rnorm(7 * 6 * 4 * 3)*10 + 273.15 + 10 +#' 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 +#' 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) +#' 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, + nolapse = FALSE, verbose = FALSE, compute_delta = FALSE, method = "bilinear", delta = NULL) { # Check/detect time_dim diff --git a/man/RFTemp.Rd b/man/RFTemp.Rd index 628faf2d..f9fd5270 100644 --- a/man/RFTemp.Rd +++ b/man/RFTemp.Rd @@ -32,9 +32,9 @@ The input array is expected to have at least two dimensions named (these default names can be changed with the \code{lon_dim} and \code{lat_dim} parameters)} -\item{lon}{Vector or array of longitudes.} +\item{lon}{Vector or array of longitudes corresponding to the fine orography.} -\item{lat}{Vector or array of latitudes.} +\item{lat}{Vector or array of latitudes corresponding to the fine orography.} \item{oro}{Array containing fine-scale orography (in m) The destination downscaling area must be contained in the orography field.} @@ -86,16 +86,16 @@ 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 +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 +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) +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: -- GitLab From db54dfc54a3bc77e59440621bbbe885da7adef86 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 18:52:24 +0200 Subject: [PATCH 37/41] fixing latoro lonoro definition --- R/CST_RFTemp.R | 4 ++-- man/RFTemp.Rd | 8 ++++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/R/CST_RFTemp.R b/R/CST_RFTemp.R index 9d8fd10b..cebac85a 100644 --- a/R/CST_RFTemp.R +++ b/R/CST_RFTemp.R @@ -120,8 +120,8 @@ CST_RFTemp <- function(data, oro, xlim = NULL, ylim = NULL, lapse = 6.5, #' \code{lat_dim} parameters) #' @param lon Vector or array of longitudes. #' @param lat Vector or array of latitudes. -#' @param lon Vector or array of longitudes corresponding to the fine orography. -#' @param lat Vector or array of latitudes corresponding to the fine orography. +#' @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; diff --git a/man/RFTemp.Rd b/man/RFTemp.Rd index f9fd5270..106ae6e2 100644 --- a/man/RFTemp.Rd +++ b/man/RFTemp.Rd @@ -32,13 +32,17 @@ The input array is expected to have at least two dimensions named (these default names can be changed with the \code{lon_dim} and \code{lat_dim} parameters)} -\item{lon}{Vector or array of longitudes corresponding to the fine orography.} +\item{lon}{Vector or array of longitudes.} -\item{lat}{Vector or array of latitudes corresponding to the fine orography.} +\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.} -- GitLab From d21512ca91d5e97ba27ac5d9577e3813b49f338c Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 18:54:08 +0200 Subject: [PATCH 38/41] RF_Temp listed in NEWS --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 208fec39..19caa068 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 -- GitLab From 25aad46ed3a94e2f8d8ad8a2f2058bca31f52c9f Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 18:55:33 +0200 Subject: [PATCH 39/41] fix the link to the package CRAN site --- vignettes/ENSclustering_vignette.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/ENSclustering_vignette.Rmd b/vignettes/ENSclustering_vignette.Rmd index 52707d73..af8d5c6f 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 -- GitLab From 0175425d907ed493f99f72fb6d7be851c8ec185b Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 19:24:44 +0200 Subject: [PATCH 40/41] tests fixed in wrong branch --- tests/testthat/test-CST_WeatherRegimes.R | 42 ++++++++++++------------ 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/tests/testthat/test-CST_WeatherRegimes.R b/tests/testthat/test-CST_WeatherRegimes.R index 33e75e24..3f0c1dd3 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')) }) -- GitLab From 2b1ef57a6d505d407b26027d9d053755cbd363f5 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 29 Jun 2020 19:40:45 +0200 Subject: [PATCH 41/41] fixed tests --- tests/testthat/test-CST_WeatherRegimes.R | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/testthat/test-CST_WeatherRegimes.R b/tests/testthat/test-CST_WeatherRegimes.R index 3f0c1dd3..e2bc6feb 100644 --- a/tests/testthat/test-CST_WeatherRegimes.R +++ b/tests/testthat/test-CST_WeatherRegimes.R @@ -41,7 +41,7 @@ test_that("Sanity checks", { 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 @@ -53,9 +53,9 @@ test_that("Sanity checks", { 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 @@ -81,7 +81,7 @@ test_that("Sanity checks", { 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 @@ -90,10 +90,10 @@ test_that("Sanity checks", { 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')) }) -- GitLab