From c4d28f6782a14fba2977683b35c90168640ba4cc Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Tue, 7 May 2019 17:22:07 +0200 Subject: [PATCH 01/24] First Adamont functions adapted to CSTools --- R/Adamont_NearestNeighbors.R | 88 ++++++++++++++++++++++++++++++++++++ R/Adamont_Quantiles_S2D.R | 87 +++++++++++++++++++++++++++++++++++ 2 files changed, 175 insertions(+) create mode 100644 R/Adamont_NearestNeighbors.R create mode 100644 R/Adamont_Quantiles_S2D.R diff --git a/R/Adamont_NearestNeighbors.R b/R/Adamont_NearestNeighbors.R new file mode 100644 index 00000000..8fc3f49f --- /dev/null +++ b/R/Adamont_NearestNeighbors.R @@ -0,0 +1,88 @@ +#' ADAMONT Nearest Neighbors computes the distance between reference data grid centroid and SF data grid +#' +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#'@description This function computes the nearest neighbor for each reference data (lon, lat) point in the experiment dataset by computing the distance between the reference dataset grid and the experiment data. This is the first step in the ADAMONT method adapted from Verfaillie et al. (2018). +#' +#'@param method a string among three options ('ADA': standard ADAMONT distance, 'simple': lon/lat straight euclidian distance, 'radius': distance on the sphere) +#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment longitudes in \code{$lon} and latitudes in \code{$lat} +#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the reference data on a different grid, with longitudes in \code{$lon} and latitudes in \code{$lat}. +#' +#'@return NN +#' +#'@import s2dverification +#'@import ncdf4 + +NearestNeighbors <- function (method='ADA', exp, obs) { + + exp_lon <- exp$lon + exp_lat <- exp$lat + dim_exp_lon <- dim(exp_lon) + dim_exp_lat <- dim(exp_lat) + obs_lon <- obs$lon + obs_lat <- obs$lat + dim_obs_lon <- dim(obs_lon) + dim_obs_lat <- dim(obs_lat) + + # Flatten longitudes and latitudes in case of 2-D longitudes and latitudes (Lambert grids, etc.) + if ((length(dim_exp_lon)==2) & (length(dim_exp_lat)==2)){ + dim(exp_lon) <- c(dim_exp_lon[1]*dim_exp_lon[2]) + dim(exp_lat) <- c(dim_exp_lat[1]*dim_exp_lat[2]) + } + if ((length(dim_obs_lon)==2) & (length(dim_obs_lat)==2)){ + dim(obs_lon) <- c(dim_obs_lon[1]*dim_obs_lon[2]) + dim(obs_lat) <- c(dim_obs_lat[1]*dim_obs_lat[2]) + } + + # Now lat and lon arrays have 1 dimension, length npt (= nlat*nlon) + + OBS_grid <- cbind(obs_lon,obs_lat) + EXP_grid <- cbind(exp_lon,exp_lat) + + dist_min<-min_lon<-min_lat<-imin_lon<-imin_lat<-array(dim=nrow(OBS_grid)) + if (method == 'ADA'){ + C<-cos(OBS_grid[,2]*pi/180)^2 + for (i in 1:nrow(OBS_grid)){ + dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+C[i]*(OBS_grid[i,1]-EXP_grid[,1])^2 + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else if (method == 'simple'){ + for (i in 1:nrow(OBS_grid)){ + dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+(OBS_grid[i,1]-EXP_grid[,1])^2 + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else if (method == 'radius'){ + R <- 6371e3 # metres, Earth radius + EXP_gridr<-EXP_grid*pi/180 + OBS_gridr<-OBS_grid*pi/180 + for (i in 1:nrow(OBS_grid)){ + a<-sin((OBS_gridr[i,2]-EXP_gridr[,2])/2)^2 + cos(OBS_gridr[i,2])*cos(EXP_gridr[,2])*sin((OBS_gridr[i,1]-EXP_gridr[,1])/2)^2 + c<-2*atan2(sqrt(a),sqrt(1-a)) + dist<-R*c + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else { + stop("Adamont_NearestNeighbors supports method = 'ADA', 'simple' or 'radius' only.") + } + + # Reshape outputs to original grid?? Not really necessary? + #dim(min_lon)=dim_obs_lon + #dim(min_lat)=dim_obs_lat + + NN=list(dist_min=dist_min, min_lon=min_lon, min_lat=min_lat, imin_lon=imin_lon, imin_lat=imin_lat) + + return(NN) +} + diff --git a/R/Adamont_Quantiles_S2D.R b/R/Adamont_Quantiles_S2D.R new file mode 100644 index 00000000..a0111782 --- /dev/null +++ b/R/Adamont_Quantiles_S2D.R @@ -0,0 +1,87 @@ +#' ADAMONT Quantiles S2D computes quantiles for seasonal or decadal forecast data +#' +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#'@description This function computes quantiles on either reference or seasonal-to-decadal ensemble DAILY forecast data following the ADAMONT method adapted from Verfaillie et al. (2018). A prerequisite is to attribute each day and member of the S2D re-forecast to a given weather type / regime. +#' +#'@param exp an object of class \code{s2dv_cube} of lonlat data as returned by \code{CST_Load} function, containing the reference or seasonal forecast experiment data in \code{$data}, dimensions c(dataset,member,sdate,ftime,lat,lon) +#'@param var variable name (string) +#'@param wt an object of same member, sdate and ftime dimensions as exp, with weather type attribution (values 1 to Nwt where Nwt is the number of weather types) +#'@param tslice list of c(tmin,tmax) is the range of forecast times for which quantiles should be computed (default: NULL, values computed over all forecast times) +#'@param method a string among two options for extrapolation of extrema (default 'ADA': ADAMONT quantile computation, 'interp': interpolated quantile values) +#'@param seqq a list of quantiles (default to c(0.005, seq(from=0.01, to=0.00, by=0.01), 0.995)) +#'@param wetonly a logical (default to FALSE); which days quantiles of precipitation are computed (all days or wet) +#' +#'@return qdat an array with dimensions c(Nwt, Nq+2, lat, lon) where the Nq+2 values are [min value; Nq quantiles; max value]. +#' +#'@import s2dverification +#'@import ncdf4 +#'@import abind +#' +#' EXAMPLE TO BE INSERTED HERE... +#' +#' +#'exp <- runif(n=1*10*20*60*6*7,min=0.,max=5.) +#'dim(exp) <- c(dataset=1,member=10,sdate=20,ftime=60,lat=6,lon=7) +#'wt <- sample(1:4,1*10*20*60,replace=T) # Random sample for weather types numbered from 1 to 4 +#'dim(wt) <- c(dataset=1,member=10,sdate=20,ftime=60) + + +Quantiles_S2D <- function(exp, var, wt, tslice=NULL, method='ADA', seqq=c(0.005,seq(from=0.01, to=0.99, by=0.01),0.995), wetonly=FALSE) { + +# +# +# CHECKS NEEDED: +# WT VALUES SHOULD BE 1:NWT +# DIMENSIONS OF EXP AND WT DATA SHOULD BE CONSISTENT +# + +interp_quant<-function(data){ + out<-quantile(data,p=seqq,na.rm=TRUE) + return(c(min(data,na.rm=TRUE),out,max(data,na.rm=TRUE))) +} + +ADA_quant<-function(data){ + out<-quantile(data,p=seqq,type=1,na.rm=TRUE) + return(c(min(data,na.rm=TRUE),out,max(data,na.rm=TRUE))) +} + +# Which quantiles method to apply +quant <- switch(method,'ADA'=ADA_quant,'interp'=interp_quant) + +# Number of different weather types +wt_f <- wt +dim(wt_f) <- c(prod(dim(wt))) +Nwt <- length(unique(wt_f)) + +# Select ftimes for quantile computation +if (is.null(tslice)){ + warning('tslice=NULL: Computing quantiles over each forecast time') + dat <- exp$data +}else if (length(tslice) != 2) { + stop("Parameter 'tslice' must have length equal to 2 (c(tmin,tmax)) ") +}else { + dat <- exp$data[,,,tslice[1]:tslice[2],,,drop=FALSE] +} + +# Wet days only if variable is precipitation/snow and wetonly=TRUE +if ((var=='prlr') && (wetonly==TRUE)){ + norain=which(exp<1.) + dat[norain]=NA +} + +# Rearrange and flatten dimensions to use apply() function +dim_dat <- dim(dat) +dat_f <- aperm(dat,c(5,6,1,2,3,4)) +dim(dat_f) <- c(dim_dat[5],dim_dat[6],prod(dim_dat[1:4])) + +qexp <- array(dim=c(type=Nwt,quantiles=dim(array(seqq))+2,dim_dat[5],dim_dat[6])) + +for (reg in 1:Nwt){ + qexp[reg,,,] <- apply(dat_f[,,which(wt == reg)],c(1,2),quant) +} + +return(qexp) +} + + -- GitLab From 2603bf8e08baae5d61db8eecb74cc09e75ca4789 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Mon, 13 May 2019 11:41:02 +0200 Subject: [PATCH 02/24] Include Adamont_QQCorr.R -- not tested yetgit add Adamont_Quantiles_S2D.R --- R/Adamont_QQCorr.R | 52 +++++++++++++++++++++++++++++++++++++++ R/Adamont_Quantiles_S2D.R | 3 ++- 2 files changed, 54 insertions(+), 1 deletion(-) create mode 100644 R/Adamont_QQCorr.R diff --git a/R/Adamont_QQCorr.R b/R/Adamont_QQCorr.R new file mode 100644 index 00000000..3049613f --- /dev/null +++ b/R/Adamont_QQCorr.R @@ -0,0 +1,52 @@ +#' ADAMONT QQ Corr computes quantile-quantile correction of seasonal or decadal forecast data +#' +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#' +#'@param exp dimensions c(dataset,member,sdate,ftime,lat,lon) +#'@param wt dimensions c(dataset,member,sdate,ftime) +#'@param qexp dimensions c(Nwt, Nq+2, nlat, nlon) +#'@param qobs dimensions c(Nwt, Nq+2, nlat_o, nlon_o) +#'@param NN list (output from Adamont_NearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o) +#' +#'@import qmap + +QQ_Slope <- function(qexp, qobs) { + + Nwt <- dim(qexp)[1] + Nq <- dim(qexp)[2] - 2 + slope1 <- (qobs[,2]-qobs[,1])/(qexp[,2]-qexp[,1]) + slope2 <- (qobs[,Nq+2]-qobs[,Nq+1])/(qexp[,Nq+2]-qexp[,Nq+1]) + + return(list(s1=slope1,s2=slope2)) +} + +QQ_Corr <- function(exp, wt, qexp, qobs, NN) { + + # Data needs to be organized for use in qmap package + # Aim of correction is to have qexp data on the qobs grid, bias adjusted using QQ method + + dim_exp <- dim(exp) + dim_qobs <- dim(qobs) # CHECKS TO BE DONE: Nwt, Nq should be consistent between exp and obs + nlat_o <- dim_qobs[3] + nlon_o <- dim_qobs[4] + Nwt <- dim_qobs[1] + Nq <- dim_qobs[2]-2 + exp_corr <- array(dim=c(dim_exp[1:4],nlat_o,nlon_o)) + for (ilat in 1:nlat_o){ + for (ilon in 1:nlon_o){ + for (reg in 1:Nwt){ + modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq-2) + fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq-2) + qq_slope <- QQ_Slope(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[reg,,ilat,ilon]) + fobjj<-list(par=list(modq=modq, + fitq=as.matrix(qobs[reg,,ilat,ilon],nrow=Nq-2), + slope=as.matrix(c(qq_slope$s1,qq_slope$s2),nrow=2) ) + ) + corr_qmap <- doQmapRQUANT(exp[which(wt==reg),NN$imin_lat[ilat],NN$imin_lon[ilon]],fobjj,type="linear2") + exp_corr[which(wt==reg),ilat,ilon] <- corr_qmap + } + } + } +} + diff --git a/R/Adamont_Quantiles_S2D.R b/R/Adamont_Quantiles_S2D.R index a0111782..1d84c892 100644 --- a/R/Adamont_Quantiles_S2D.R +++ b/R/Adamont_Quantiles_S2D.R @@ -12,7 +12,7 @@ #'@param seqq a list of quantiles (default to c(0.005, seq(from=0.01, to=0.00, by=0.01), 0.995)) #'@param wetonly a logical (default to FALSE); which days quantiles of precipitation are computed (all days or wet) #' -#'@return qdat an array with dimensions c(Nwt, Nq+2, lat, lon) where the Nq+2 values are [min value; Nq quantiles; max value]. +#'@return qexp an array with dimensions c(Nwt, Nq+2, lat, lon) where the Nq+2 values are [min value; Nq quantiles; max value]. #' #'@import s2dverification #'@import ncdf4 @@ -26,6 +26,7 @@ #'wt <- sample(1:4,1*10*20*60,replace=T) # Random sample for weather types numbered from 1 to 4 #'dim(wt) <- c(dataset=1,member=10,sdate=20,ftime=60) +##### Input order should be changed... Quantiles_S2D <- function(exp, var, wt, tslice=NULL, method='ADA', seqq=c(0.005,seq(from=0.01, to=0.99, by=0.01),0.995), wetonly=FALSE) { -- GitLab From 30581e5ac6d29c3ced118601d5b0c4123f7c7840 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Mon, 13 May 2019 15:14:28 +0200 Subject: [PATCH 03/24] Output in Adamont_QQCorr.R --- R/Adamont_QQCorr.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/Adamont_QQCorr.R b/R/Adamont_QQCorr.R index 3049613f..1592be73 100644 --- a/R/Adamont_QQCorr.R +++ b/R/Adamont_QQCorr.R @@ -48,5 +48,7 @@ QQ_Corr <- function(exp, wt, qexp, qobs, NN) { } } } + + return(exp_corr) } -- GitLab From 8c66e0c86f211e01fbac452245542e5c3af1609d Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Tue, 14 May 2019 18:36:05 +0200 Subject: [PATCH 04/24] minor changes to Adamont_QQCorr.R --- R/Adamont_QQCorr.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/Adamont_QQCorr.R b/R/Adamont_QQCorr.R index 1592be73..99064633 100644 --- a/R/Adamont_QQCorr.R +++ b/R/Adamont_QQCorr.R @@ -16,7 +16,7 @@ QQ_Slope <- function(qexp, qobs) { Nwt <- dim(qexp)[1] Nq <- dim(qexp)[2] - 2 slope1 <- (qobs[,2]-qobs[,1])/(qexp[,2]-qexp[,1]) - slope2 <- (qobs[,Nq+2]-qobs[,Nq+1])/(qexp[,Nq+2]-qexp[,Nq+1]) + slope2 <- (qobs[,Nq+2]-qobs[,Nq+1])/(qexp[,Nq+2]-qexp[,Nq+1]) # CHANGE TO ONES[NWT, NQ+2] return(list(s1=slope1,s2=slope2)) } @@ -36,11 +36,11 @@ QQ_Corr <- function(exp, wt, qexp, qobs, NN) { for (ilat in 1:nlat_o){ for (ilon in 1:nlon_o){ for (reg in 1:Nwt){ - modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq-2) - fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq-2) + modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq+2) + fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq+2) qq_slope <- QQ_Slope(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[reg,,ilat,ilon]) fobjj<-list(par=list(modq=modq, - fitq=as.matrix(qobs[reg,,ilat,ilon],nrow=Nq-2), + fitq=fitq, slope=as.matrix(c(qq_slope$s1,qq_slope$s2),nrow=2) ) ) corr_qmap <- doQmapRQUANT(exp[which(wt==reg),NN$imin_lat[ilat],NN$imin_lon[ilon]],fobjj,type="linear2") -- GitLab From 4126b6415c011c415c6728260a855c9c2f7a4c23 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Thu, 16 May 2019 13:03:34 +0200 Subject: [PATCH 05/24] Very preliminary skeleton for Adamont_Analog.R --- R/Adamont_Analog.R | 83 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 R/Adamont_Analog.R diff --git a/R/Adamont_Analog.R b/R/Adamont_Analog.R new file mode 100644 index 00000000..79729327 --- /dev/null +++ b/R/Adamont_Analog.R @@ -0,0 +1,83 @@ +#' ADAMONT Analog finds analogous days in the reference dataset to experiment data based on +#' precipitation patterns +#' +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#' +#'@param exp_corr output from QQ_Corr, dimensions c(dataset=1,member,sdate,ftime,nlat_o,nlon_o) +#'@param wt corresponding weather types, dimensions c(dataset=1,member,sdate,ftime) +#'@param obs reference data, dimensions c(dataset=1,member,sdate,ftime,nlat_o,nlon_o) +#'@param wt_obs corresponding weather types, dimensions c(dataset=1,member,sdate,ftime) +#'@param var variable +#'@param method string indicating the method used for analog definition +#' This needs to be detailed more... Options rain1, rain01, patt_corr +#'@param exp_pr, needed if var is not precip +#'@param obs_pr, needed if var is not precip +#' +#'@import ncdf4, abind + + +Search_Analog_rain1 <- function() { + + return() +} + +Search_Analog_rain01 <- function() { + + return() +} + +Search_Analog_pattcorr <- function() { + + return() +} + + +Analog <- function(exp_corr, wt, obs, wt_obs, var, method='rain1', exp_pr=NULL, obs_pr=NULL) { + + # CHECKS NEEDED BEFORE RUNNING + # Do the dimensions correspond between obs_pr and obs? exp_pr and exp? + # INCLUDING CHECKS ON DIMENSIONS + # ONE OPTION IS TO GIVE THE NAME OF PRECIP VARIABLE IN THE WRAPPER FUNCTION + + + if (method == 'rain1'){ + + + } else if (method == 'rain01'){ + + } else if (method == 'pattcorr'){ + + } else { + stop("Adamont Analog function supports methods 'rain1', 'rain01', 'pattcorr' only") + } + + + + +# dim_exp <- dim(exp) +# dim_qobs <- dim(qobs) # CHECKS TO BE DONE: Nwt, Nq should be consistent between exp and obs +# nlat_o <- dim_qobs[3] +# nlon_o <- dim_qobs[4] +# Nwt <- dim_qobs[1] +# Nq <- dim_qobs[2]-2 +# exp_corr <- array(dim=c(dim_exp[1:4],nlat_o,nlon_o)) +# for (ilat in 1:nlat_o){ +# for (ilon in 1:nlon_o){ +# for (reg in 1:Nwt){ +# modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq+2) +# fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq+2) +# qq_slope <- QQ_Slope(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[reg,,ilat,ilon]) +# fobjj<-list(par=list(modq=modq, +# fitq=fitq, +# slope=as.matrix(c(qq_slope$s1,qq_slope$s2),nrow=2) ) +# ) +# corr_qmap <- doQmapRQUANT(exp[which(wt==reg),NN$imin_lat[ilat],NN$imin_lon[ilon]],fobjj,type="linear2") +# exp_corr[which(wt==reg),ilat,ilon] <- corr_qmap +# } +# } +# } + + return(analog) +} + -- GitLab From 7343e3a1fba86bc7a7c869ef7710dd9eec160b66 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Wed, 31 Jul 2019 14:47:57 +0200 Subject: [PATCH 06/24] First version of Adamont_Analog function --- R/Adamont_Analog.R | 140 +++++++++++++++++++++++++++++---------------- 1 file changed, 92 insertions(+), 48 deletions(-) diff --git a/R/Adamont_Analog.R b/R/Adamont_Analog.R index 79729327..8399409f 100644 --- a/R/Adamont_Analog.R +++ b/R/Adamont_Analog.R @@ -6,77 +6,121 @@ #' #'@param exp_corr output from QQ_Corr, dimensions c(dataset=1,member,sdate,ftime,nlat_o,nlon_o) #'@param wt corresponding weather types, dimensions c(dataset=1,member,sdate,ftime) -#'@param obs reference data, dimensions c(dataset=1,member,sdate,ftime,nlat_o,nlon_o) -#'@param wt_obs corresponding weather types, dimensions c(dataset=1,member,sdate,ftime) +#'@param obs reference data, dimensions c(dataset=1,member_obs=1,sdate,ftime,nlat_o,nlon_o) +#'@param wt_obs corresponding weather types, dimensions c(dataset=1,member_obs=1,sdate,ftime) #'@param var variable #'@param method string indicating the method used for analog definition #' This needs to be detailed more... Options rain1, rain01, patt_corr +#'@param thres real number indicating the precipitation threshold to define a rainy day #'@param exp_pr, needed if var is not precip #'@param obs_pr, needed if var is not precip #' #'@import ncdf4, abind -Search_Analog_rain1 <- function() { - - return() +Search_Analog_rain1 <- function(exp_day,obs,iwt,thres) { + # exp_day: exp_day[1,imb,sdate,ftime,,] + # obs: full obs array (all sdates and ftimes) + # iwt: indices corresponding to the weather type of exp_day + # thres: threshold for precipitation (rainy day or not) + analog5 <- array(dim=c(5)) + accuracy <- array(dim=length(iwt)) + for (i in iwt){ + accuracy[which(iwt==i)]<-sum((obs[i,,]>=thres)*(exp_day>=thres)) + } + best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] + vals_accuracy <- unique(best_accuracy) + nunique <- length(vals_accuracy) + count <- array(dim=nunique+1) + count[1] <- 0 + for (ss in 1:nunique){ + count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) + analog5[(count[ss]+1):count[ss+1]] <- sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) + } + return(analog5) } -Search_Analog_rain01 <- function() { - - return() +Search_Analog_rain01 <- function(exp_day,obs,iwt,thres) { + # exp_day: exp_day[1,imb,sdate,ftime,,] + # obs: full obs array (all sdates and ftimes) + # iwt: indices corresponding to the weather type of exp_day + # thres: threshold for precipitation (rainy day or not) + # Output: analog5 -> 5 indices corresponding to analog days of exp_day + analog5 <- array(dim=c(5)) + accuracy <- array(dim=length(iwt)) + for (i in iwt){ + accuracy[which(iwt==i)]<-sum(diag(table((obs[i,,]>=thres),(exp_day>=thres)))) + } + best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] + vals_accuracy <- unique(best_accuracy) + nunique <- length(vals_accuracy) + count <- array(dim=nunique+1) + count[1] <- 0 + for (ss in 1:nunique){ + count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) + analog5[(count[ss]+1):count[ss+1]] <- sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) + } + return(analog5) } -Search_Analog_pattcorr <- function() { - - return() +Search_Analog_pattcorr <- function(exp_day,obs,iwt,thres) { + # exp_day: exp_day[1,imb,sdate,ftime,,] + # obs: full obs array (all sdates and ftimes) + # iwt: indices corresponding to the weather type of exp_day + # thres: threshold for precipitation (rainy day or not) + # Output: analog5 -> 5 indices corresponding to analog days of exp_day + analog5 <- array(dim=c(5)) + accuracy <- array(dim=length(iwt)) + for (i in iwt){ + accuracy[which(iwt==i)]<-cor(as.vector(obs[i,,]),as.vector(exp_day)) + } + best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] + vals_accuracy <- unique(best_accuracy) + nunique <- length(vals_accuracy) + count <- array(dim=nunique+1) + count[1] <- 0 + for (ss in 1:nunique){ + count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) + analog5[(count[ss]+1):count[ss+1]] <- sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) + } + return(analog5) } -Analog <- function(exp_corr, wt, obs, wt_obs, var, method='rain1', exp_pr=NULL, obs_pr=NULL) { +Analog <- function(exp_corr, wt, obs, wt_obs, var, method='rain1', thres=1., exp_pr=NULL, obs_pr=NULL) { + + # WORK IN PROGRESS... + # CHECKS NEEDED BEFORE RUNNING # Do the dimensions correspond between obs_pr and obs? exp_pr and exp? # INCLUDING CHECKS ON DIMENSIONS + # CHECK ON VALUE OF PRECIPITATION # ONE OPTION IS TO GIVE THE NAME OF PRECIP VARIABLE IN THE WRAPPER FUNCTION - - if (method == 'rain1'){ - - - } else if (method == 'rain01'){ - - } else if (method == 'pattcorr'){ - - } else { - stop("Adamont Analog function supports methods 'rain1', 'rain01', 'pattcorr' only") + search_analog <- switch(method, 'rain1'=Search_Analog_rain1, 'rain01'=Search_Analog_rain01, 'pattcorr'=Search_Analog_pattcorr, stop("Adamont Analog function supports methods 'rain1', 'rain01', 'pattcorr' only")) + + analog <- array(dim=c(1,member=dim(exp_corr)[2],sdate=dim(exp_corr)[3],ftime=dim(exp_corr)[4],ndays=5)) + + for (imb in 1:dim(exp_corr)[2]){ + for (isdate in 1:dim(exp_corr)[3]){ + for (itime in 1:dim(exp_corr)[4]){ + iwt <- wt[1,imb,isdate,itime] + index <- which(wt_obs == iwt) + # No rain day case : pick five random days among corresponding weather type + if (sum(exp_pr[1,imb,isdate,itime,,]>=thres) == 0) { + analog[1,imb,isdate,itime,] <- obs[sample(index,5)] + } else { + # Find five closest analogs among corresponding weather type + if (length(index)!=0){ + obs_pr2 <- obs_pr + dim(obs_pr2) <- c(1*1*dim(obs_pr)[3]*dim(obs_pr[4]),dim(obs_pr[5]),dim(obs_pr[6])) + analog[1,imb,isdate,itime,] <- search_analog(exp_pr[1,imb,isdate,itime,,],obs_pr2,index,thres) + } + } + } + } } - - - - -# dim_exp <- dim(exp) -# dim_qobs <- dim(qobs) # CHECKS TO BE DONE: Nwt, Nq should be consistent between exp and obs -# nlat_o <- dim_qobs[3] -# nlon_o <- dim_qobs[4] -# Nwt <- dim_qobs[1] -# Nq <- dim_qobs[2]-2 -# exp_corr <- array(dim=c(dim_exp[1:4],nlat_o,nlon_o)) -# for (ilat in 1:nlat_o){ -# for (ilon in 1:nlon_o){ -# for (reg in 1:Nwt){ -# modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq+2) -# fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq+2) -# qq_slope <- QQ_Slope(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[reg,,ilat,ilon]) -# fobjj<-list(par=list(modq=modq, -# fitq=fitq, -# slope=as.matrix(c(qq_slope$s1,qq_slope$s2),nrow=2) ) -# ) -# corr_qmap <- doQmapRQUANT(exp[which(wt==reg),NN$imin_lat[ilat],NN$imin_lon[ilon]],fobjj,type="linear2") -# exp_corr[which(wt==reg),ilat,ilon] <- corr_qmap -# } -# } -# } return(analog) } -- GitLab From fdc083c86f2dc0ed24c46334406a98f78697efd9 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Wed, 2 Oct 2019 15:02:10 +0200 Subject: [PATCH 07/24] New update of Adamont functions for CSTools (Adamont_Analog.R still needs testing) --- R/Adamont_Analog.R | 17 +++++++++++++---- R/Adamont_QQCorr.R | 27 ++++++++++++++++++++------- R/Adamont_Quantiles_S2D.R | 11 ++++++----- 3 files changed, 39 insertions(+), 16 deletions(-) diff --git a/R/Adamont_Analog.R b/R/Adamont_Analog.R index 8399409f..b79cc32a 100644 --- a/R/Adamont_Analog.R +++ b/R/Adamont_Analog.R @@ -18,6 +18,15 @@ #'@import ncdf4, abind +My_Sample <- function(what,howmany){ + if (length(what) > 1){ + ff <- sample(what,howmany) + } else { + ff <- what + } + return(ff) +} + Search_Analog_rain1 <- function(exp_day,obs,iwt,thres) { # exp_day: exp_day[1,imb,sdate,ftime,,] # obs: full obs array (all sdates and ftimes) @@ -35,7 +44,7 @@ Search_Analog_rain1 <- function(exp_day,obs,iwt,thres) { count[1] <- 0 for (ss in 1:nunique){ count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) - analog5[(count[ss]+1):count[ss+1]] <- sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) + analog5[(count[ss]+1):count[ss+1]] <- My_Sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) } return(analog5) } @@ -58,7 +67,7 @@ Search_Analog_rain01 <- function(exp_day,obs,iwt,thres) { count[1] <- 0 for (ss in 1:nunique){ count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) - analog5[(count[ss]+1):count[ss+1]] <- sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) + analog5[(count[ss]+1):count[ss+1]] <- My_Sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) } return(analog5) } @@ -81,7 +90,7 @@ Search_Analog_pattcorr <- function(exp_day,obs,iwt,thres) { count[1] <- 0 for (ss in 1:nunique){ count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) - analog5[(count[ss]+1):count[ss+1]] <- sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) + analog5[(count[ss]+1):count[ss+1]] <- My_Sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) } return(analog5) } @@ -109,7 +118,7 @@ Analog <- function(exp_corr, wt, obs, wt_obs, var, method='rain1', thres=1., exp index <- which(wt_obs == iwt) # No rain day case : pick five random days among corresponding weather type if (sum(exp_pr[1,imb,isdate,itime,,]>=thres) == 0) { - analog[1,imb,isdate,itime,] <- obs[sample(index,5)] + analog[1,imb,isdate,itime,] <- obs[My_Sample(index,5)] } else { # Find five closest analogs among corresponding weather type if (length(index)!=0){ diff --git a/R/Adamont_QQCorr.R b/R/Adamont_QQCorr.R index 99064633..208326ce 100644 --- a/R/Adamont_QQCorr.R +++ b/R/Adamont_QQCorr.R @@ -16,12 +16,12 @@ QQ_Slope <- function(qexp, qobs) { Nwt <- dim(qexp)[1] Nq <- dim(qexp)[2] - 2 slope1 <- (qobs[,2]-qobs[,1])/(qexp[,2]-qexp[,1]) - slope2 <- (qobs[,Nq+2]-qobs[,Nq+1])/(qexp[,Nq+2]-qexp[,Nq+1]) # CHANGE TO ONES[NWT, NQ+2] + slope2 <- matrix(1,nrow=Nwt,ncol=Nq+2) return(list(s1=slope1,s2=slope2)) } -QQ_Corr <- function(exp, wt, qexp, qobs, NN) { +QQ_Corr <- function(exp, wt, qexp, qobs, NN, wetonly=FALSE, thres=1.) { # Data needs to be organized for use in qmap package # Aim of correction is to have qexp data on the qobs grid, bias adjusted using QQ method @@ -32,23 +32,36 @@ QQ_Corr <- function(exp, wt, qexp, qobs, NN) { nlon_o <- dim_qobs[4] Nwt <- dim_qobs[1] Nq <- dim_qobs[2]-2 - exp_corr <- array(dim=c(dim_exp[1:4],nlat_o,nlon_o)) + + exp_corr <- array(dim=c(nlat_o,nlon_o,dim_exp[1:4])) + exp <- aperm(exp,c(5,6,1,2,3,4)) + + # REORDER DATA FOR USE OF APPLY FUNCTION AND WHICH WITHOUT PROBLEMS... + for (ilat in 1:nlat_o){ + print(paste("ilat: ",ilat,sep="")) for (ilon in 1:nlon_o){ + # Test outside loop on reg + qq_slope <- QQ_Slope(qexp[,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[,,ilat,ilon]) for (reg in 1:Nwt){ modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq+2) fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq+2) - qq_slope <- QQ_Slope(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[reg,,ilat,ilon]) + #qq_slope <- QQ_Slope(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[reg,,ilat,ilon]) fobjj<-list(par=list(modq=modq, fitq=fitq, - slope=as.matrix(c(qq_slope$s1,qq_slope$s2),nrow=2) ) + slope=as.matrix(c(qq_slope$s1[reg],qq_slope$s2[reg]),nrow=2) ) ) - corr_qmap <- doQmapRQUANT(exp[which(wt==reg),NN$imin_lat[ilat],NN$imin_lon[ilon]],fobjj,type="linear2") - exp_corr[which(wt==reg),ilat,ilon] <- corr_qmap + attr(fobjj$par$slope,"dimnames")<-list(c("lower","upper")) + class(fobjj) <- "fitQmapRQUANT" + indices <- which(wt==reg,arr.ind=TRUE) + for (ind in 1:dim(indices)[1]){ + exp_corr[ilat,ilon,indices[ind,1],indices[ind,2],indices[ind,3],indices[ind,4]] <- doQmapRQUANT(exp[NN$imin_lat[ilat],NN$imin_lon[ilon],indices[ind,1],indices[ind,2],indices[ind,3],indices[ind,4]],fobjj,type="linear",wet.day=ifelse(wetonly==TRUE,thres,FALSE)) + } } } } + exp_corr <- aperm(exp_corr,c(3,4,5,6,1,2)) return(exp_corr) } diff --git a/R/Adamont_Quantiles_S2D.R b/R/Adamont_Quantiles_S2D.R index 1d84c892..7a072d5a 100644 --- a/R/Adamont_Quantiles_S2D.R +++ b/R/Adamont_Quantiles_S2D.R @@ -50,11 +50,6 @@ ADA_quant<-function(data){ # Which quantiles method to apply quant <- switch(method,'ADA'=ADA_quant,'interp'=interp_quant) -# Number of different weather types -wt_f <- wt -dim(wt_f) <- c(prod(dim(wt))) -Nwt <- length(unique(wt_f)) - # Select ftimes for quantile computation if (is.null(tslice)){ warning('tslice=NULL: Computing quantiles over each forecast time') @@ -63,8 +58,14 @@ if (is.null(tslice)){ stop("Parameter 'tslice' must have length equal to 2 (c(tmin,tmax)) ") }else { dat <- exp$data[,,,tslice[1]:tslice[2],,,drop=FALSE] + wt <- wt[,,,tslice[1]:tslice[2],drop=FALSE] } +# Number of different weather types +wt_f <- wt +dim(wt_f) <- c(prod(dim(wt))) +Nwt <- length(unique(wt_f)) + # Wet days only if variable is precipitation/snow and wetonly=TRUE if ((var=='prlr') && (wetonly==TRUE)){ norain=which(exp<1.) -- GitLab From d83fe47e3c3f1f309a08c7fc8fa85239483cbe8d Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Thu, 9 Jul 2020 11:52:53 +0200 Subject: [PATCH 08/24] Changes in Adamont_QQCorr.R function --- R/Adamont_QQCorr.R | 4 ---- 1 file changed, 4 deletions(-) diff --git a/R/Adamont_QQCorr.R b/R/Adamont_QQCorr.R index 208326ce..177f4d63 100644 --- a/R/Adamont_QQCorr.R +++ b/R/Adamont_QQCorr.R @@ -36,17 +36,13 @@ QQ_Corr <- function(exp, wt, qexp, qobs, NN, wetonly=FALSE, thres=1.) { exp_corr <- array(dim=c(nlat_o,nlon_o,dim_exp[1:4])) exp <- aperm(exp,c(5,6,1,2,3,4)) - # REORDER DATA FOR USE OF APPLY FUNCTION AND WHICH WITHOUT PROBLEMS... - for (ilat in 1:nlat_o){ - print(paste("ilat: ",ilat,sep="")) for (ilon in 1:nlon_o){ # Test outside loop on reg qq_slope <- QQ_Slope(qexp[,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[,,ilat,ilon]) for (reg in 1:Nwt){ modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq+2) fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq+2) - #qq_slope <- QQ_Slope(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[reg,,ilat,ilon]) fobjj<-list(par=list(modq=modq, fitq=fitq, slope=as.matrix(c(qq_slope$s1[reg],qq_slope$s2[reg]),nrow=2) ) -- GitLab From 90cfc0a3cdf8c5dd3f9bdd821ae7215a7993c5c2 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Tue, 1 Sep 2020 16:51:25 +0200 Subject: [PATCH 09/24] Adjustments to Adamont_Quantiles_S2D.R --- R/Adamont_Quantiles_S2D.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/Adamont_Quantiles_S2D.R b/R/Adamont_Quantiles_S2D.R index 7a072d5a..179ba34d 100644 --- a/R/Adamont_Quantiles_S2D.R +++ b/R/Adamont_Quantiles_S2D.R @@ -6,7 +6,7 @@ #' #'@param exp an object of class \code{s2dv_cube} of lonlat data as returned by \code{CST_Load} function, containing the reference or seasonal forecast experiment data in \code{$data}, dimensions c(dataset,member,sdate,ftime,lat,lon) #'@param var variable name (string) -#'@param wt an object of same member, sdate and ftime dimensions as exp, with weather type attribution (values 1 to Nwt where Nwt is the number of weather types) +#'@param wt an object of same member, sdate and ftime dimensions as exp, with weather type attribution (values 1 to Nwt where Nwt is the number of weather types); can be $cluster output of \code{CST_WeatherRegimes} #'@param tslice list of c(tmin,tmax) is the range of forecast times for which quantiles should be computed (default: NULL, values computed over all forecast times) #'@param method a string among two options for extrapolation of extrema (default 'ADA': ADAMONT quantile computation, 'interp': interpolated quantile values) #'@param seqq a list of quantiles (default to c(0.005, seq(from=0.01, to=0.00, by=0.01), 0.995)) -- GitLab From a54dc08eb321756edbd555b27bf0aed5cf4bbcad Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Thu, 10 Sep 2020 12:35:02 +0200 Subject: [PATCH 10/24] Included function Adamont_TimeDisaggregate.R + bugfixes on Adamont_Analog.R --- R/Adamont_Analog.R | 9 +- R/Adamont_TimeDisaggregate.R | 155 +++++++++++++++++++++++++++++++++++ 2 files changed, 162 insertions(+), 2 deletions(-) create mode 100644 R/Adamont_TimeDisaggregate.R diff --git a/R/Adamont_Analog.R b/R/Adamont_Analog.R index b79cc32a..6d998db4 100644 --- a/R/Adamont_Analog.R +++ b/R/Adamont_Analog.R @@ -100,7 +100,12 @@ Analog <- function(exp_corr, wt, obs, wt_obs, var, method='rain1', thres=1., exp # WORK IN PROGRESS... - + if ((var != "prec") & (is.null(exp_pr) | is.null(obs_pr))) { + stop("Analog is selected according to precipitation pattern: exp_pr and obs_pr must be provided") + } else { + exp_pr <- exp_corr + obs_pr <- obs + } # CHECKS NEEDED BEFORE RUNNING # Do the dimensions correspond between obs_pr and obs? exp_pr and exp? # INCLUDING CHECKS ON DIMENSIONS @@ -123,7 +128,7 @@ Analog <- function(exp_corr, wt, obs, wt_obs, var, method='rain1', thres=1., exp # Find five closest analogs among corresponding weather type if (length(index)!=0){ obs_pr2 <- obs_pr - dim(obs_pr2) <- c(1*1*dim(obs_pr)[3]*dim(obs_pr[4]),dim(obs_pr[5]),dim(obs_pr[6])) + dim(obs_pr2) <- c(1*1*dim(obs_pr)[3]*dim(obs_pr)[4],dim(obs_pr)[5],dim(obs_pr)[6]) analog[1,imb,isdate,itime,] <- search_analog(exp_pr[1,imb,isdate,itime,,],obs_pr2,index,thres) } } diff --git a/R/Adamont_TimeDisaggregate.R b/R/Adamont_TimeDisaggregate.R new file mode 100644 index 00000000..a2ab7c22 --- /dev/null +++ b/R/Adamont_TimeDisaggregate.R @@ -0,0 +1,155 @@ +#' ADAMONT TimeDisaggregate uses analogous days from ADAMONT Analog and higher frequency time series in the reference dataset to experiment data based on +#' precipitation patterns +#' +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#' +#'@param exp_corr output from QQ_Corr, dimensions c(dataset=1,member,sdate,ftime,nlat_o,nlon_o) +#'@param analog output from Analog, dimensions c(dataset=1,member,sdate,ftime,5) (5 analogs) +#'@param obs reference subdaily data, dimensions c(dataset=1,member_obs=1,sdate,ftime,subd,nlat_o,nlon_o) +#' subd is number of subdaily time steps in obs data (24 for hourly, 4 for 6-hourly, etc.) +#'@param obs_daily reference daily data, dimensions c(dataset=1,member_obs=1,sdate,ftime,nlat_o,nlon_o) +#'@param var variable name +#'@param indx index of analog day used for precip correction, dimensions c(dataset=1,member,sdate,ftime) +#' (if var=prec, NULL required, if not, should be provided) +#'@param thres precipitation threshold for definition of rainy day +#' +#'@import ncdf4, abind + + +TimeDisaggregate <- function(exp_corr, analog, obs, obs_daily, var, indx, thres) { + + Disagg <- function(mod, ref, var){ + tres <- length(ref) # Length of subdaily data (24: hourly, 4: 6-hourly, etc.) + D_mod <- array(dim=c(tres)) + if (mean(ref) > 10e-10){ + correction1 = c("Qair", "SWdown", "prec", "LWdown") + correction2 = c("uas", "vas") + # CHANGE TO STANDARD VARIABLE NAMES ABOVE... + if (var %in% correction1){ + a <- mod/mean(ref) + b <- 0 + } else if (var %in% correction2) { + a <- mod/mean(ref) + if (a <= 1) { + b = 0 + } else { + a = 1 + b = mod - mean(ref) + } + } else { # Standard case such as temperature + a <- 1 + b <- mod-mean(ref) + } + } else { # mean(ref) <= 10e-10 + correction1 = c("LWdown") + correction2 = c("uas", "vas") + if (var %in% correction2){ + a = 0 + b = rep(mod, length(ref)) + } else if (var %in% correction1) { + a <- mod/mean(ref) + b <- 0 + } else { + a <- 0 + b <- 0 + } + } + D_mod <- a*ref + b + return(D_mod) + } + + SelPrecAnalog <- function(mod, ref, analog5, thres) { + # Function to select which of the 5 possible analogs is chosen + # mod: array with corrected daily precip data at grid point + # ref: corresponding obs array with subdaily precip data at grid point for five analog days + # analog5: 5 analog days + # thres: threshold for precipitation + day_indx <- 0 + repeat { + day_indx <- day_indx+1 + ind <- analog5[day_indx] + if (day_indx == 5 | (sum(is.na(ref[day_indx,]))==0 & ((mod*86400 < thres ) | (mod*86400 >= thres & mean(ref[day_indx,])*86400 >= thres)))) { + break + } + } + #return(day_indx) + return(ind) + } + + dim_exp <- dim(exp_corr) + exp_dt <- array(dim=c(subd=dim(obs)[5],dim_exp[1:6])) + # obs dimensions: c(dataset=1,member_obs=1,sdate,ftime,subd,nlat_o,nlon_o) + obs2 <- aperm(obs,c(5,1,2,3,4,6,7)) # Start with subdaily dimension + # Flatten obs2 array + dim(obs2) <- c(dim(obs2)[1],prod(dim(obs2)[2:5]),dim(obs2)[6],dim(obs2)[7]) + + if (var != "prec"){ + if (is.null(indx)){ + stop("TimeDisaggregate MUST first be run on precipitation to select appropriate analog index indx") + } else { + # RUN Disagg function on data using analog day indx + for (imb in 1:dim(exp_corr)[2]){ + for (sd in 1:dim(exp_corr)[3]){ + for (ltime in 1:dim(exp_corr)[4]){ + for (ilat in 1:dim(exp_corr)[5]){ + for (ilon in 1:dim(exp_corr)[6]){ + exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) + } + } + } + } + } + } + } else { + if (is.null(indx)){ + # Select appropriate analog among 5 possible -> indx array + indx <- array(dim=dim_exp) +# for (imb in 1:dim(exp_corr)[2]){ +# indx[1,imb,,,,] <- mapply(SelPrecAnalog,exp_corr[1,imb,,,,],obs2,analog[1,imb,,,],thres) +# } + for (imb in 1:dim(exp_corr)[2]){ + for (sd in 1:dim(exp_corr)[3]){ + for (ltime in 1:dim(exp_corr)[4]){ + for (ilat in 1:dim(exp_corr)[5]){ + for (ilon in 1:dim(exp_corr)[6]){ + indx[1,imb,sd,ltime,ilat,ilon] <- SelPrecAnalog(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,,ilat,ilon],analog[1,imb,sd,ltime,],thres) + } + } + } + } + } + # RUN Disagg function on data + for (imb in 1:dim(exp_corr)[2]){ + for (sd in 1:dim(exp_corr)[3]){ + for (ltime in 1:dim(exp_corr)[4]){ + for (ilat in 1:dim(exp_corr)[5]){ + for (ilon in 1:dim(exp_corr)[6]){ + exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) + } + } + } + } + } + # HERE should be same logic as above: apply Disagg function on all dimensions except sub-daily. + + } else { + warning("Using previously defined indx array for TimeDisaggregate: precipitation disaggregation is run again") + for (imb in 1:dim(exp_corr)[2]){ + for (sd in 1:dim(exp_corr)[3]){ + for (ltime in 1:dim(exp_corr)[4]){ + for (ilat in 1:dim(exp_corr)[5]){ + for (ilon in 1:dim(exp_corr)[6]){ + exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) + } + } + } + } + } + } + + } + + return(list(disagg=exp_dt, indx=indx)) +} + -- GitLab From e0794577beed4d9aff84a80cc2723cfc43e90214 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Fri, 11 Sep 2020 12:23:07 +0200 Subject: [PATCH 11/24] Edits to conform to R code style --- R/Adamont_Analog.R | 207 ++++++++++++++-------------- R/Adamont_QQCorr.R | 86 ++++++------ R/Adamont_Quantiles_S2D.R | 2 - R/Adamont_TimeDisaggregate.R | 259 +++++++++++++++++------------------ 4 files changed, 266 insertions(+), 288 deletions(-) diff --git a/R/Adamont_Analog.R b/R/Adamont_Analog.R index 6d998db4..b3aa8f25 100644 --- a/R/Adamont_Analog.R +++ b/R/Adamont_Analog.R @@ -18,124 +18,117 @@ #'@import ncdf4, abind -My_Sample <- function(what,howmany){ - if (length(what) > 1){ - ff <- sample(what,howmany) - } else { - ff <- what - } - return(ff) +.MySample <- function(what,howmany){ + if (length(what) > 1){ + ff <- sample(what,howmany) + } else { + ff <- what + } + return(ff) } -Search_Analog_rain1 <- function(exp_day,obs,iwt,thres) { - # exp_day: exp_day[1,imb,sdate,ftime,,] - # obs: full obs array (all sdates and ftimes) - # iwt: indices corresponding to the weather type of exp_day - # thres: threshold for precipitation (rainy day or not) - analog5 <- array(dim=c(5)) - accuracy <- array(dim=length(iwt)) - for (i in iwt){ - accuracy[which(iwt==i)]<-sum((obs[i,,]>=thres)*(exp_day>=thres)) - } - best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] - vals_accuracy <- unique(best_accuracy) - nunique <- length(vals_accuracy) - count <- array(dim=nunique+1) - count[1] <- 0 - for (ss in 1:nunique){ - count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) - analog5[(count[ss]+1):count[ss+1]] <- My_Sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) - } - return(analog5) +.SearchAnalogRain1 <- function(exp_day,obs,iwt,thres) { + ## exp_day: exp_day[1,imb,sdate,ftime,,] + ## obs: full obs array (all sdates and ftimes) + ## iwt: indices corresponding to the weather type of exp_day + ## thres: threshold for precipitation (rainy day or not) + analog5 <- array(dim=c(5)) + accuracy <- array(dim=length(iwt)) + for (i in iwt){ + accuracy[which(iwt==i)]<-sum((obs[i,,]>=thres)*(exp_day>=thres)) + } + best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] + vals_accuracy <- unique(best_accuracy) + nunique <- length(vals_accuracy) + count <- array(dim=nunique+1) + count[1] <- 0 + for (ss in 1:nunique){ + count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) + analog5[(count[ss]+1):count[ss+1]] <- MySample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) + } + return(analog5) } -Search_Analog_rain01 <- function(exp_day,obs,iwt,thres) { - # exp_day: exp_day[1,imb,sdate,ftime,,] - # obs: full obs array (all sdates and ftimes) - # iwt: indices corresponding to the weather type of exp_day - # thres: threshold for precipitation (rainy day or not) - # Output: analog5 -> 5 indices corresponding to analog days of exp_day - analog5 <- array(dim=c(5)) - accuracy <- array(dim=length(iwt)) - for (i in iwt){ - accuracy[which(iwt==i)]<-sum(diag(table((obs[i,,]>=thres),(exp_day>=thres)))) - } - best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] - vals_accuracy <- unique(best_accuracy) - nunique <- length(vals_accuracy) - count <- array(dim=nunique+1) - count[1] <- 0 - for (ss in 1:nunique){ - count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) - analog5[(count[ss]+1):count[ss+1]] <- My_Sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) - } - return(analog5) +.SearchAnalogRain01 <- function(exp_day,obs,iwt,thres) { + ## exp_day: exp_day[1,imb,sdate,ftime,,] + ## obs: full obs array (all sdates and ftimes) + ## iwt: indices corresponding to the weather type of exp_day + ## thres: threshold for precipitation (rainy day or not) + ## Output: analog5 -> 5 indices corresponding to analog days of exp_day + analog5 <- array(dim=c(5)) + accuracy <- array(dim=length(iwt)) + for (i in iwt){ + accuracy[which(iwt==i)]<-sum(diag(table((obs[i,,]>=thres),(exp_day>=thres)))) + } + best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] + vals_accuracy <- unique(best_accuracy) + nunique <- length(vals_accuracy) + count <- array(dim=nunique+1) + count[1] <- 0 + for (ss in 1:nunique){ + count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) + analog5[(count[ss]+1):count[ss+1]] <- MySample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) + } + return(analog5) } -Search_Analog_pattcorr <- function(exp_day,obs,iwt,thres) { - # exp_day: exp_day[1,imb,sdate,ftime,,] - # obs: full obs array (all sdates and ftimes) - # iwt: indices corresponding to the weather type of exp_day - # thres: threshold for precipitation (rainy day or not) - # Output: analog5 -> 5 indices corresponding to analog days of exp_day - analog5 <- array(dim=c(5)) - accuracy <- array(dim=length(iwt)) - for (i in iwt){ - accuracy[which(iwt==i)]<-cor(as.vector(obs[i,,]),as.vector(exp_day)) - } - best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] - vals_accuracy <- unique(best_accuracy) - nunique <- length(vals_accuracy) - count <- array(dim=nunique+1) - count[1] <- 0 - for (ss in 1:nunique){ - count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) - analog5[(count[ss]+1):count[ss+1]] <- My_Sample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) - } - return(analog5) +.SearchAnalogPattcorr <- function(exp_day,obs,iwt,thres) { + ## exp_day: exp_day[1,imb,sdate,ftime,,] + ## obs: full obs array (all sdates and ftimes) + ## iwt: indices corresponding to the weather type of exp_day + ## thres: threshold for precipitation (rainy day or not) + ## Output: analog5 -> 5 indices corresponding to analog days of exp_day + analog5 <- array(dim=c(5)) + accuracy <- array(dim=length(iwt)) + for (i in iwt){ + accuracy[which(iwt==i)]<-cor(as.vector(obs[i,,]),as.vector(exp_day)) + } + best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] + vals_accuracy <- unique(best_accuracy) + nunique <- length(vals_accuracy) + count <- array(dim=nunique+1) + count[1] <- 0 + for (ss in 1:nunique){ + count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) + analog5[(count[ss]+1):count[ss+1]] <- MySample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) + } + return(analog5) } -Analog <- function(exp_corr, wt, obs, wt_obs, var, method='rain1', thres=1., exp_pr=NULL, obs_pr=NULL) { - - - # WORK IN PROGRESS... - if ((var != "prec") & (is.null(exp_pr) | is.null(obs_pr))) { - stop("Analog is selected according to precipitation pattern: exp_pr and obs_pr must be provided") - } else { - exp_pr <- exp_corr - obs_pr <- obs - } - # CHECKS NEEDED BEFORE RUNNING - # Do the dimensions correspond between obs_pr and obs? exp_pr and exp? - # INCLUDING CHECKS ON DIMENSIONS - # CHECK ON VALUE OF PRECIPITATION - # ONE OPTION IS TO GIVE THE NAME OF PRECIP VARIABLE IN THE WRAPPER FUNCTION - - search_analog <- switch(method, 'rain1'=Search_Analog_rain1, 'rain01'=Search_Analog_rain01, 'pattcorr'=Search_Analog_pattcorr, stop("Adamont Analog function supports methods 'rain1', 'rain01', 'pattcorr' only")) - - analog <- array(dim=c(1,member=dim(exp_corr)[2],sdate=dim(exp_corr)[3],ftime=dim(exp_corr)[4],ndays=5)) - - for (imb in 1:dim(exp_corr)[2]){ - for (isdate in 1:dim(exp_corr)[3]){ - for (itime in 1:dim(exp_corr)[4]){ - iwt <- wt[1,imb,isdate,itime] - index <- which(wt_obs == iwt) - # No rain day case : pick five random days among corresponding weather type - if (sum(exp_pr[1,imb,isdate,itime,,]>=thres) == 0) { - analog[1,imb,isdate,itime,] <- obs[My_Sample(index,5)] - } else { - # Find five closest analogs among corresponding weather type - if (length(index)!=0){ +AdamontAnalog <- function(exp_corr, wt, obs, wt_obs, var, method='rain1', thres=1., exp_pr=NULL, obs_pr=NULL) { + ## CHECKS TO BE ADDED + ## Do the dimensions correspond between obs_pr and obs? exp_pr and exp? + ## INCLUDING CHECKS ON DIMENSIONS + ## CHECK ON VALUE OF PRECIPITATION + ## ONE OPTION IS TO GIVE THE NAME OF PRECIP VARIABLE IN THE WRAPPER FUNCTION + if ((var != "prec") & (is.null(exp_pr) | is.null(obs_pr))) { + stop("Analog is selected according to precipitation pattern: exp_pr and obs_pr must be provided") + } else { + exp_pr <- exp_corr + obs_pr <- obs + } + search_analog <- switch(method, 'rain1'=SearchAnalogRain1, 'rain01'=SearchAnalogRain01, 'pattcorr'=SearchAnalogPattcorr, stop("Adamont Analog function supports methods 'rain1', 'rain01', 'pattcorr' only")) + analog <- array(dim=c(1,member=dim(exp_corr)[2],sdate=dim(exp_corr)[3],ftime=dim(exp_corr)[4],ndays=5)) + for (imb in 1:dim(exp_corr)[2]){ + for (isdate in 1:dim(exp_corr)[3]){ + for (itime in 1:dim(exp_corr)[4]){ + iwt <- wt[1,imb,isdate,itime] + index <- which(wt_obs == iwt) + ## No rain day case : pick five random days among corresponding weather type + if (sum(exp_pr[1,imb,isdate,itime,,]>=thres) == 0) { + analog[1,imb,isdate,itime,] <- obs[MySample(index,5)] + } else { + ## Find five closest analogs among corresponding weather type + if (length(index)!=0){ obs_pr2 <- obs_pr dim(obs_pr2) <- c(1*1*dim(obs_pr)[3]*dim(obs_pr)[4],dim(obs_pr)[5],dim(obs_pr)[6]) analog[1,imb,isdate,itime,] <- search_analog(exp_pr[1,imb,isdate,itime,,],obs_pr2,index,thres) - } - } - } - } - } - - return(analog) + } + } + } + } + } + return(analog) } diff --git a/R/Adamont_QQCorr.R b/R/Adamont_QQCorr.R index 177f4d63..c6e3eeac 100644 --- a/R/Adamont_QQCorr.R +++ b/R/Adamont_QQCorr.R @@ -1,4 +1,4 @@ -#' ADAMONT QQ Corr computes quantile-quantile correction of seasonal or decadal forecast data +#' ADAMONT QQCorr computes quantile-quantile correction of seasonal or decadal forecast data #' #'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version #'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation @@ -11,53 +11,45 @@ #' #'@import qmap -QQ_Slope <- function(qexp, qobs) { - - Nwt <- dim(qexp)[1] - Nq <- dim(qexp)[2] - 2 - slope1 <- (qobs[,2]-qobs[,1])/(qexp[,2]-qexp[,1]) - slope2 <- matrix(1,nrow=Nwt,ncol=Nq+2) - - return(list(s1=slope1,s2=slope2)) +.QQSlope <- function(qexp, qobs) { + Nwt <- dim(qexp)[1] + Nq <- dim(qexp)[2] - 2 + slope1 <- (qobs[,2]-qobs[,1])/(qexp[,2]-qexp[,1]) + slope2 <- matrix(1,nrow=Nwt,ncol=Nq+2) + return(list(s1=slope1,s2=slope2)) } -QQ_Corr <- function(exp, wt, qexp, qobs, NN, wetonly=FALSE, thres=1.) { - - # Data needs to be organized for use in qmap package - # Aim of correction is to have qexp data on the qobs grid, bias adjusted using QQ method - - dim_exp <- dim(exp) - dim_qobs <- dim(qobs) # CHECKS TO BE DONE: Nwt, Nq should be consistent between exp and obs - nlat_o <- dim_qobs[3] - nlon_o <- dim_qobs[4] - Nwt <- dim_qobs[1] - Nq <- dim_qobs[2]-2 - - exp_corr <- array(dim=c(nlat_o,nlon_o,dim_exp[1:4])) - exp <- aperm(exp,c(5,6,1,2,3,4)) - - for (ilat in 1:nlat_o){ - for (ilon in 1:nlon_o){ - # Test outside loop on reg - qq_slope <- QQ_Slope(qexp[,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[,,ilat,ilon]) - for (reg in 1:Nwt){ - modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq+2) - fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq+2) - fobjj<-list(par=list(modq=modq, - fitq=fitq, - slope=as.matrix(c(qq_slope$s1[reg],qq_slope$s2[reg]),nrow=2) ) - ) - attr(fobjj$par$slope,"dimnames")<-list(c("lower","upper")) - class(fobjj) <- "fitQmapRQUANT" - indices <- which(wt==reg,arr.ind=TRUE) - for (ind in 1:dim(indices)[1]){ - exp_corr[ilat,ilon,indices[ind,1],indices[ind,2],indices[ind,3],indices[ind,4]] <- doQmapRQUANT(exp[NN$imin_lat[ilat],NN$imin_lon[ilon],indices[ind,1],indices[ind,2],indices[ind,3],indices[ind,4]],fobjj,type="linear",wet.day=ifelse(wetonly==TRUE,thres,FALSE)) - } - } - } - } - - exp_corr <- aperm(exp_corr,c(3,4,5,6,1,2)) - return(exp_corr) +QQCorr <- function(exp, wt, qexp, qobs, NN, wetonly=FALSE, thres=1.) { + ## Data needs to be organized for use in qmap package + ## Aim of correction is to have qexp data on the qobs grid, bias adjusted using QQ method + dim_exp <- dim(exp) + dim_qobs <- dim(qobs) # CHECKS TO BE DONE: Nwt, Nq should be consistent between exp and obs + nlat_o <- dim_qobs[3] + nlon_o <- dim_qobs[4] + Nwt <- dim_qobs[1] + Nq <- dim_qobs[2]-2 + exp_corr <- array(dim=c(nlat_o,nlon_o,dim_exp[1:4])) + exp <- aperm(exp,c(5,6,1,2,3,4)) + for (ilat in 1:nlat_o){ + for (ilon in 1:nlon_o){ + qq_slope <- QQSlope(qexp[,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[,,ilat,ilon]) + for (reg in 1:Nwt){ + modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq+2) + fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq+2) + fobjj<-list(par=list(modq=modq, + fitq=fitq, + slope=as.matrix(c(qq_slope$s1[reg],qq_slope$s2[reg]),nrow=2)) + ) + attr(fobjj$par$slope,"dimnames")<-list(c("lower","upper")) + class(fobjj) <- "fitQmapRQUANT" + indices <- which(wt==reg,arr.ind=TRUE) + for (ind in 1:dim(indices)[1]){ + exp_corr[ilat,ilon,indices[ind,1],indices[ind,2],indices[ind,3],indices[ind,4]] <- doQmapRQUANT(exp[NN$imin_lat[ilat],NN$imin_lon[ilon],indices[ind,1],indices[ind,2],indices[ind,3],indices[ind,4]],fobjj,type="linear",wet.day=ifelse(wetonly==TRUE,thres,FALSE)) + } + } + } + } + exp_corr <- aperm(exp_corr,c(3,4,5,6,1,2)) + return(exp_corr) } diff --git a/R/Adamont_Quantiles_S2D.R b/R/Adamont_Quantiles_S2D.R index 179ba34d..04292e20 100644 --- a/R/Adamont_Quantiles_S2D.R +++ b/R/Adamont_Quantiles_S2D.R @@ -26,8 +26,6 @@ #'wt <- sample(1:4,1*10*20*60,replace=T) # Random sample for weather types numbered from 1 to 4 #'dim(wt) <- c(dataset=1,member=10,sdate=20,ftime=60) -##### Input order should be changed... - Quantiles_S2D <- function(exp, var, wt, tslice=NULL, method='ADA', seqq=c(0.005,seq(from=0.01, to=0.99, by=0.01),0.995), wetonly=FALSE) { # diff --git a/R/Adamont_TimeDisaggregate.R b/R/Adamont_TimeDisaggregate.R index a2ab7c22..b6eec30a 100644 --- a/R/Adamont_TimeDisaggregate.R +++ b/R/Adamont_TimeDisaggregate.R @@ -16,140 +16,135 @@ #' #'@import ncdf4, abind +.Disagg <- function(mod, ref, var){ + tres <- length(ref) # Length of subdaily data (24: hourly, 4: 6-hourly, etc.) + D_mod <- array(dim=c(tres)) + if (mean(ref) > 10e-10){ + correction1 <- c("Qair", "SWdown", "prec", "LWdown") + correction2 <- c("uas", "vas") + ## CHANGE TO STANDARD VARIABLE NAMES ABOVE... + if (var %in% correction1){ + a <- mod/mean(ref) + b <- 0 + } else if (var %in% correction2) { + a <- mod/mean(ref) + if (a <= 1) { + b <- 0 + } else { + a <- 1 + b <- mod - mean(ref) + } + } else { # Standard case such as temperature + a <- 1 + b <- mod-mean(ref) + } + } else { # mean(ref) <= 10e-10 + correction1 <- c("LWdown") + correction2 <- c("uas", "vas") + if (var %in% correction2){ + a <- 0 + b <- rep(mod, length(ref)) + } else if (var %in% correction1) { + a <- mod/mean(ref) + b <- 0 + } else { + a <- 0 + b <- 0 + } + } + D_mod <- a*ref + b + return(D_mod) +} -TimeDisaggregate <- function(exp_corr, analog, obs, obs_daily, var, indx, thres) { - - Disagg <- function(mod, ref, var){ - tres <- length(ref) # Length of subdaily data (24: hourly, 4: 6-hourly, etc.) - D_mod <- array(dim=c(tres)) - if (mean(ref) > 10e-10){ - correction1 = c("Qair", "SWdown", "prec", "LWdown") - correction2 = c("uas", "vas") - # CHANGE TO STANDARD VARIABLE NAMES ABOVE... - if (var %in% correction1){ - a <- mod/mean(ref) - b <- 0 - } else if (var %in% correction2) { - a <- mod/mean(ref) - if (a <= 1) { - b = 0 - } else { - a = 1 - b = mod - mean(ref) - } - } else { # Standard case such as temperature - a <- 1 - b <- mod-mean(ref) - } - } else { # mean(ref) <= 10e-10 - correction1 = c("LWdown") - correction2 = c("uas", "vas") - if (var %in% correction2){ - a = 0 - b = rep(mod, length(ref)) - } else if (var %in% correction1) { - a <- mod/mean(ref) - b <- 0 - } else { - a <- 0 - b <- 0 - } - } - D_mod <- a*ref + b - return(D_mod) - } - - SelPrecAnalog <- function(mod, ref, analog5, thres) { - # Function to select which of the 5 possible analogs is chosen - # mod: array with corrected daily precip data at grid point - # ref: corresponding obs array with subdaily precip data at grid point for five analog days - # analog5: 5 analog days - # thres: threshold for precipitation - day_indx <- 0 - repeat { - day_indx <- day_indx+1 - ind <- analog5[day_indx] - if (day_indx == 5 | (sum(is.na(ref[day_indx,]))==0 & ((mod*86400 < thres ) | (mod*86400 >= thres & mean(ref[day_indx,])*86400 >= thres)))) { - break - } - } - #return(day_indx) - return(ind) - } - - dim_exp <- dim(exp_corr) - exp_dt <- array(dim=c(subd=dim(obs)[5],dim_exp[1:6])) - # obs dimensions: c(dataset=1,member_obs=1,sdate,ftime,subd,nlat_o,nlon_o) - obs2 <- aperm(obs,c(5,1,2,3,4,6,7)) # Start with subdaily dimension - # Flatten obs2 array - dim(obs2) <- c(dim(obs2)[1],prod(dim(obs2)[2:5]),dim(obs2)[6],dim(obs2)[7]) - - if (var != "prec"){ - if (is.null(indx)){ - stop("TimeDisaggregate MUST first be run on precipitation to select appropriate analog index indx") - } else { - # RUN Disagg function on data using analog day indx - for (imb in 1:dim(exp_corr)[2]){ - for (sd in 1:dim(exp_corr)[3]){ - for (ltime in 1:dim(exp_corr)[4]){ - for (ilat in 1:dim(exp_corr)[5]){ - for (ilon in 1:dim(exp_corr)[6]){ - exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) - } - } - } - } - } - } - } else { - if (is.null(indx)){ - # Select appropriate analog among 5 possible -> indx array - indx <- array(dim=dim_exp) -# for (imb in 1:dim(exp_corr)[2]){ -# indx[1,imb,,,,] <- mapply(SelPrecAnalog,exp_corr[1,imb,,,,],obs2,analog[1,imb,,,],thres) -# } - for (imb in 1:dim(exp_corr)[2]){ - for (sd in 1:dim(exp_corr)[3]){ - for (ltime in 1:dim(exp_corr)[4]){ - for (ilat in 1:dim(exp_corr)[5]){ - for (ilon in 1:dim(exp_corr)[6]){ - indx[1,imb,sd,ltime,ilat,ilon] <- SelPrecAnalog(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,,ilat,ilon],analog[1,imb,sd,ltime,],thres) - } - } - } - } - } - # RUN Disagg function on data - for (imb in 1:dim(exp_corr)[2]){ - for (sd in 1:dim(exp_corr)[3]){ - for (ltime in 1:dim(exp_corr)[4]){ - for (ilat in 1:dim(exp_corr)[5]){ - for (ilon in 1:dim(exp_corr)[6]){ - exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) - } - } - } - } - } - # HERE should be same logic as above: apply Disagg function on all dimensions except sub-daily. +.SelPrecAnalog <- function(mod, ref, analog5, thres) { + ## Function to select which of the 5 possible analogs is chosen + ## mod: array with corrected daily precip data at grid point + ## ref: corresponding obs array with subdaily precip data at grid point for five analog days + ## analog5: 5 analog days + ## thres: threshold for precipitation + day_indx <- 0 + repeat { + day_indx <- day_indx+1 + ind <- analog5[day_indx] + if (day_indx == 5 | (sum(is.na(ref[day_indx,]))==0 & ((mod*86400 < thres ) | (mod*86400 >= thres & mean(ref[day_indx,])*86400 >= thres)))) { + break + } + } + return(ind) +} - } else { - warning("Using previously defined indx array for TimeDisaggregate: precipitation disaggregation is run again") - for (imb in 1:dim(exp_corr)[2]){ - for (sd in 1:dim(exp_corr)[3]){ - for (ltime in 1:dim(exp_corr)[4]){ - for (ilat in 1:dim(exp_corr)[5]){ - for (ilon in 1:dim(exp_corr)[6]){ - exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) - } - } - } - } - } - } - - } - return(list(disagg=exp_dt, indx=indx)) +TimeDisaggregate <- function(exp_corr, analog, obs, obs_daily, var, indx, thres) { + dim_exp <- dim(exp_corr) + exp_dt <- array(dim=c(subd=dim(obs)[5],dim_exp[1:6])) + ## obs dimensions: c(dataset=1,member_obs=1,sdate,ftime,subd,nlat_o,nlon_o) + obs2 <- aperm(obs,c(5,1,2,3,4,6,7)) # Start with subdaily dimension + ## Flatten obs2 array + dim(obs2) <- c(dim(obs2)[1],prod(dim(obs2)[2:5]),dim(obs2)[6],dim(obs2)[7]) + if (var != "prec"){ + if (is.null(indx)){ + stop("TimeDisaggregate MUST first be run on precipitation to select appropriate analog index indx") + } else { + ## RUN Disagg function on data using analog day indx + for (imb in 1:dim(exp_corr)[2]){ + for (sd in 1:dim(exp_corr)[3]){ + for (ltime in 1:dim(exp_corr)[4]){ + for (ilat in 1:dim(exp_corr)[5]){ + for (ilon in 1:dim(exp_corr)[6]){ + exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) + } + } + } + } + } + } + } else { + if (is.null(indx)){ + ## Select appropriate analog among 5 possible -> indx array + indx <- array(dim=dim_exp) + ##for (imb in 1:dim(exp_corr)[2]){ + ## indx[1,imb,,,,] <- mapply(SelPrecAnalog,exp_corr[1,imb,,,,],obs2,analog[1,imb,,,],thres) + ##} + for (imb in 1:dim(exp_corr)[2]){ + for (sd in 1:dim(exp_corr)[3]){ + for (ltime in 1:dim(exp_corr)[4]){ + for (ilat in 1:dim(exp_corr)[5]){ + for (ilon in 1:dim(exp_corr)[6]){ + indx[1,imb,sd,ltime,ilat,ilon] <- SelPrecAnalog(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,,ilat,ilon],analog[1,imb,sd,ltime,],thres) + } + } + } + } + } + ## RUN Disagg function on data + for (imb in 1:dim(exp_corr)[2]){ + for (sd in 1:dim(exp_corr)[3]){ + for (ltime in 1:dim(exp_corr)[4]){ + for (ilat in 1:dim(exp_corr)[5]){ + for (ilon in 1:dim(exp_corr)[6]){ + exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) + } + } + } + } + } + } else { + warning("Using previously defined indx array for TimeDisaggregate: precipitation disaggregation is run again") + for (imb in 1:dim(exp_corr)[2]){ + for (sd in 1:dim(exp_corr)[3]){ + for (ltime in 1:dim(exp_corr)[4]){ + for (ilat in 1:dim(exp_corr)[5]){ + for (ilon in 1:dim(exp_corr)[6]){ + exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) + } + } + } + } + } + } + } + ##Rearrange dimension order in exp_dt before output + exp_dt <- aperm(exp_dt,c(2,3,4,5,1,6,7)) # Put subdaily dimension back in place + return(list(disagg=exp_dt, indx=indx)) } -- GitLab From b5586cd4f6ea0c3aab60c72b8561ec69b87f1342 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Thu, 15 Oct 2020 19:24:24 +0200 Subject: [PATCH 12/24] Changed Adamont_Analog.R to AdamontAnalog.R using the Apply() function, named dimensions and generalized input/output s2dvcube data --- R/AdamontAnalog.R | 94 +++++++++++++++++++++++++++++++ R/Adamont_Analog.R | 134 --------------------------------------------- R/Adamont_QQCorr.R | 4 +- 3 files changed, 97 insertions(+), 135 deletions(-) create mode 100644 R/AdamontAnalog.R delete mode 100644 R/Adamont_Analog.R diff --git a/R/AdamontAnalog.R b/R/AdamontAnalog.R new file mode 100644 index 00000000..7d830ede --- /dev/null +++ b/R/AdamontAnalog.R @@ -0,0 +1,94 @@ +#' AdamontAnalog finds analogous days in the reference dataset to experiment data based on +#' weather types +#' +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#' +#'@param exp experiment data, can be output from quantile correction +#'@param wt corresponding weather types (same dimensions as exp but lat/lon) +#'@param obs reference data, lat/lon dimensions need to be the same as exp +#'@param wt_obs corresponding weather types (same dimensions as obs but lat/lon) +#'@param method string indicating the method used for analog definition +#' Currently coded are pattcorr: pattern correlation +#' rain1 (for precip patterns): rain occurrence consistency +#' rain01 (for precip patterns): rain occurrence/non occurrence consistency +#'@param thres real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1 +#' +#'@import CSTools, multiApply, ClimProjDiags + + +AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, method = 'pattcorr', thres = NULL,search_obsdims = c('member', 'sdate', 'ftime'), londim = 'lon', latdim = 'lat') { + # exp: lat, lon, sdate, ftime, member + # obs: lat, lon, dims for searching 'sdate' 'ftime'... + # wt_exp: sdate, ftime, member + # wt_obs: the dims for searching + # ADD TESTS: IF METHOD=PATTCORR, THRES=NULL IS FINE, IF NOT, THRES IS NEEDED! + + # Reshaping obs: + ## The dimensions were to search in a single dim + if (length(search_obsdims) > 1) { + for (i in 1:(length(search_obsdims) - 1)) { + obs <- MergeDims(obs, search_obsdims[i:(i + 1)], rename = search_obsdims[i + 1]) + wt_obs <- MergeDims(wt_obs, search_obsdims[i:(i + 1)], rename = search_obsdims[i + 1]) + } + } + names(dim(obs))[which(names(dim(obs)) == search_obsdims[length(search_obsdims)])] <- 'time' + names(dim(wt_obs))[which(names(dim(wt_obs)) == search_obsdims[length(search_obsdims)])] <- 'time' + # Split 'time' dim in weather types + obs <- SplitDim(obs, split_dim = 'time', indices = as.vector(wt_obs)) + names(dim(obs))[which(names(dim(obs)) == 'monthly')] <- 'type' + + analog_vals <- Apply(list(exp, obs, wt_exp), target_dims = list(c(londim, latdim), + c(londim, latdim, 'time', 'type'), NULL), + .analogs, method = method, thres = thres)$output1 + + # Reshaping output: + analog_vals <- Subset(analog_vals,along='type',indices=1,drop='selected') + poslat <- which(names(dim(analog_vals)) == latdim) + poslon <- which(names(dim(analog_vals)) == londim) + postime <- which(names(dim(analog_vals)) == 'time') # Dimension with N analogs + pos <- 1:length(dim(analog_vals)) + if (poslat > poslon){ + analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslon,poslat)) + } else { + analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslat,poslon)) + } + # Renaming 'time' dim to 'analog' + names(dim(analog_vals))[which(names(dim(analog_vals)) == 'time')] <- 'analog' + return(analog_vals) +} + + +.analogs <- function(exp, obs, wt_exp, method = 'pattcorr', thres = NULL, + londimexp = 'lon', latdimexp = 'lat', londimobs = 'lon', latdimobs = 'lat') { + # exp: lon, lat + # obs: lon, lat, time, wt + # wt_exp: wt single scalar + + search_analog <- switch(method, 'rain1' = .rain1, 'rain01' = .rain01, + 'pattcorr' = .pattcor, + stop(paste0("Adamont Analog function supports methods ", + "'rain1', 'rain01', 'pattcorr' only"))) + + obs <- Subset(obs, along = 'type', indices = wt_exp) + accuracy <- Apply(list(exp, obs), target_dims = list(c(londimexp, latdimexp), + c(londimobs, latdimobs)), + search_analog, thres = thres)$output1 + obs <- Subset(obs, along = 'time', indices = order(accuracy, decreasing = TRUE)[1:5]) + return(obs) +} + +.rain1 <- function(exp_day, obs_day, thres) { + accuracy <- sum((obs_day >= thres) * (exp_day >= thres)) + return(accuracy) +} +.rain01 <- function(exp_day, obs_day, thres) { + accuracy <- sum(diag(table((obs_day >= thres),(exp_day >= thres)))) + return(accuracy) +} +.pattcor <- function(exp_day, obs_day, thres = NULL) { + accuracy <- cor(as.vector(obs_day),as.vector(exp_day)) + return(accuracy) +} + + diff --git a/R/Adamont_Analog.R b/R/Adamont_Analog.R deleted file mode 100644 index b3aa8f25..00000000 --- a/R/Adamont_Analog.R +++ /dev/null @@ -1,134 +0,0 @@ -#' ADAMONT Analog finds analogous days in the reference dataset to experiment data based on -#' precipitation patterns -#' -#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version -#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation -#' -#'@param exp_corr output from QQ_Corr, dimensions c(dataset=1,member,sdate,ftime,nlat_o,nlon_o) -#'@param wt corresponding weather types, dimensions c(dataset=1,member,sdate,ftime) -#'@param obs reference data, dimensions c(dataset=1,member_obs=1,sdate,ftime,nlat_o,nlon_o) -#'@param wt_obs corresponding weather types, dimensions c(dataset=1,member_obs=1,sdate,ftime) -#'@param var variable -#'@param method string indicating the method used for analog definition -#' This needs to be detailed more... Options rain1, rain01, patt_corr -#'@param thres real number indicating the precipitation threshold to define a rainy day -#'@param exp_pr, needed if var is not precip -#'@param obs_pr, needed if var is not precip -#' -#'@import ncdf4, abind - - -.MySample <- function(what,howmany){ - if (length(what) > 1){ - ff <- sample(what,howmany) - } else { - ff <- what - } - return(ff) -} - -.SearchAnalogRain1 <- function(exp_day,obs,iwt,thres) { - ## exp_day: exp_day[1,imb,sdate,ftime,,] - ## obs: full obs array (all sdates and ftimes) - ## iwt: indices corresponding to the weather type of exp_day - ## thres: threshold for precipitation (rainy day or not) - analog5 <- array(dim=c(5)) - accuracy <- array(dim=length(iwt)) - for (i in iwt){ - accuracy[which(iwt==i)]<-sum((obs[i,,]>=thres)*(exp_day>=thres)) - } - best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] - vals_accuracy <- unique(best_accuracy) - nunique <- length(vals_accuracy) - count <- array(dim=nunique+1) - count[1] <- 0 - for (ss in 1:nunique){ - count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) - analog5[(count[ss]+1):count[ss+1]] <- MySample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) - } - return(analog5) -} - -.SearchAnalogRain01 <- function(exp_day,obs,iwt,thres) { - ## exp_day: exp_day[1,imb,sdate,ftime,,] - ## obs: full obs array (all sdates and ftimes) - ## iwt: indices corresponding to the weather type of exp_day - ## thres: threshold for precipitation (rainy day or not) - ## Output: analog5 -> 5 indices corresponding to analog days of exp_day - analog5 <- array(dim=c(5)) - accuracy <- array(dim=length(iwt)) - for (i in iwt){ - accuracy[which(iwt==i)]<-sum(diag(table((obs[i,,]>=thres),(exp_day>=thres)))) - } - best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] - vals_accuracy <- unique(best_accuracy) - nunique <- length(vals_accuracy) - count <- array(dim=nunique+1) - count[1] <- 0 - for (ss in 1:nunique){ - count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) - analog5[(count[ss]+1):count[ss+1]] <- MySample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) - } - return(analog5) -} - -.SearchAnalogPattcorr <- function(exp_day,obs,iwt,thres) { - ## exp_day: exp_day[1,imb,sdate,ftime,,] - ## obs: full obs array (all sdates and ftimes) - ## iwt: indices corresponding to the weather type of exp_day - ## thres: threshold for precipitation (rainy day or not) - ## Output: analog5 -> 5 indices corresponding to analog days of exp_day - analog5 <- array(dim=c(5)) - accuracy <- array(dim=length(iwt)) - for (i in iwt){ - accuracy[which(iwt==i)]<-cor(as.vector(obs[i,,]),as.vector(exp_day)) - } - best_accuracy <- sort(accuracy,decreasing=TRUE)[1:5] - vals_accuracy <- unique(best_accuracy) - nunique <- length(vals_accuracy) - count <- array(dim=nunique+1) - count[1] <- 0 - for (ss in 1:nunique){ - count[ss+1] <- count[ss]+sum(best_accuracy==vals_accuracy[ss]) - analog5[(count[ss]+1):count[ss+1]] <- MySample(iwt[which(accuracy == vals_accuracy[ss])],sum(best_accuracy == vals_accuracy[ss])) - } - return(analog5) -} - - -AdamontAnalog <- function(exp_corr, wt, obs, wt_obs, var, method='rain1', thres=1., exp_pr=NULL, obs_pr=NULL) { - ## CHECKS TO BE ADDED - ## Do the dimensions correspond between obs_pr and obs? exp_pr and exp? - ## INCLUDING CHECKS ON DIMENSIONS - ## CHECK ON VALUE OF PRECIPITATION - ## ONE OPTION IS TO GIVE THE NAME OF PRECIP VARIABLE IN THE WRAPPER FUNCTION - if ((var != "prec") & (is.null(exp_pr) | is.null(obs_pr))) { - stop("Analog is selected according to precipitation pattern: exp_pr and obs_pr must be provided") - } else { - exp_pr <- exp_corr - obs_pr <- obs - } - search_analog <- switch(method, 'rain1'=SearchAnalogRain1, 'rain01'=SearchAnalogRain01, 'pattcorr'=SearchAnalogPattcorr, stop("Adamont Analog function supports methods 'rain1', 'rain01', 'pattcorr' only")) - analog <- array(dim=c(1,member=dim(exp_corr)[2],sdate=dim(exp_corr)[3],ftime=dim(exp_corr)[4],ndays=5)) - for (imb in 1:dim(exp_corr)[2]){ - for (isdate in 1:dim(exp_corr)[3]){ - for (itime in 1:dim(exp_corr)[4]){ - iwt <- wt[1,imb,isdate,itime] - index <- which(wt_obs == iwt) - ## No rain day case : pick five random days among corresponding weather type - if (sum(exp_pr[1,imb,isdate,itime,,]>=thres) == 0) { - analog[1,imb,isdate,itime,] <- obs[MySample(index,5)] - } else { - ## Find five closest analogs among corresponding weather type - if (length(index)!=0){ - obs_pr2 <- obs_pr - dim(obs_pr2) <- c(1*1*dim(obs_pr)[3]*dim(obs_pr)[4],dim(obs_pr)[5],dim(obs_pr)[6]) - analog[1,imb,isdate,itime,] <- search_analog(exp_pr[1,imb,isdate,itime,,],obs_pr2,index,thres) - } - } - } - } - } - return(analog) -} - diff --git a/R/Adamont_QQCorr.R b/R/Adamont_QQCorr.R index c6e3eeac..f9d94766 100644 --- a/R/Adamont_QQCorr.R +++ b/R/Adamont_QQCorr.R @@ -10,6 +10,8 @@ #'@param NN list (output from Adamont_NearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o) #' #'@import qmap +#'@import multiApply +#'@import abind .QQSlope <- function(qexp, qobs) { Nwt <- dim(qexp)[1] @@ -32,7 +34,7 @@ QQCorr <- function(exp, wt, qexp, qobs, NN, wetonly=FALSE, thres=1.) { exp <- aperm(exp,c(5,6,1,2,3,4)) for (ilat in 1:nlat_o){ for (ilon in 1:nlon_o){ - qq_slope <- QQSlope(qexp[,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[,,ilat,ilon]) + qq_slope <- .QQSlope(qexp[,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[,,ilat,ilon]) for (reg in 1:Nwt){ modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq+2) fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq+2) -- GitLab From f2880b974371f86ac9628e381d9e224862628318 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Tue, 10 Nov 2020 12:58:11 +0100 Subject: [PATCH 13/24] Adjusting AdamontAnalog to new version of SplitDims function --- R/AdamontAnalog.R | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/R/AdamontAnalog.R b/R/AdamontAnalog.R index 7d830ede..21a8aeb8 100644 --- a/R/AdamontAnalog.R +++ b/R/AdamontAnalog.R @@ -24,8 +24,11 @@ AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, method = 'pattcorr', thres # wt_obs: the dims for searching # ADD TESTS: IF METHOD=PATTCORR, THRES=NULL IS FINE, IF NOT, THRES IS NEEDED! + # Position of lat/lon dimensions in exp data + poslatexp <- which(names(dim(exp)) == latdim) + poslonexp <- which(names(dim(exp)) == londim) # Reshaping obs: - ## The dimensions were to search in a single dim + ## The dimensions where to search in a single dim if (length(search_obsdims) > 1) { for (i in 1:(length(search_obsdims) - 1)) { obs <- MergeDims(obs, search_obsdims[i:(i + 1)], rename = search_obsdims[i + 1]) @@ -35,8 +38,7 @@ AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, method = 'pattcorr', thres names(dim(obs))[which(names(dim(obs)) == search_obsdims[length(search_obsdims)])] <- 'time' names(dim(wt_obs))[which(names(dim(wt_obs)) == search_obsdims[length(search_obsdims)])] <- 'time' # Split 'time' dim in weather types - obs <- SplitDim(obs, split_dim = 'time', indices = as.vector(wt_obs)) - names(dim(obs))[which(names(dim(obs)) == 'monthly')] <- 'type' + obs <- SplitDim(obs, split_dim = 'time', indices = as.vector(wt_obs), new_dim_name='type') analog_vals <- Apply(list(exp, obs, wt_exp), target_dims = list(c(londim, latdim), c(londim, latdim, 'time', 'type'), NULL), @@ -48,7 +50,7 @@ AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, method = 'pattcorr', thres poslon <- which(names(dim(analog_vals)) == londim) postime <- which(names(dim(analog_vals)) == 'time') # Dimension with N analogs pos <- 1:length(dim(analog_vals)) - if (poslat > poslon){ + if (poslatexp > poslonexp){ analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslon,poslat)) } else { analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslat,poslon)) -- GitLab From 48c62ebd78b349626383c256ed1f398b78a90d28 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Fri, 13 Nov 2020 19:05:28 +0100 Subject: [PATCH 14/24] Finalizing AdamontAnalog function, preparing AdamontNearestNeighbors.R and AdamontQQCorr (2nd one still in progress) --- R/AdamontNearestNeighbors.R | 94 +++++++++++++++++ R/{Adamont_QQCorr.R => AdamontQQCorr.R} | 35 ++++++- R/CST_AdamontAnalog.R | 129 ++++++++++++++++++++++++ 3 files changed, 256 insertions(+), 2 deletions(-) create mode 100644 R/AdamontNearestNeighbors.R rename R/{Adamont_QQCorr.R => AdamontQQCorr.R} (64%) create mode 100644 R/CST_AdamontAnalog.R diff --git a/R/AdamontNearestNeighbors.R b/R/AdamontNearestNeighbors.R new file mode 100644 index 00000000..d9b86e30 --- /dev/null +++ b/R/AdamontNearestNeighbors.R @@ -0,0 +1,94 @@ +#' ADAMONT Nearest Neighbors computes the distance between reference data grid centroid and SF data grid +#' +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#'@description This function computes the nearest neighbor for each reference data (lon, lat) point in the experiment dataset by computing the distance between the reference dataset grid and the experiment data. This is the first step in the ADAMONT method adapted from Verfaillie et al. (2018). +#' +#'@param method a string among three options ('ADA': standard ADAMONT distance, 'simple': lon/lat straight euclidian distance, 'radius': distance on the sphere) +#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment longitudes in \code{$lon} and latitudes in \code{$lat} +#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the reference data on a different grid, with longitudes in \code{$lon} and latitudes in \code{$lat}. +#' +#'@return NN +#' +#'@import s2dverification +#'@import ncdf4 + +NearestNeighbors <- function (exp, obs, method='ADA') { + + if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { + stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + exp_lon <- exp$lon + exp_lat <- exp$lat + obs_lon <- obs$lon + obs_lat <- obs$lat + dim_exp_lon <- dim(exp_lon) + dim_exp_lat <- dim(exp_lat) + dim_obs_lon <- dim(obs_lon) + dim_obs_lat <- dim(obs_lat) + + # Flatten longitudes and latitudes in case of 2-D longitudes and latitudes (Lambert grids, etc.) + if ((length(dim_exp_lon)==2) & (length(dim_exp_lat)==2)){ + dim(exp_lon) <- c(dim_exp_lon[1]*dim_exp_lon[2]) + dim(exp_lat) <- c(dim_exp_lat[1]*dim_exp_lat[2]) + } + if ((length(dim_obs_lon)==2) & (length(dim_obs_lat)==2)){ + dim(obs_lon) <- c(dim_obs_lon[1]*dim_obs_lon[2]) + dim(obs_lat) <- c(dim_obs_lat[1]*dim_obs_lat[2]) + } + + # Now lat and lon arrays have 1 dimension, length npt (= nlat*nlon) + + OBS_grid <- cbind(obs_lon,obs_lat) + EXP_grid <- cbind(exp_lon,exp_lat) + + dist_min<-min_lon<-min_lat<-imin_lon<-imin_lat<-array(dim=nrow(OBS_grid)) + if (method == 'ADA'){ + C<-cos(OBS_grid[,2]*pi/180)^2 + for (i in 1:nrow(OBS_grid)){ + dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+C[i]*(OBS_grid[i,1]-EXP_grid[,1])^2 + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else if (method == 'simple'){ + for (i in 1:nrow(OBS_grid)){ + dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+(OBS_grid[i,1]-EXP_grid[,1])^2 + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else if (method == 'radius'){ + R <- 6371e3 # metres, Earth radius + EXP_gridr<-EXP_grid*pi/180 + OBS_gridr<-OBS_grid*pi/180 + for (i in 1:nrow(OBS_grid)){ + a<-sin((OBS_gridr[i,2]-EXP_gridr[,2])/2)^2 + cos(OBS_gridr[i,2])*cos(EXP_gridr[,2])*sin((OBS_gridr[i,1]-EXP_gridr[,1])/2)^2 + c<-2*atan2(sqrt(a),sqrt(1-a)) + dist<-R*c + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else { + stop("Adamont_NearestNeighbors supports method = 'ADA', 'simple' or 'radius' only.") + } + + # Reshape outputs to original grid + dim(min_lon)=dim_obs_lon + dim(min_lat)=dim_obs_lat + dim(imin_lon)=dim_obs_lon + dim(imin_lat)=dim_obs_lat + + NN=list(dist_min=dist_min, min_lon=min_lon, min_lat=min_lat, imin_lon=imin_lon, imin_lat=imin_lat) + + return(NN) +} + diff --git a/R/Adamont_QQCorr.R b/R/AdamontQQCorr.R similarity index 64% rename from R/Adamont_QQCorr.R rename to R/AdamontQQCorr.R index f9d94766..aefc6a4b 100644 --- a/R/Adamont_QQCorr.R +++ b/R/AdamontQQCorr.R @@ -1,4 +1,4 @@ -#' ADAMONT QQCorr computes quantile-quantile correction of seasonal or decadal forecast data +#'ADAMONT QQCorr computes quantile-quantile correction of seasonal or decadal forecast data #' #'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version #'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation @@ -21,7 +21,34 @@ return(list(s1=slope1,s2=slope2)) } -QQCorr <- function(exp, wt, qexp, qobs, NN, wetonly=FALSE, thres=1.) { +AdamontQQCorr <- function(exp, wt, obs, wt_obs, NN, wetonly=FALSE, thres=1.,corrdims = c('member', 'sdate', 'ftime'), londim='lon', latdim='lat') { + + obsdims <- names(dim(obs$data)) + poslat <- which(obsdims == latdim) + poslon <- which(obsdims == londim) + nlat_o <- dim(obs$data)[poslat] + nlon_o <- dim(obs$data)[poslon] + ilat_o <- array(c(1:nlat_o)) + dim(ilat_o) <- c('lat') + ilon_o <- array(c(1:nlon_o)) + dim(ilon_o) <- c('lon') + # First step if obs data is higher resolution than exp data is to use nearest neighbor + # to compute downscaling of exp data + exp_corr <- Apply(list(exp$data,ilat_o,ilon_o),target_dims=list(c(londim,latdim),'lat','lon'),.getNN,NN=NN)$output1 + + # Reorder exp_corr dimensions to match exp dimensions + + + # Both exp and obs data are now on the same grid + # Use CST_QuantileMapping function for quantile mapping, depending on weather type + + # Reshape obs and exp data + # Apply QuantileMapping with exp_cor output for each weather type + + + # Reshaping obs and exp_corr data + + ## Data needs to be organized for use in qmap package ## Aim of correction is to have qexp data on the qobs grid, bias adjusted using QQ method dim_exp <- dim(exp) @@ -55,3 +82,7 @@ QQCorr <- function(exp, wt, qexp, qobs, NN, wetonly=FALSE, thres=1.) { return(exp_corr) } +.getNN <- function(exp,ilat,ilon,NN){ + return(exp[NN$imin_lat[ilat,ilon],NN$imin_lon[ilat,ilon]]) +} + diff --git a/R/CST_AdamontAnalog.R b/R/CST_AdamontAnalog.R new file mode 100644 index 00000000..76578026 --- /dev/null +++ b/R/CST_AdamontAnalog.R @@ -0,0 +1,129 @@ +#'AdamontAnalog finds analogous data in the reference dataset to experiment data based on +#'weather types +#' +#'@description This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data. +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#' +#'@param exp experiment data an object of class \code{s2dv_cube}, can be output from quantile correction +#'@param wt corresponding weather types (same dimensions as \code{exp} but lat/lon), also of class \code{s2dv_cube} +#'@param obs reference data, also of class \code{s2dv_cube}. Note that lat/lon dimensions need to be the same as \code{exp} +#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but lat/lon) +#'@param nanalogs integer defining the number of analog values to return (default: 5) +#'@param method a character string indicating the method used for analog definition +#' Coded are 'pattcorr': pattern correlation +#' 'rain1' (for precip patterns): rain occurrence consistency +#' 'rain01' (for precip patterns): rain occurrence/non occurrence consistency +#'@param thres real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1 +#'@param search_obsdims list of dimensions in \code{obs} along which analogs are searched for +#'@param londim name of longitude dimension +#'@param latdim name of latitude dimension +#'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data +#'@import CSTools, multiApply, ClimProjDiags + +CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, method = 'pattcorr', thres = NULL, search_obsdims = c('member', 'sdate', 'ftime'), londim = 'lon', latdim = 'lat') { + + if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { + stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (!(method %in% c('pattcorr','rain1','rain01'))) { + stop("Input parameter 'method' must be 'pattcorr', 'rain1', or 'rain01'") + } + if (!is.null(nanalogs)){ + nanalogs <- 5 + } + dimnames <- names(dim(obs$data)) + if (!is.null(which(dimnames == latdim)) || !is.null(which(dimnames == londim))) { + stop("'londim' or 'latdim' inputs should correspond to 'obs' dimension names") + } + # More sanity checks needed? Same dimensions in obs and exp; wt_obs and wt_exp + analog_vals <- AdamontAnalog(exp = exp$data, obs = obs$data, wt_exp = wt_exp, + nanalogs = nanalogs, method = method, thres = thres, + search_obsdims = search_obsdims, londim = londim, + latdim = latdim ) + + return(analog_vals) +} + +AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs=5, method = 'pattcorr', thres = NULL,search_obsdims = c('member', 'sdate', 'ftime'), londim = 'lon', latdim = 'lat') { + # exp: lat, lon, sdate, ftime, member + # obs: lat, lon, dims for searching 'sdate' 'ftime'... + # wt_exp: sdate, ftime, member + # wt_obs: the dims for searching + if (method %in% c('rain1','rain01') & is.null(thres)){ + stop("Threshold 'thres' must be defined with methods 'rain1' and 'rain01'") + } + if (method == 'pattcorr' & !is.null(thres)){ + warning("Parameter 'thres' is not used with method 'pattcorr'.") + } + # Position of lat/lon dimensions in exp data + poslatexp <- which(names(dim(exp)) == latdim) + poslonexp <- which(names(dim(exp)) == londim) + # Reshaping obs: + ## The dimensions where to search in a single dim + if (length(search_obsdims) > 1) { + for (i in 1:(length(search_obsdims) - 1)) { + obs <- MergeDims(obs, search_obsdims[i:(i + 1)], rename = search_obsdims[i + 1]) + wt_obs <- MergeDims(wt_obs, search_obsdims[i:(i + 1)], rename = search_obsdims[i + 1]) + } + } + names(dim(obs))[which(names(dim(obs)) == search_obsdims[length(search_obsdims)])] <- 'time' + names(dim(wt_obs))[which(names(dim(wt_obs)) == search_obsdims[length(search_obsdims)])] <- 'time' + # Split 'time' dim in weather types + obs <- SplitDim(obs, split_dim = 'time', indices = as.vector(wt_obs), new_dim_name='type') + + analog_vals <- Apply(list(exp, obs, wt_exp), target_dims = list(c(londim, latdim), + c(londim, latdim, 'time', 'type'), NULL), + .analogs, method = method, thres = thres)$output1 + + # Reshaping output: + analog_vals <- Subset(analog_vals,along='type',indices=1,drop='selected') + poslat <- which(names(dim(analog_vals)) == latdim) + poslon <- which(names(dim(analog_vals)) == londim) + postime <- which(names(dim(analog_vals)) == 'time') # Dimension with N analogs + pos <- 1:length(dim(analog_vals)) + if (poslatexp > poslonexp){ + analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslon,poslat)) + } else { + analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslat,poslon)) + } + # Renaming 'time' dim to 'analog' + names(dim(analog_vals))[which(names(dim(analog_vals)) == 'time')] <- 'analog' + return(analog_vals) +} + + +.analogs <- function(exp, obs, wt_exp, nanalogs = 5, method = 'pattcorr', thres = NULL, + londimexp = 'lon', latdimexp = 'lat', londimobs = 'lon', latdimobs = 'lat') { + # exp: lon, lat + # obs: lon, lat, time, wt + # wt_exp: wt single scalar + + search_analog <- switch(method, 'rain1' = .rain1, 'rain01' = .rain01, + 'pattcorr' = .pattcor, + stop(paste0("Adamont Analog function supports methods ", + "'rain1', 'rain01', 'pattcorr' only"))) + + obs <- Subset(obs, along = 'type', indices = wt_exp) + accuracy <- Apply(list(exp, obs), target_dims = list(c(londimexp, latdimexp), + c(londimobs, latdimobs)), + search_analog, thres = thres)$output1 + obs <- Subset(obs, along = 'time', indices = order(accuracy, decreasing = TRUE)[1:nanalogs]) + return(obs) +} + +.rain1 <- function(exp_day, obs_day, thres) { + accuracy <- sum((obs_day >= thres) * (exp_day >= thres)) + return(accuracy) +} +.rain01 <- function(exp_day, obs_day, thres) { + accuracy <- sum(diag(table((obs_day >= thres),(exp_day >= thres)))) + return(accuracy) +} +.pattcor <- function(exp_day, obs_day, thres = NULL) { + accuracy <- cor(as.vector(obs_day),as.vector(exp_day)) + return(accuracy) +} + + -- GitLab From 9c6347a9974d3a9361f320c83ccaac84d578742b Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Thu, 19 Nov 2020 17:29:12 +0100 Subject: [PATCH 15/24] Work on AdamontQQCorr function to handle extra NAs if needed for QuantileMapping function --- R/AdamontNearestNeighbors.R | 113 +++++++++++++-------------- R/AdamontQQCorr.R | 32 +------- R/CST_AdamontAnalog.R | 28 ++++--- R/CST_AdamontQQCorr.R | 148 ++++++++++++++++++++++++++++++++++++ 4 files changed, 225 insertions(+), 96 deletions(-) create mode 100644 R/CST_AdamontQQCorr.R diff --git a/R/AdamontNearestNeighbors.R b/R/AdamontNearestNeighbors.R index d9b86e30..1873dbb8 100644 --- a/R/AdamontNearestNeighbors.R +++ b/R/AdamontNearestNeighbors.R @@ -27,66 +27,69 @@ NearestNeighbors <- function (exp, obs, method='ADA') { dim_exp_lat <- dim(exp_lat) dim_obs_lon <- dim(obs_lon) dim_obs_lat <- dim(obs_lat) - + # Check if one of the grids is non-regular: + if ((length(dim_exp_lon)==2) || (length(dim_obs_lon)==2)){ # Flatten longitudes and latitudes in case of 2-D longitudes and latitudes (Lambert grids, etc.) - if ((length(dim_exp_lon)==2) & (length(dim_exp_lat)==2)){ - dim(exp_lon) <- c(dim_exp_lon[1]*dim_exp_lon[2]) - dim(exp_lat) <- c(dim_exp_lat[1]*dim_exp_lat[2]) - } - if ((length(dim_obs_lon)==2) & (length(dim_obs_lat)==2)){ - dim(obs_lon) <- c(dim_obs_lon[1]*dim_obs_lon[2]) - dim(obs_lat) <- c(dim_obs_lat[1]*dim_obs_lat[2]) - } - - # Now lat and lon arrays have 1 dimension, length npt (= nlat*nlon) - - OBS_grid <- cbind(obs_lon,obs_lat) - EXP_grid <- cbind(exp_lon,exp_lat) + if ((length(dim_exp_lon)==2) & (length(dim_exp_lat)==2)){ + dim(exp_lon) <- c(dim_exp_lon[1]*dim_exp_lon[2]) + dim(exp_lat) <- c(dim_exp_lat[1]*dim_exp_lat[2]) + } + if ((length(dim_obs_lon)==2) & (length(dim_obs_lat)==2)){ + dim(obs_lon) <- c(dim_obs_lon[1]*dim_obs_lon[2]) + dim(obs_lat) <- c(dim_obs_lat[1]*dim_obs_lat[2]) + } + # Now lat and lon arrays have 1 dimension, length npt (= nlat*nlon) + OBS_grid <- cbind(obs_lon,obs_lat) + EXP_grid <- cbind(exp_lon,exp_lat) + dist_min<-min_lon<-min_lat<-imin_lon<-imin_lat<-array(dim=nrow(OBS_grid)) + if (method == 'ADA'){ + C<-cos(OBS_grid[,2]*pi/180)^2 + for (i in 1:nrow(OBS_grid)){ + dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+C[i]*(OBS_grid[i,1]-EXP_grid[,1])^2 + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else if (method == 'simple'){ + for (i in 1:nrow(OBS_grid)){ + dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+(OBS_grid[i,1]-EXP_grid[,1])^2 + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else if (method == 'radius'){ + R <- 6371e3 # metres, Earth radius + EXP_gridr<-EXP_grid*pi/180 + OBS_gridr<-OBS_grid*pi/180 + for (i in 1:nrow(OBS_grid)){ + a<-sin((OBS_gridr[i,2]-EXP_gridr[,2])/2)^2 + cos(OBS_gridr[i,2])*cos(EXP_gridr[,2])*sin((OBS_gridr[i,1]-EXP_gridr[,1])/2)^2 + c<-2*atan2(sqrt(a),sqrt(1-a)) + dist<-R*c + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else { + stop("AdamontNearestNeighbors supports method = 'ADA', 'simple' or 'radius' only.") + } + + # Reshape outputs to original grid + dim(min_lon)=dim_obs_lon + dim(min_lat)=dim_obs_lat + dim(imin_lon)=dim_obs_lon + dim(imin_lat)=dim_obs_lat - dist_min<-min_lon<-min_lat<-imin_lon<-imin_lat<-array(dim=nrow(OBS_grid)) - if (method == 'ADA'){ - C<-cos(OBS_grid[,2]*pi/180)^2 - for (i in 1:nrow(OBS_grid)){ - dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+C[i]*(OBS_grid[i,1]-EXP_grid[,1])^2 - dist_min[i]<-min(dist) - min_lon[i]<-EXP_grid[which.min(dist),1] - min_lat[i]<-EXP_grid[which.min(dist),2] - imin_lon[i]<-which(exp_lon==min_lon[i]) - imin_lat[i]<-which(exp_lat==min_lat[i]) - } - } else if (method == 'simple'){ - for (i in 1:nrow(OBS_grid)){ - dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+(OBS_grid[i,1]-EXP_grid[,1])^2 - dist_min[i]<-min(dist) - min_lon[i]<-EXP_grid[which.min(dist),1] - min_lat[i]<-EXP_grid[which.min(dist),2] - imin_lon[i]<-which(exp_lon==min_lon[i]) - imin_lat[i]<-which(exp_lat==min_lat[i]) - } - } else if (method == 'radius'){ - R <- 6371e3 # metres, Earth radius - EXP_gridr<-EXP_grid*pi/180 - OBS_gridr<-OBS_grid*pi/180 - for (i in 1:nrow(OBS_grid)){ - a<-sin((OBS_gridr[i,2]-EXP_gridr[,2])/2)^2 + cos(OBS_gridr[i,2])*cos(EXP_gridr[,2])*sin((OBS_gridr[i,1]-EXP_gridr[,1])/2)^2 - c<-2*atan2(sqrt(a),sqrt(1-a)) - dist<-R*c - dist_min[i]<-min(dist) - min_lon[i]<-EXP_grid[which.min(dist),1] - min_lat[i]<-EXP_grid[which.min(dist),2] - imin_lon[i]<-which(exp_lon==min_lon[i]) - imin_lat[i]<-which(exp_lat==min_lat[i]) - } } else { - stop("Adamont_NearestNeighbors supports method = 'ADA', 'simple' or 'radius' only.") + # Regular lon/lat grid case: has been handled by CST_Load() + stop("AdamontNearestNeighbors is meant for non-regular lat/lon grids; use e.g. CST_Load to interpolate exp onto obs grid") } - # Reshape outputs to original grid - dim(min_lon)=dim_obs_lon - dim(min_lat)=dim_obs_lat - dim(imin_lon)=dim_obs_lon - dim(imin_lat)=dim_obs_lat - NN=list(dist_min=dist_min, min_lon=min_lon, min_lat=min_lat, imin_lon=imin_lon, imin_lat=imin_lat) return(NN) diff --git a/R/AdamontQQCorr.R b/R/AdamontQQCorr.R index aefc6a4b..4767d8f2 100644 --- a/R/AdamontQQCorr.R +++ b/R/AdamontQQCorr.R @@ -41,7 +41,7 @@ AdamontQQCorr <- function(exp, wt, obs, wt_obs, NN, wetonly=FALSE, thres=1.,corr # Both exp and obs data are now on the same grid # Use CST_QuantileMapping function for quantile mapping, depending on weather type - + #test_qmap <- CST_QuantileMapping(exp=mod,obs=obs,sample_dims=c('dataset','member','sdate','ftime'),method='RQUANT') # Reshape obs and exp data # Apply QuantileMapping with exp_cor output for each weather type @@ -49,36 +49,6 @@ AdamontQQCorr <- function(exp, wt, obs, wt_obs, NN, wetonly=FALSE, thres=1.,corr # Reshaping obs and exp_corr data - ## Data needs to be organized for use in qmap package - ## Aim of correction is to have qexp data on the qobs grid, bias adjusted using QQ method - dim_exp <- dim(exp) - dim_qobs <- dim(qobs) # CHECKS TO BE DONE: Nwt, Nq should be consistent between exp and obs - nlat_o <- dim_qobs[3] - nlon_o <- dim_qobs[4] - Nwt <- dim_qobs[1] - Nq <- dim_qobs[2]-2 - exp_corr <- array(dim=c(nlat_o,nlon_o,dim_exp[1:4])) - exp <- aperm(exp,c(5,6,1,2,3,4)) - for (ilat in 1:nlat_o){ - for (ilon in 1:nlon_o){ - qq_slope <- .QQSlope(qexp[,,NN$imin_lat[ilat],NN$imin_lon[ilon]],qobs[,,ilat,ilon]) - for (reg in 1:Nwt){ - modq <- as.matrix(qexp[reg,,NN$imin_lat[ilat],NN$imin_lon[ilon]],nrow=Nq+2) - fitq <- as.matrix(qobs[reg,,ilat,ilon],nrow=Nq+2) - fobjj<-list(par=list(modq=modq, - fitq=fitq, - slope=as.matrix(c(qq_slope$s1[reg],qq_slope$s2[reg]),nrow=2)) - ) - attr(fobjj$par$slope,"dimnames")<-list(c("lower","upper")) - class(fobjj) <- "fitQmapRQUANT" - indices <- which(wt==reg,arr.ind=TRUE) - for (ind in 1:dim(indices)[1]){ - exp_corr[ilat,ilon,indices[ind,1],indices[ind,2],indices[ind,3],indices[ind,4]] <- doQmapRQUANT(exp[NN$imin_lat[ilat],NN$imin_lon[ilon],indices[ind,1],indices[ind,2],indices[ind,3],indices[ind,4]],fobjj,type="linear",wet.day=ifelse(wetonly==TRUE,thres,FALSE)) - } - } - } - } - exp_corr <- aperm(exp_corr,c(3,4,5,6,1,2)) return(exp_corr) } diff --git a/R/CST_AdamontAnalog.R b/R/CST_AdamontAnalog.R index 76578026..47f2d381 100644 --- a/R/CST_AdamontAnalog.R +++ b/R/CST_AdamontAnalog.R @@ -5,10 +5,10 @@ #'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version #'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation #' -#'@param exp experiment data an object of class \code{s2dv_cube}, can be output from quantile correction -#'@param wt corresponding weather types (same dimensions as \code{exp} but lat/lon), also of class \code{s2dv_cube} +#'@param exp experiment data an object of class \code{s2dv_cube}, can be output from quantile correction using CST_AdamontQQCorr.R +#'@param wt_exp corresponding weather types (same dimensions as \code{exp$data} but lat/lon) #'@param obs reference data, also of class \code{s2dv_cube}. Note that lat/lon dimensions need to be the same as \code{exp} -#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but lat/lon) +#'@param wt_obs corresponding weather types (same dimensions as \code{obs$data} but lat/lon) #'@param nanalogs integer defining the number of analog values to return (default: 5) #'@param method a character string indicating the method used for analog definition #' Coded are 'pattcorr': pattern correlation @@ -20,7 +20,14 @@ #'@param latdim name of latitude dimension #'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data #'@import CSTools, multiApply, ClimProjDiags - +#'@examples +#'\dontrun{ +#'wt_exp <- sample(1:3, 15*6*3, replace=T) +#'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) +#'wt_obs <- sample(1:3, 6*3, replace=T) +#'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) +# analog_vals <- CST_AdamontAnalog(exp=lonlat_data$exp, obs=lonlat_data$obs, wt_exp=wt_exp, wt_obs=wt_obs, nanalogs=2) +#'} CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, method = 'pattcorr', thres = NULL, search_obsdims = c('member', 'sdate', 'ftime'), londim = 'lon', latdim = 'lat') { if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { @@ -30,18 +37,19 @@ CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, method = 'patt if (!(method %in% c('pattcorr','rain1','rain01'))) { stop("Input parameter 'method' must be 'pattcorr', 'rain1', or 'rain01'") } - if (!is.null(nanalogs)){ + if (is.null(nanalogs)){ nanalogs <- 5 } dimnames <- names(dim(obs$data)) - if (!is.null(which(dimnames == latdim)) || !is.null(which(dimnames == londim))) { + if (is.null(which(dimnames == latdim)) || is.null(which(dimnames == londim))) { stop("'londim' or 'latdim' inputs should correspond to 'obs' dimension names") } # More sanity checks needed? Same dimensions in obs and exp; wt_obs and wt_exp - analog_vals <- AdamontAnalog(exp = exp$data, obs = obs$data, wt_exp = wt_exp, - nanalogs = nanalogs, method = method, thres = thres, - search_obsdims = search_obsdims, londim = londim, - latdim = latdim ) + analog_vals <- AdamontAnalog(exp = exp$data, obs = obs$data, wt_exp = wt_exp, + wt_obs=wt_obs, nanalogs = nanalogs, + method = method, thres = thres, + search_obsdims = search_obsdims, londim = londim, + latdim = latdim ) return(analog_vals) } diff --git a/R/CST_AdamontQQCorr.R b/R/CST_AdamontQQCorr.R new file mode 100644 index 00000000..1916d9a5 --- /dev/null +++ b/R/CST_AdamontQQCorr.R @@ -0,0 +1,148 @@ +#'ADAMONT QQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types +#' +#'@description This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference \code{obs}, typically provided by reanalysis data. +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#' +#'@param exp experiment data an object of class \code{s2dv_cube} +#'@param wt_exp corresponding weather types (same dimensions as \code{exp$data} but lat/lon) +#'@param obs reference data, also of class \code{s2dv_cube}. lat/lon dimensions can differ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and AdamontNearestNeighbors \code{NN} output should be provided +#'@param corrdims list of dimensions in \code{exp} for which quantile mapping correction is applied +#'@param londim character name of longitude dimension in \code{exp} and \code{obs} +#'@param latdim character name of latitude dimension in \code{exp} and \code{obs} +#'@param regrid boolean indicating whether Nearest Neighbors regridding is needed +#'@param NN list (output from Adamont_NearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o) +#' +#'@import qmap +#'@import multiApply +#'@import abind +#'@examples +#'\dontrun{ +#'wt_exp <- sample(1:3, 15*6*3, replace=T) +#'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) +#'wt_obs <- sample(1:3, 6*3, replace=T) +#'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) +#'exp_corr <- CST_AdamontQQCorr(exp=lonlat_data$exp, wt_exp=wt_exp, obs=lonlat_data$obs, wt_obs=wt_obs, corrdims = c('dataset','member','sdate','ftime')) +#'} +CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member','sdate','ftime'), londim='lon', latdim='lat') { + + if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')){ + stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + dimnames <- names(dim(obs$data)) + if (is.null(which(dimnames == latdim)) || is.null(which(dimnames==londim))){ + stop("'londim' or 'latdim' input doesn't match with 'obs' dimension names") + } + if ((length(dim(exp$lon))==2) || (length(dim(obs$lon))==2)){ + myNN <- NearestNeighbors(exp=exp, obs=obs, method='ADA') + exp_corr <- AdamontQQCorr(exp=exp$data, wt_exp=wt_exp, obs=obs$data, wt_obs=wt_obs, corrdims=corrdims, londim=londim, latdim=latdim, regrid=TRUE, NN=myNN) + } else { + # If not (standard case) + exp_corr <- AdamontQQCorr(exp=exp$data, wt_exp=wt_exp, obs=obs$data, wt_obs=wt_obs, corrdims=corrdims, londim=londim, latdim=latdim, regrid=FALSE) + } + return(exp_corr) +} + +AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member', 'sdate', 'ftime'), londim='lon', latdim='lat', regrid=FALSE, NN=NULL) { + + # The regridding part should only be done if lat/lon dimensions of obs and + # exp differ. + # INCLUDE SANITY CHECKS: LATDIM AND LONDIM MUST ALSO MATCH EXP LAT/LON DIM NAMES + # IF REGRID = TRUE AND NN = NULL: error !! + if (regrid == 'TRUE'){ + obsdims <- names(dim(obs)) + poslat <- which(obsdims == latdim) + poslon <- which(obsdims == londim) + nlat_o <- dim(obs)[poslat] + nlon_o <- dim(obs)[poslon] + ilat_o <- array(c(1:nlat_o)) + names(dim(ilat_o))[1] <- latdim + ilon_o <- array(c(1:nlon_o)) + names(dim(ilon_o))[1] <- londim + # First step if obs data is higher resolution than exp data is to use nearest neighbor + # to compute downscaling of exp data + exp_corr <- Apply(list(exp,ilat_o,ilon_o),target_dims=list(c(latdim,londim),latdim,londim),.getNN,NN=NN)$output1 + + # Reorder exp_corr dimensions to match exp dimensions + dexpc <- match(names(dim(exp)), names(dim(exp_corr))) + exp_corr <- aperm(exp_corr,dexpc) + dimnames(exp_corr) <- dimnames(exp)[dexpc] + # Keep original wt_exp for remapping data + wt_exp2 <- wt_exp + + # Both exp and obs data are now on the same grid + } else { + exp_corr <- exp + # Keep original wt_exp for remapping data + wt_exp2 <- wt_exp + } + + # Use CST_QuantileMapping function for quantile mapping, depending on weather type + for (i in 1:(length(corrdims) - 1)) { + obs <- MergeDims(obs, corrdims[i:(i+1)], rename=corrdims[i+1]) + wt_obs <- MergeDims(wt_obs, corrdims[i:(i+1)], rename=corrdims[i+1]) + exp_corr <- MergeDims(exp_corr, corrdims[i:(i+1)], rename=corrdims[i+1]) + wt_exp2 <- MergeDims(wt_exp2, corrdims[i:(i+1)], rename=corrdims[i+1]) + } + names(dim(obs))[which(names(dim(obs)) == corrdims[length(corrdims)])] <- 'time' + names(dim(wt_obs))[which(names(dim(wt_obs)) == corrdims[length(corrdims)])] <- 'time' + names(dim(exp_corr))[which(names(dim(exp_corr)) == corrdims[length(corrdims)])] <- 'time' + names(dim(wt_exp2))[which(names(dim(wt_exp2)) == corrdims[length(corrdims)])] <- 'time' + # Split 'time' dim in weather types + obs <- SplitDim(obs, split_dim='time',indices=as.vector(wt_obs), new_dim_name='type') + exp_corr <- SplitDim(exp_corr, split_dim='time',indices=as.vector(wt_exp2), new_dim_name='type') + # Add NAs to exp_corr if needed to have compatible sample dimensions + numtobs <- dim(obs)[which(names(dim(obs))=='time')] + numtexp <- dim(exp_corr)[which(names(dim(exp_corr))=='time')] + if (numtexp%%numtobs > 0){ + # Create extra dimension and include NAs + ndimexp <- names(dim(exp_corr)) + ndimobs <- names(dim(obs)) + postime <- which(ndimexp=='time') + dimadd <- dim(exp_corr) + dimadd[postime] <- ceiling(numtexp/numtobs)*numtobs-numtexp + exp_corr <- abind::abind(exp_corr,array(NA,dimadd),along=postime) + names(dim(exp_corr)) <- ndimexp + exp_corr <- SplitDim(exp_corr,'time',freq=numtobs,indices=NULL) + dimobs <- c(dim(obs),1) + dim(obs) <- dimobs + names(dim(obs)) <- c(ndimobs,'index') + res <- QuantileMapping(exp=exp_corr,obs=obs,sample_dims=c('time','index'),method='RQUANT') + res <- MergeDims(res,c('time','index')) + # Remove the extra NA values added previously + res <- Subset(res,along='time',indices=1:numtexp) + } else { + # Apply QuantileMapping to exp_corr depending on weather type + res <- QuantileMapping(exp=exp_corr,obs=obs,sample_dims='time',samplemethod='RQUANT') + } + rm(exp_corr) # Save space in memory + # Reshape exp_corr data onto time dimension before 'Split' + rep_pos <- array(NA,c(time=length(wt_exp2))) + pos_time <- which(names(dim(res)) == 'time') + pos_type <- which(names(dim(res)) == 'type') + for (x in unique(wt_exp2)){ + rep_pos[which(wt_exp2==x)]<-1:length(which(wt_exp2==x)) + } + exp_corr <- .unsplit_wtype(exp=res,wt_exp=wt_exp2,rep_pos=rep_pos,pos_time=pos_time) + # Now reshape exp_corr data onto original dimensions + dim(exp_corr) <- c(dim(wt_exp), dim(exp_corr)[-c(pos_time,pos_type)]) + return(exp_corr) +} + +.getNN <- function(exp,ilat,ilon,NN){ + return(exp[NN$imin_lat[ilat,ilon],NN$imin_lon[ilat,ilon]]) +} + +.unsplit_wtype <- function(exp=exp,dim_wt='type',wt_exp=wt_exp,dim_time='time',rep_pos=rep_pos,pos_time=1){ + # Initiate output + new <- Subset(Subset(exp, along=dim_wt, indices=wt_exp[1]), along=dim_time, indices=rep_pos[1]) + dimnames <- names(dim(new)) + for (x in 2:length(wt_exp)){ + dat <- Subset(Subset(exp, along=dim_wt, indices=wt_exp[x]), along=dim_time, indices=rep_pos[x]) + new <- abind::abind(new,dat,along=pos_time) + } + names(dim(new)) <- dimnames + return(new) +} + -- GitLab From b077ab69a14cffef95c49673b554fdf2f60b1290 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Thu, 19 Nov 2020 17:30:18 +0100 Subject: [PATCH 16/24] Remove intermediate (unneeded functions) --- R/AdamontAnalog.R | 96 ---------------------- R/AdamontQQCorr.R | 58 -------------- R/Adamont_NearestNeighbors.R | 88 -------------------- R/Adamont_Quantiles_S2D.R | 87 -------------------- R/Adamont_TimeDisaggregate.R | 150 ----------------------------------- 5 files changed, 479 deletions(-) delete mode 100644 R/AdamontAnalog.R delete mode 100644 R/AdamontQQCorr.R delete mode 100644 R/Adamont_NearestNeighbors.R delete mode 100644 R/Adamont_Quantiles_S2D.R delete mode 100644 R/Adamont_TimeDisaggregate.R diff --git a/R/AdamontAnalog.R b/R/AdamontAnalog.R deleted file mode 100644 index 21a8aeb8..00000000 --- a/R/AdamontAnalog.R +++ /dev/null @@ -1,96 +0,0 @@ -#' AdamontAnalog finds analogous days in the reference dataset to experiment data based on -#' weather types -#' -#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version -#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation -#' -#'@param exp experiment data, can be output from quantile correction -#'@param wt corresponding weather types (same dimensions as exp but lat/lon) -#'@param obs reference data, lat/lon dimensions need to be the same as exp -#'@param wt_obs corresponding weather types (same dimensions as obs but lat/lon) -#'@param method string indicating the method used for analog definition -#' Currently coded are pattcorr: pattern correlation -#' rain1 (for precip patterns): rain occurrence consistency -#' rain01 (for precip patterns): rain occurrence/non occurrence consistency -#'@param thres real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1 -#' -#'@import CSTools, multiApply, ClimProjDiags - - -AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, method = 'pattcorr', thres = NULL,search_obsdims = c('member', 'sdate', 'ftime'), londim = 'lon', latdim = 'lat') { - # exp: lat, lon, sdate, ftime, member - # obs: lat, lon, dims for searching 'sdate' 'ftime'... - # wt_exp: sdate, ftime, member - # wt_obs: the dims for searching - # ADD TESTS: IF METHOD=PATTCORR, THRES=NULL IS FINE, IF NOT, THRES IS NEEDED! - - # Position of lat/lon dimensions in exp data - poslatexp <- which(names(dim(exp)) == latdim) - poslonexp <- which(names(dim(exp)) == londim) - # Reshaping obs: - ## The dimensions where to search in a single dim - if (length(search_obsdims) > 1) { - for (i in 1:(length(search_obsdims) - 1)) { - obs <- MergeDims(obs, search_obsdims[i:(i + 1)], rename = search_obsdims[i + 1]) - wt_obs <- MergeDims(wt_obs, search_obsdims[i:(i + 1)], rename = search_obsdims[i + 1]) - } - } - names(dim(obs))[which(names(dim(obs)) == search_obsdims[length(search_obsdims)])] <- 'time' - names(dim(wt_obs))[which(names(dim(wt_obs)) == search_obsdims[length(search_obsdims)])] <- 'time' - # Split 'time' dim in weather types - obs <- SplitDim(obs, split_dim = 'time', indices = as.vector(wt_obs), new_dim_name='type') - - analog_vals <- Apply(list(exp, obs, wt_exp), target_dims = list(c(londim, latdim), - c(londim, latdim, 'time', 'type'), NULL), - .analogs, method = method, thres = thres)$output1 - - # Reshaping output: - analog_vals <- Subset(analog_vals,along='type',indices=1,drop='selected') - poslat <- which(names(dim(analog_vals)) == latdim) - poslon <- which(names(dim(analog_vals)) == londim) - postime <- which(names(dim(analog_vals)) == 'time') # Dimension with N analogs - pos <- 1:length(dim(analog_vals)) - if (poslatexp > poslonexp){ - analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslon,poslat)) - } else { - analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslat,poslon)) - } - # Renaming 'time' dim to 'analog' - names(dim(analog_vals))[which(names(dim(analog_vals)) == 'time')] <- 'analog' - return(analog_vals) -} - - -.analogs <- function(exp, obs, wt_exp, method = 'pattcorr', thres = NULL, - londimexp = 'lon', latdimexp = 'lat', londimobs = 'lon', latdimobs = 'lat') { - # exp: lon, lat - # obs: lon, lat, time, wt - # wt_exp: wt single scalar - - search_analog <- switch(method, 'rain1' = .rain1, 'rain01' = .rain01, - 'pattcorr' = .pattcor, - stop(paste0("Adamont Analog function supports methods ", - "'rain1', 'rain01', 'pattcorr' only"))) - - obs <- Subset(obs, along = 'type', indices = wt_exp) - accuracy <- Apply(list(exp, obs), target_dims = list(c(londimexp, latdimexp), - c(londimobs, latdimobs)), - search_analog, thres = thres)$output1 - obs <- Subset(obs, along = 'time', indices = order(accuracy, decreasing = TRUE)[1:5]) - return(obs) -} - -.rain1 <- function(exp_day, obs_day, thres) { - accuracy <- sum((obs_day >= thres) * (exp_day >= thres)) - return(accuracy) -} -.rain01 <- function(exp_day, obs_day, thres) { - accuracy <- sum(diag(table((obs_day >= thres),(exp_day >= thres)))) - return(accuracy) -} -.pattcor <- function(exp_day, obs_day, thres = NULL) { - accuracy <- cor(as.vector(obs_day),as.vector(exp_day)) - return(accuracy) -} - - diff --git a/R/AdamontQQCorr.R b/R/AdamontQQCorr.R deleted file mode 100644 index 4767d8f2..00000000 --- a/R/AdamontQQCorr.R +++ /dev/null @@ -1,58 +0,0 @@ -#'ADAMONT QQCorr computes quantile-quantile correction of seasonal or decadal forecast data -#' -#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version -#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation -#' -#'@param exp dimensions c(dataset,member,sdate,ftime,lat,lon) -#'@param wt dimensions c(dataset,member,sdate,ftime) -#'@param qexp dimensions c(Nwt, Nq+2, nlat, nlon) -#'@param qobs dimensions c(Nwt, Nq+2, nlat_o, nlon_o) -#'@param NN list (output from Adamont_NearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o) -#' -#'@import qmap -#'@import multiApply -#'@import abind - -.QQSlope <- function(qexp, qobs) { - Nwt <- dim(qexp)[1] - Nq <- dim(qexp)[2] - 2 - slope1 <- (qobs[,2]-qobs[,1])/(qexp[,2]-qexp[,1]) - slope2 <- matrix(1,nrow=Nwt,ncol=Nq+2) - return(list(s1=slope1,s2=slope2)) -} - -AdamontQQCorr <- function(exp, wt, obs, wt_obs, NN, wetonly=FALSE, thres=1.,corrdims = c('member', 'sdate', 'ftime'), londim='lon', latdim='lat') { - - obsdims <- names(dim(obs$data)) - poslat <- which(obsdims == latdim) - poslon <- which(obsdims == londim) - nlat_o <- dim(obs$data)[poslat] - nlon_o <- dim(obs$data)[poslon] - ilat_o <- array(c(1:nlat_o)) - dim(ilat_o) <- c('lat') - ilon_o <- array(c(1:nlon_o)) - dim(ilon_o) <- c('lon') - # First step if obs data is higher resolution than exp data is to use nearest neighbor - # to compute downscaling of exp data - exp_corr <- Apply(list(exp$data,ilat_o,ilon_o),target_dims=list(c(londim,latdim),'lat','lon'),.getNN,NN=NN)$output1 - - # Reorder exp_corr dimensions to match exp dimensions - - - # Both exp and obs data are now on the same grid - # Use CST_QuantileMapping function for quantile mapping, depending on weather type - #test_qmap <- CST_QuantileMapping(exp=mod,obs=obs,sample_dims=c('dataset','member','sdate','ftime'),method='RQUANT') - # Reshape obs and exp data - # Apply QuantileMapping with exp_cor output for each weather type - - - # Reshaping obs and exp_corr data - - - return(exp_corr) -} - -.getNN <- function(exp,ilat,ilon,NN){ - return(exp[NN$imin_lat[ilat,ilon],NN$imin_lon[ilat,ilon]]) -} - diff --git a/R/Adamont_NearestNeighbors.R b/R/Adamont_NearestNeighbors.R deleted file mode 100644 index 8fc3f49f..00000000 --- a/R/Adamont_NearestNeighbors.R +++ /dev/null @@ -1,88 +0,0 @@ -#' ADAMONT Nearest Neighbors computes the distance between reference data grid centroid and SF data grid -#' -#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version -#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation -#'@description This function computes the nearest neighbor for each reference data (lon, lat) point in the experiment dataset by computing the distance between the reference dataset grid and the experiment data. This is the first step in the ADAMONT method adapted from Verfaillie et al. (2018). -#' -#'@param method a string among three options ('ADA': standard ADAMONT distance, 'simple': lon/lat straight euclidian distance, 'radius': distance on the sphere) -#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment longitudes in \code{$lon} and latitudes in \code{$lat} -#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the reference data on a different grid, with longitudes in \code{$lon} and latitudes in \code{$lat}. -#' -#'@return NN -#' -#'@import s2dverification -#'@import ncdf4 - -NearestNeighbors <- function (method='ADA', exp, obs) { - - exp_lon <- exp$lon - exp_lat <- exp$lat - dim_exp_lon <- dim(exp_lon) - dim_exp_lat <- dim(exp_lat) - obs_lon <- obs$lon - obs_lat <- obs$lat - dim_obs_lon <- dim(obs_lon) - dim_obs_lat <- dim(obs_lat) - - # Flatten longitudes and latitudes in case of 2-D longitudes and latitudes (Lambert grids, etc.) - if ((length(dim_exp_lon)==2) & (length(dim_exp_lat)==2)){ - dim(exp_lon) <- c(dim_exp_lon[1]*dim_exp_lon[2]) - dim(exp_lat) <- c(dim_exp_lat[1]*dim_exp_lat[2]) - } - if ((length(dim_obs_lon)==2) & (length(dim_obs_lat)==2)){ - dim(obs_lon) <- c(dim_obs_lon[1]*dim_obs_lon[2]) - dim(obs_lat) <- c(dim_obs_lat[1]*dim_obs_lat[2]) - } - - # Now lat and lon arrays have 1 dimension, length npt (= nlat*nlon) - - OBS_grid <- cbind(obs_lon,obs_lat) - EXP_grid <- cbind(exp_lon,exp_lat) - - dist_min<-min_lon<-min_lat<-imin_lon<-imin_lat<-array(dim=nrow(OBS_grid)) - if (method == 'ADA'){ - C<-cos(OBS_grid[,2]*pi/180)^2 - for (i in 1:nrow(OBS_grid)){ - dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+C[i]*(OBS_grid[i,1]-EXP_grid[,1])^2 - dist_min[i]<-min(dist) - min_lon[i]<-EXP_grid[which.min(dist),1] - min_lat[i]<-EXP_grid[which.min(dist),2] - imin_lon[i]<-which(exp_lon==min_lon[i]) - imin_lat[i]<-which(exp_lat==min_lat[i]) - } - } else if (method == 'simple'){ - for (i in 1:nrow(OBS_grid)){ - dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+(OBS_grid[i,1]-EXP_grid[,1])^2 - dist_min[i]<-min(dist) - min_lon[i]<-EXP_grid[which.min(dist),1] - min_lat[i]<-EXP_grid[which.min(dist),2] - imin_lon[i]<-which(exp_lon==min_lon[i]) - imin_lat[i]<-which(exp_lat==min_lat[i]) - } - } else if (method == 'radius'){ - R <- 6371e3 # metres, Earth radius - EXP_gridr<-EXP_grid*pi/180 - OBS_gridr<-OBS_grid*pi/180 - for (i in 1:nrow(OBS_grid)){ - a<-sin((OBS_gridr[i,2]-EXP_gridr[,2])/2)^2 + cos(OBS_gridr[i,2])*cos(EXP_gridr[,2])*sin((OBS_gridr[i,1]-EXP_gridr[,1])/2)^2 - c<-2*atan2(sqrt(a),sqrt(1-a)) - dist<-R*c - dist_min[i]<-min(dist) - min_lon[i]<-EXP_grid[which.min(dist),1] - min_lat[i]<-EXP_grid[which.min(dist),2] - imin_lon[i]<-which(exp_lon==min_lon[i]) - imin_lat[i]<-which(exp_lat==min_lat[i]) - } - } else { - stop("Adamont_NearestNeighbors supports method = 'ADA', 'simple' or 'radius' only.") - } - - # Reshape outputs to original grid?? Not really necessary? - #dim(min_lon)=dim_obs_lon - #dim(min_lat)=dim_obs_lat - - NN=list(dist_min=dist_min, min_lon=min_lon, min_lat=min_lat, imin_lon=imin_lon, imin_lat=imin_lat) - - return(NN) -} - diff --git a/R/Adamont_Quantiles_S2D.R b/R/Adamont_Quantiles_S2D.R deleted file mode 100644 index 04292e20..00000000 --- a/R/Adamont_Quantiles_S2D.R +++ /dev/null @@ -1,87 +0,0 @@ -#' ADAMONT Quantiles S2D computes quantiles for seasonal or decadal forecast data -#' -#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version -#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation -#'@description This function computes quantiles on either reference or seasonal-to-decadal ensemble DAILY forecast data following the ADAMONT method adapted from Verfaillie et al. (2018). A prerequisite is to attribute each day and member of the S2D re-forecast to a given weather type / regime. -#' -#'@param exp an object of class \code{s2dv_cube} of lonlat data as returned by \code{CST_Load} function, containing the reference or seasonal forecast experiment data in \code{$data}, dimensions c(dataset,member,sdate,ftime,lat,lon) -#'@param var variable name (string) -#'@param wt an object of same member, sdate and ftime dimensions as exp, with weather type attribution (values 1 to Nwt where Nwt is the number of weather types); can be $cluster output of \code{CST_WeatherRegimes} -#'@param tslice list of c(tmin,tmax) is the range of forecast times for which quantiles should be computed (default: NULL, values computed over all forecast times) -#'@param method a string among two options for extrapolation of extrema (default 'ADA': ADAMONT quantile computation, 'interp': interpolated quantile values) -#'@param seqq a list of quantiles (default to c(0.005, seq(from=0.01, to=0.00, by=0.01), 0.995)) -#'@param wetonly a logical (default to FALSE); which days quantiles of precipitation are computed (all days or wet) -#' -#'@return qexp an array with dimensions c(Nwt, Nq+2, lat, lon) where the Nq+2 values are [min value; Nq quantiles; max value]. -#' -#'@import s2dverification -#'@import ncdf4 -#'@import abind -#' -#' EXAMPLE TO BE INSERTED HERE... -#' -#' -#'exp <- runif(n=1*10*20*60*6*7,min=0.,max=5.) -#'dim(exp) <- c(dataset=1,member=10,sdate=20,ftime=60,lat=6,lon=7) -#'wt <- sample(1:4,1*10*20*60,replace=T) # Random sample for weather types numbered from 1 to 4 -#'dim(wt) <- c(dataset=1,member=10,sdate=20,ftime=60) - -Quantiles_S2D <- function(exp, var, wt, tslice=NULL, method='ADA', seqq=c(0.005,seq(from=0.01, to=0.99, by=0.01),0.995), wetonly=FALSE) { - -# -# -# CHECKS NEEDED: -# WT VALUES SHOULD BE 1:NWT -# DIMENSIONS OF EXP AND WT DATA SHOULD BE CONSISTENT -# - -interp_quant<-function(data){ - out<-quantile(data,p=seqq,na.rm=TRUE) - return(c(min(data,na.rm=TRUE),out,max(data,na.rm=TRUE))) -} - -ADA_quant<-function(data){ - out<-quantile(data,p=seqq,type=1,na.rm=TRUE) - return(c(min(data,na.rm=TRUE),out,max(data,na.rm=TRUE))) -} - -# Which quantiles method to apply -quant <- switch(method,'ADA'=ADA_quant,'interp'=interp_quant) - -# Select ftimes for quantile computation -if (is.null(tslice)){ - warning('tslice=NULL: Computing quantiles over each forecast time') - dat <- exp$data -}else if (length(tslice) != 2) { - stop("Parameter 'tslice' must have length equal to 2 (c(tmin,tmax)) ") -}else { - dat <- exp$data[,,,tslice[1]:tslice[2],,,drop=FALSE] - wt <- wt[,,,tslice[1]:tslice[2],drop=FALSE] -} - -# Number of different weather types -wt_f <- wt -dim(wt_f) <- c(prod(dim(wt))) -Nwt <- length(unique(wt_f)) - -# Wet days only if variable is precipitation/snow and wetonly=TRUE -if ((var=='prlr') && (wetonly==TRUE)){ - norain=which(exp<1.) - dat[norain]=NA -} - -# Rearrange and flatten dimensions to use apply() function -dim_dat <- dim(dat) -dat_f <- aperm(dat,c(5,6,1,2,3,4)) -dim(dat_f) <- c(dim_dat[5],dim_dat[6],prod(dim_dat[1:4])) - -qexp <- array(dim=c(type=Nwt,quantiles=dim(array(seqq))+2,dim_dat[5],dim_dat[6])) - -for (reg in 1:Nwt){ - qexp[reg,,,] <- apply(dat_f[,,which(wt == reg)],c(1,2),quant) -} - -return(qexp) -} - - diff --git a/R/Adamont_TimeDisaggregate.R b/R/Adamont_TimeDisaggregate.R deleted file mode 100644 index b6eec30a..00000000 --- a/R/Adamont_TimeDisaggregate.R +++ /dev/null @@ -1,150 +0,0 @@ -#' ADAMONT TimeDisaggregate uses analogous days from ADAMONT Analog and higher frequency time series in the reference dataset to experiment data based on -#' precipitation patterns -#' -#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version -#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation -#' -#'@param exp_corr output from QQ_Corr, dimensions c(dataset=1,member,sdate,ftime,nlat_o,nlon_o) -#'@param analog output from Analog, dimensions c(dataset=1,member,sdate,ftime,5) (5 analogs) -#'@param obs reference subdaily data, dimensions c(dataset=1,member_obs=1,sdate,ftime,subd,nlat_o,nlon_o) -#' subd is number of subdaily time steps in obs data (24 for hourly, 4 for 6-hourly, etc.) -#'@param obs_daily reference daily data, dimensions c(dataset=1,member_obs=1,sdate,ftime,nlat_o,nlon_o) -#'@param var variable name -#'@param indx index of analog day used for precip correction, dimensions c(dataset=1,member,sdate,ftime) -#' (if var=prec, NULL required, if not, should be provided) -#'@param thres precipitation threshold for definition of rainy day -#' -#'@import ncdf4, abind - -.Disagg <- function(mod, ref, var){ - tres <- length(ref) # Length of subdaily data (24: hourly, 4: 6-hourly, etc.) - D_mod <- array(dim=c(tres)) - if (mean(ref) > 10e-10){ - correction1 <- c("Qair", "SWdown", "prec", "LWdown") - correction2 <- c("uas", "vas") - ## CHANGE TO STANDARD VARIABLE NAMES ABOVE... - if (var %in% correction1){ - a <- mod/mean(ref) - b <- 0 - } else if (var %in% correction2) { - a <- mod/mean(ref) - if (a <= 1) { - b <- 0 - } else { - a <- 1 - b <- mod - mean(ref) - } - } else { # Standard case such as temperature - a <- 1 - b <- mod-mean(ref) - } - } else { # mean(ref) <= 10e-10 - correction1 <- c("LWdown") - correction2 <- c("uas", "vas") - if (var %in% correction2){ - a <- 0 - b <- rep(mod, length(ref)) - } else if (var %in% correction1) { - a <- mod/mean(ref) - b <- 0 - } else { - a <- 0 - b <- 0 - } - } - D_mod <- a*ref + b - return(D_mod) -} - -.SelPrecAnalog <- function(mod, ref, analog5, thres) { - ## Function to select which of the 5 possible analogs is chosen - ## mod: array with corrected daily precip data at grid point - ## ref: corresponding obs array with subdaily precip data at grid point for five analog days - ## analog5: 5 analog days - ## thres: threshold for precipitation - day_indx <- 0 - repeat { - day_indx <- day_indx+1 - ind <- analog5[day_indx] - if (day_indx == 5 | (sum(is.na(ref[day_indx,]))==0 & ((mod*86400 < thres ) | (mod*86400 >= thres & mean(ref[day_indx,])*86400 >= thres)))) { - break - } - } - return(ind) -} - - -TimeDisaggregate <- function(exp_corr, analog, obs, obs_daily, var, indx, thres) { - dim_exp <- dim(exp_corr) - exp_dt <- array(dim=c(subd=dim(obs)[5],dim_exp[1:6])) - ## obs dimensions: c(dataset=1,member_obs=1,sdate,ftime,subd,nlat_o,nlon_o) - obs2 <- aperm(obs,c(5,1,2,3,4,6,7)) # Start with subdaily dimension - ## Flatten obs2 array - dim(obs2) <- c(dim(obs2)[1],prod(dim(obs2)[2:5]),dim(obs2)[6],dim(obs2)[7]) - if (var != "prec"){ - if (is.null(indx)){ - stop("TimeDisaggregate MUST first be run on precipitation to select appropriate analog index indx") - } else { - ## RUN Disagg function on data using analog day indx - for (imb in 1:dim(exp_corr)[2]){ - for (sd in 1:dim(exp_corr)[3]){ - for (ltime in 1:dim(exp_corr)[4]){ - for (ilat in 1:dim(exp_corr)[5]){ - for (ilon in 1:dim(exp_corr)[6]){ - exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) - } - } - } - } - } - } - } else { - if (is.null(indx)){ - ## Select appropriate analog among 5 possible -> indx array - indx <- array(dim=dim_exp) - ##for (imb in 1:dim(exp_corr)[2]){ - ## indx[1,imb,,,,] <- mapply(SelPrecAnalog,exp_corr[1,imb,,,,],obs2,analog[1,imb,,,],thres) - ##} - for (imb in 1:dim(exp_corr)[2]){ - for (sd in 1:dim(exp_corr)[3]){ - for (ltime in 1:dim(exp_corr)[4]){ - for (ilat in 1:dim(exp_corr)[5]){ - for (ilon in 1:dim(exp_corr)[6]){ - indx[1,imb,sd,ltime,ilat,ilon] <- SelPrecAnalog(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,,ilat,ilon],analog[1,imb,sd,ltime,],thres) - } - } - } - } - } - ## RUN Disagg function on data - for (imb in 1:dim(exp_corr)[2]){ - for (sd in 1:dim(exp_corr)[3]){ - for (ltime in 1:dim(exp_corr)[4]){ - for (ilat in 1:dim(exp_corr)[5]){ - for (ilon in 1:dim(exp_corr)[6]){ - exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) - } - } - } - } - } - } else { - warning("Using previously defined indx array for TimeDisaggregate: precipitation disaggregation is run again") - for (imb in 1:dim(exp_corr)[2]){ - for (sd in 1:dim(exp_corr)[3]){ - for (ltime in 1:dim(exp_corr)[4]){ - for (ilat in 1:dim(exp_corr)[5]){ - for (ilon in 1:dim(exp_corr)[6]){ - exp_dt[,1,imb,sd,ltime,ilat,ilon] <- Disagg(exp_corr[1,imb,sd,ltime,ilat,ilon],obs2[,indx[1,imb,sd,ltime,ilat,ilon],ilat,ilon],var) - } - } - } - } - } - } - } - ##Rearrange dimension order in exp_dt before output - exp_dt <- aperm(exp_dt,c(2,3,4,5,1,6,7)) # Put subdaily dimension back in place - return(list(disagg=exp_dt, indx=indx)) -} - -- GitLab From 46244d94eaf29a2270c69ac16fa1a8416a4ae18e Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Fri, 20 Nov 2020 12:30:21 +0100 Subject: [PATCH 17/24] fixing automatic test for CST_MultiMetric.R --- tests/testthat/test-CST_MultiMetric.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/testthat/test-CST_MultiMetric.R b/tests/testthat/test-CST_MultiMetric.R index b08b68ba..1a41a83a 100644 --- a/tests/testthat/test-CST_MultiMetric.R +++ b/tests/testthat/test-CST_MultiMetric.R @@ -11,8 +11,7 @@ test_that("basic use case", { attr(exp, 'class') <- 's2dv_cube' attr(obs, 'class') <- 's2dv_cube' - result <- list(data = array(rep(c(rep(1, 9), -rep(0.89999999999999991118215802998747676610950, 3)), 48), + result <- list(data = array(rep(c(rep(1, 9), rep(0, 3)), 48), dim = c(dataset = 3, dataset = 1, statistics = 4, lat = 6, lon = 8)), lat = lat, lon = lon) -- GitLab From c6b1cc4a763bbaf680e7cf089e99d164a3b57a33 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Mon, 30 Nov 2020 14:09:01 +0100 Subject: [PATCH 18/24] Update documentation of Adamont functions and minor edits to AdamontNearestNeighbors --- DESCRIPTION | 2 +- NAMESPACE | 3 ++ R/AdamontNearestNeighbors.R | 8 +++-- R/CST_AdamontAnalog.R | 33 ++++++++++++++++++- R/CST_AdamontQQCorr.R | 29 ++++++++++++++++- man/AdamontAnalog.Rd | 63 +++++++++++++++++++++++++++++++++++++ man/AdamontQQCorr.Rd | 52 ++++++++++++++++++++++++++++++ man/CST_AdamontAnalog.Rd | 63 +++++++++++++++++++++++++++++++++++++ man/CST_AdamontQQCorr.Rd | 50 +++++++++++++++++++++++++++++ man/NearestNeighbors.Rd | 30 ++++++++++++++++++ 10 files changed, 328 insertions(+), 5 deletions(-) create mode 100644 man/AdamontAnalog.Rd create mode 100644 man/AdamontQQCorr.Rd create mode 100644 man/CST_AdamontAnalog.Rd create mode 100644 man/CST_AdamontQQCorr.Rd create mode 100644 man/NearestNeighbors.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 07433f27..4c4935c9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -78,4 +78,4 @@ VignetteBuilder: knitr License: Apache License 2.0 Encoding: UTF-8 LazyData: true -RoxygenNote: 7.0.2 +RoxygenNote: 7.1.1 diff --git a/NAMESPACE b/NAMESPACE index 9f9a0e34..72a4d9a4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -42,6 +42,9 @@ export(SplitDim) export(WeatherRegime) export(as.s2dv_cube) export(s2dv_cube) +import("CSTools,") +import("multiApply,") +import(ClimProjDiags) import(abind) import(ggplot2) import(multiApply) diff --git a/R/AdamontNearestNeighbors.R b/R/AdamontNearestNeighbors.R index 1873dbb8..996a4dfc 100644 --- a/R/AdamontNearestNeighbors.R +++ b/R/AdamontNearestNeighbors.R @@ -8,7 +8,11 @@ #'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment longitudes in \code{$lon} and latitudes in \code{$lat} #'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the reference data on a different grid, with longitudes in \code{$lon} and latitudes in \code{$lat}. #' -#'@return NN +#'@return NN a list, containing the following: +#' min_lon: array of dimensions \code{obs$lon} giving the longitude of closest gridpoint in exp +#' min_lat: array of dimensions \code{obs$lat} giving the latitude of closest gridpoint in exp +#' imin_lon: array of dimensions \code{obs$lon} giving the longitude index of closest gridpoint in exp +#' imin_lat: array of dimensions \code{obs$lat} giving the latitude index of closest gridpoint in exp #' #'@import s2dverification #'@import ncdf4 @@ -90,7 +94,7 @@ NearestNeighbors <- function (exp, obs, method='ADA') { stop("AdamontNearestNeighbors is meant for non-regular lat/lon grids; use e.g. CST_Load to interpolate exp onto obs grid") } - NN=list(dist_min=dist_min, min_lon=min_lon, min_lat=min_lat, imin_lon=imin_lon, imin_lat=imin_lat) + NN=list(min_lon=min_lon, min_lat=min_lat, imin_lon=imin_lon, imin_lat=imin_lat) return(NN) } diff --git a/R/CST_AdamontAnalog.R b/R/CST_AdamontAnalog.R index 47f2d381..34d21a42 100644 --- a/R/CST_AdamontAnalog.R +++ b/R/CST_AdamontAnalog.R @@ -1,4 +1,4 @@ -#'AdamontAnalog finds analogous data in the reference dataset to experiment data based on +#'CST_AdamontAnalog finds analogous data in the reference dataset to experiment data based on #'weather types #' #'@description This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data. @@ -54,6 +54,37 @@ CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, method = 'patt return(analog_vals) } +#'AdamontAnalog finds analogous data in the reference dataset to experiment data based on +#'weather types +#' +#'@description This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data. +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#' +#'@param exp experiment data, such as \code{$data} array of an object of class \code{s2dv_cube}, can be output from quantile correction using CST_AdamontQQCorr.R +#'@param wt_exp corresponding weather types (same dimensions as \code{exp} but lat/lon) +#'@param obs reference data, can also be \code{$data} array of class \code{s2dv_cube}. Note that lat/lon dimensions need to be the same as \code{exp} +#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but lat/lon) +#'@param nanalogs integer defining the number of analog values to return (default: 5) +#'@param method a character string indicating the method used for analog definition +#' Coded are 'pattcorr': pattern correlation +#' 'rain1' (for precip patterns): rain occurrence consistency +#' 'rain01' (for precip patterns): rain occurrence/non occurrence consistency +#'@param thres real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1 +#'@param search_obsdims list of dimensions in \code{obs} along which analogs are searched for +#'@param londim name of longitude dimension +#'@param latdim name of latitude dimension +#'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data +#'@import CSTools, multiApply, ClimProjDiags +#'@examples +#'\dontrun{ +#'wt_exp <- sample(1:3, 15*6*3, replace=T) +#'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) +#'wt_obs <- sample(1:3, 6*3, replace=T) +#'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) +# analog_vals <- AdamontAnalog(exp=lonlat_data$exp$data, obs=lonlat_data$obs$data, wt_exp=wt_exp, wt_obs=wt_obs, nanalogs=2) +#'} + AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs=5, method = 'pattcorr', thres = NULL,search_obsdims = c('member', 'sdate', 'ftime'), londim = 'lon', latdim = 'lat') { # exp: lat, lon, sdate, ftime, member # obs: lat, lon, dims for searching 'sdate' 'ftime'... diff --git a/R/CST_AdamontQQCorr.R b/R/CST_AdamontQQCorr.R index 1916d9a5..2e3017ba 100644 --- a/R/CST_AdamontQQCorr.R +++ b/R/CST_AdamontQQCorr.R @@ -1,4 +1,4 @@ -#'ADAMONT QQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types +#'CST_AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types #' #'@description This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference \code{obs}, typically provided by reanalysis data. #'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version @@ -44,6 +44,33 @@ CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member','s return(exp_corr) } + +#'AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types +#' +#'@description This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference \code{obs}, typically provided by reanalysis data. +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#' +#'@param exp array with named dimensions (such as \code{$data} array of experiment data from an object of class \code{s2dv_cube}) +#'@param wt_exp corresponding weather types (same dimensions as \code{exp} but lat/lon) +#'@param obs array with named dimensions with reference data (can also be \code{$data} array of class \code{s2dv_cube}). lat/lon dimensions can differ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and AdamontNearestNeighbors \code{NN} output should be provided +#'@param corrdims list of dimensions in \code{exp} for which quantile mapping correction is applied +#'@param londim character name of longitude dimension in \code{exp} and \code{obs} +#'@param latdim character name of latitude dimension in \code{exp} and \code{obs} +#'@param regrid boolean indicating whether Nearest Neighbors regridding is needed +#'@param NN list (output from AdamontNearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o) +#' +#'@import qmap +#'@import multiApply +#'@import abind +#'@examples +#'\dontrun{ +#'wt_exp <- sample(1:3, 15*6*3, replace=T) +#'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) +#'wt_obs <- sample(1:3, 6*3, replace=T) +#'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) +#'exp_corr <- AdamontQQCorr(exp=lonlat_data$exp$data, wt_exp=wt_exp, obs=lonlat_data$obs$data, wt_obs=wt_obs, corrdims = c('dataset','member','sdate','ftime')) +#'} AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member', 'sdate', 'ftime'), londim='lon', latdim='lat', regrid=FALSE, NN=NULL) { # The regridding part should only be done if lat/lon dimensions of obs and diff --git a/man/AdamontAnalog.Rd b/man/AdamontAnalog.Rd new file mode 100644 index 00000000..576f9ce0 --- /dev/null +++ b/man/AdamontAnalog.Rd @@ -0,0 +1,63 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_AdamontAnalog.R +\name{AdamontAnalog} +\alias{AdamontAnalog} +\title{AdamontAnalog finds analogous data in the reference dataset to experiment data based on +weather types} +\usage{ +AdamontAnalog( + exp, + obs, + wt_exp, + wt_obs, + nanalogs = 5, + method = "pattcorr", + thres = NULL, + search_obsdims = c("member", "sdate", "ftime"), + londim = "lon", + latdim = "lat" +) +} +\arguments{ +\item{exp}{experiment data, such as \code{$data} array of an object of class \code{s2dv_cube}, can be output from quantile correction using CST_AdamontQQCorr.R} + +\item{obs}{reference data, can also be \code{$data} array of class \code{s2dv_cube}. Note that lat/lon dimensions need to be the same as \code{exp}} + +\item{wt_exp}{corresponding weather types (same dimensions as \code{exp} but lat/lon)} + +\item{wt_obs}{corresponding weather types (same dimensions as \code{obs} but lat/lon)} + +\item{nanalogs}{integer defining the number of analog values to return (default: 5)} + +\item{method}{a character string indicating the method used for analog definition +Coded are 'pattcorr': pattern correlation + 'rain1' (for precip patterns): rain occurrence consistency + 'rain01' (for precip patterns): rain occurrence/non occurrence consistency} + +\item{thres}{real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1} + +\item{search_obsdims}{list of dimensions in \code{obs} along which analogs are searched for} + +\item{londim}{name of longitude dimension} + +\item{latdim}{name of latitude dimension} +} +\value{ +analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data +} +\description{ +This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data. +} +\examples{ +\dontrun{ +wt_exp <- sample(1:3, 15*6*3, replace=T) +dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) +wt_obs <- sample(1:3, 6*3, replace=T) +dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) +} +} +\author{ +Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version + +Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +} diff --git a/man/AdamontQQCorr.Rd b/man/AdamontQQCorr.Rd new file mode 100644 index 00000000..d37a7086 --- /dev/null +++ b/man/AdamontQQCorr.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_AdamontQQCorr.R +\name{AdamontQQCorr} +\alias{AdamontQQCorr} +\title{AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types} +\usage{ +AdamontQQCorr( + exp, + wt_exp, + obs, + wt_obs, + corrdims = c("member", "sdate", "ftime"), + londim = "lon", + latdim = "lat", + regrid = FALSE, + NN = NULL +) +} +\arguments{ +\item{exp}{array with named dimensions (such as \code{$data} array of experiment data from an object of class \code{s2dv_cube})} + +\item{wt_exp}{corresponding weather types (same dimensions as \code{exp} but lat/lon)} + +\item{obs}{array with named dimensions with reference data (can also be \code{$data} array of class \code{s2dv_cube}). lat/lon dimensions can differ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and AdamontNearestNeighbors \code{NN} output should be provided} + +\item{corrdims}{list of dimensions in \code{exp} for which quantile mapping correction is applied} + +\item{londim}{character name of longitude dimension in \code{exp} and \code{obs}} + +\item{latdim}{character name of latitude dimension in \code{exp} and \code{obs}} + +\item{regrid}{boolean indicating whether Nearest Neighbors regridding is needed} + +\item{NN}{list (output from AdamontNearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o)} +} +\description{ +This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference \code{obs}, typically provided by reanalysis data. +} +\examples{ +\dontrun{ +wt_exp <- sample(1:3, 15*6*3, replace=T) +dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) +wt_obs <- sample(1:3, 6*3, replace=T) +dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) +exp_corr <- AdamontQQCorr(exp=lonlat_data$exp$data, wt_exp=wt_exp, obs=lonlat_data$obs$data, wt_obs=wt_obs, corrdims = c('dataset','member','sdate','ftime')) +} +} +\author{ +Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version + +Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +} diff --git a/man/CST_AdamontAnalog.Rd b/man/CST_AdamontAnalog.Rd new file mode 100644 index 00000000..579e8304 --- /dev/null +++ b/man/CST_AdamontAnalog.Rd @@ -0,0 +1,63 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_AdamontAnalog.R +\name{CST_AdamontAnalog} +\alias{CST_AdamontAnalog} +\title{CST_AdamontAnalog finds analogous data in the reference dataset to experiment data based on +weather types} +\usage{ +CST_AdamontAnalog( + exp, + obs, + wt_exp, + wt_obs, + nanalogs, + method = "pattcorr", + thres = NULL, + search_obsdims = c("member", "sdate", "ftime"), + londim = "lon", + latdim = "lat" +) +} +\arguments{ +\item{exp}{experiment data an object of class \code{s2dv_cube}, can be output from quantile correction using CST_AdamontQQCorr.R} + +\item{obs}{reference data, also of class \code{s2dv_cube}. Note that lat/lon dimensions need to be the same as \code{exp}} + +\item{wt_exp}{corresponding weather types (same dimensions as \code{exp$data} but lat/lon)} + +\item{wt_obs}{corresponding weather types (same dimensions as \code{obs$data} but lat/lon)} + +\item{nanalogs}{integer defining the number of analog values to return (default: 5)} + +\item{method}{a character string indicating the method used for analog definition +Coded are 'pattcorr': pattern correlation + 'rain1' (for precip patterns): rain occurrence consistency + 'rain01' (for precip patterns): rain occurrence/non occurrence consistency} + +\item{thres}{real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1} + +\item{search_obsdims}{list of dimensions in \code{obs} along which analogs are searched for} + +\item{londim}{name of longitude dimension} + +\item{latdim}{name of latitude dimension} +} +\value{ +analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data +} +\description{ +This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data. +} +\examples{ +\dontrun{ +wt_exp <- sample(1:3, 15*6*3, replace=T) +dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) +wt_obs <- sample(1:3, 6*3, replace=T) +dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) +} +} +\author{ +Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version + +Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +} diff --git a/man/CST_AdamontQQCorr.Rd b/man/CST_AdamontQQCorr.Rd new file mode 100644 index 00000000..dcd31017 --- /dev/null +++ b/man/CST_AdamontQQCorr.Rd @@ -0,0 +1,50 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_AdamontQQCorr.R +\name{CST_AdamontQQCorr} +\alias{CST_AdamontQQCorr} +\title{CST_AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types} +\usage{ +CST_AdamontQQCorr( + exp, + wt_exp, + obs, + wt_obs, + corrdims = c("member", "sdate", "ftime"), + londim = "lon", + latdim = "lat" +) +} +\arguments{ +\item{exp}{experiment data an object of class \code{s2dv_cube}} + +\item{wt_exp}{corresponding weather types (same dimensions as \code{exp$data} but lat/lon)} + +\item{obs}{reference data, also of class \code{s2dv_cube}. lat/lon dimensions can differ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and AdamontNearestNeighbors \code{NN} output should be provided} + +\item{corrdims}{list of dimensions in \code{exp} for which quantile mapping correction is applied} + +\item{londim}{character name of longitude dimension in \code{exp} and \code{obs}} + +\item{latdim}{character name of latitude dimension in \code{exp} and \code{obs}} + +\item{regrid}{boolean indicating whether Nearest Neighbors regridding is needed} + +\item{NN}{list (output from Adamont_NearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o)} +} +\description{ +This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference \code{obs}, typically provided by reanalysis data. +} +\examples{ +\dontrun{ +wt_exp <- sample(1:3, 15*6*3, replace=T) +dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) +wt_obs <- sample(1:3, 6*3, replace=T) +dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) +exp_corr <- CST_AdamontQQCorr(exp=lonlat_data$exp, wt_exp=wt_exp, obs=lonlat_data$obs, wt_obs=wt_obs, corrdims = c('dataset','member','sdate','ftime')) +} +} +\author{ +Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version + +Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +} diff --git a/man/NearestNeighbors.Rd b/man/NearestNeighbors.Rd new file mode 100644 index 00000000..11e6b2ff --- /dev/null +++ b/man/NearestNeighbors.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AdamontNearestNeighbors.R +\name{NearestNeighbors} +\alias{NearestNeighbors} +\title{ADAMONT Nearest Neighbors computes the distance between reference data grid centroid and SF data grid} +\usage{ +NearestNeighbors(exp, obs, method = "ADA") +} +\arguments{ +\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment longitudes in \code{$lon} and latitudes in \code{$lat}} + +\item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the reference data on a different grid, with longitudes in \code{$lon} and latitudes in \code{$lat}.} + +\item{method}{a string among three options ('ADA': standard ADAMONT distance, 'simple': lon/lat straight euclidian distance, 'radius': distance on the sphere)} +} +\value{ +NN a list, containing the following: + min_lon: array of dimensions \code{obs$lon} giving the longitude of closest gridpoint in exp + min_lat: array of dimensions \code{obs$lat} giving the latitude of closest gridpoint in exp + imin_lon: array of dimensions \code{obs$lon} giving the longitude index of closest gridpoint in exp + imin_lat: array of dimensions \code{obs$lat} giving the latitude index of closest gridpoint in exp +} +\description{ +This function computes the nearest neighbor for each reference data (lon, lat) point in the experiment dataset by computing the distance between the reference dataset grid and the experiment data. This is the first step in the ADAMONT method adapted from Verfaillie et al. (2018). +} +\author{ +Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version + +Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +} -- GitLab From 107db33896316bc1988388a8763174e2f136d740 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Mon, 30 Nov 2020 14:21:06 +0100 Subject: [PATCH 19/24] Fixing errors in pipeline --- NAMESPACE | 2 -- R/CST_AdamontAnalog.R | 6 ++++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 72a4d9a4..d8ba17fa 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -42,8 +42,6 @@ export(SplitDim) export(WeatherRegime) export(as.s2dv_cube) export(s2dv_cube) -import("CSTools,") -import("multiApply,") import(ClimProjDiags) import(abind) import(ggplot2) diff --git a/R/CST_AdamontAnalog.R b/R/CST_AdamontAnalog.R index 34d21a42..a2036878 100644 --- a/R/CST_AdamontAnalog.R +++ b/R/CST_AdamontAnalog.R @@ -19,7 +19,8 @@ #'@param londim name of longitude dimension #'@param latdim name of latitude dimension #'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data -#'@import CSTools, multiApply, ClimProjDiags +#'@import multiApply +#'@import ClimProjDiags #'@examples #'\dontrun{ #'wt_exp <- sample(1:3, 15*6*3, replace=T) @@ -75,7 +76,8 @@ CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, method = 'patt #'@param londim name of longitude dimension #'@param latdim name of latitude dimension #'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data -#'@import CSTools, multiApply, ClimProjDiags +#'@import multiApply +#'@import ClimProjDiags #'@examples #'\dontrun{ #'wt_exp <- sample(1:3, 15*6*3, replace=T) -- GitLab From 62f0994c29784a4380316f682bb15983918f2945 Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Mon, 15 Feb 2021 11:57:38 +0100 Subject: [PATCH 20/24] Minor edits to make NearestNeighbors invisible to users --- R/AdamontNearestNeighbors.R | 2 +- R/CST_AdamontQQCorr.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/AdamontNearestNeighbors.R b/R/AdamontNearestNeighbors.R index 996a4dfc..4c6dbe07 100644 --- a/R/AdamontNearestNeighbors.R +++ b/R/AdamontNearestNeighbors.R @@ -17,7 +17,7 @@ #'@import s2dverification #'@import ncdf4 -NearestNeighbors <- function (exp, obs, method='ADA') { +.NearestNeighbors <- function (exp, obs, method='ADA') { if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ", diff --git a/R/CST_AdamontQQCorr.R b/R/CST_AdamontQQCorr.R index 2e3017ba..0c810e4c 100644 --- a/R/CST_AdamontQQCorr.R +++ b/R/CST_AdamontQQCorr.R @@ -35,7 +35,7 @@ CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member','s stop("'londim' or 'latdim' input doesn't match with 'obs' dimension names") } if ((length(dim(exp$lon))==2) || (length(dim(obs$lon))==2)){ - myNN <- NearestNeighbors(exp=exp, obs=obs, method='ADA') + myNN <- .NearestNeighbors(exp=exp, obs=obs, method='ADA') exp_corr <- AdamontQQCorr(exp=exp$data, wt_exp=wt_exp, obs=obs$data, wt_obs=wt_obs, corrdims=corrdims, londim=londim, latdim=latdim, regrid=TRUE, NN=myNN) } else { # If not (standard case) -- GitLab From 43f273cf2d57b2098de28e91a269783d44b7986e Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Thu, 18 Feb 2021 16:17:16 +0100 Subject: [PATCH 21/24] Revised documentation and new sanity checks in functions CST_AdamontAnalog, AdamontAnalog, CST_AdamontQQCorr and AdamontQQCorr following expert review --- DESCRIPTION | 5 +- R/CST_AdamontAnalog.R | 182 ++++++++++++++++++++++--------- R/CST_AdamontQQCorr.R | 228 +++++++++++++++++++++++++++++---------- man/AdamontAnalog.Rd | 45 +++++--- man/AdamontQQCorr.Rd | 44 ++++++-- man/CST_AdamontAnalog.Rd | 42 +++++--- man/CST_AdamontQQCorr.Rd | 45 +++++--- 7 files changed, 434 insertions(+), 157 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 4c4935c9..12fe8f75 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,16 +13,18 @@ Authors@R: c( person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "aut"), person("Veronica", "Torralba", , "veronica.torralba@bsc.es", role = "aut"), person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut"), - person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "ctb"), + person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "aut"), person("Filippo", "Cali Quaglia", , "filippo.caliquaglia@gmail.com", role = "ctb"), person("Chihchung", "Chou", , "chihchung.chou@bsc.es", role = "ctb"), person("Nicola", "Cortesi", , "nicola.cortesi@bsc.es", role = "ctb"), person("Susanna", "Corti", , "s.corti@isac.cnr.it", role = "ctb"), person("Paolo", "Davini", , "p.davini@isac.cnr.it", role = "ctb"), + person("Gildas", "Dayon", , "gildas.dayon@meteo.fr", role = "ctb"), person("Marta", "Dominguez", , "mdomingueza@aemet.es", role = "ctb"), person("Federico", "Fabiano", , "f.fabiano@isac.cnr.it", role = "ctb"), person("Ignazio", "Giuntoli", , "i.giuntoli@isac.cnr.it", role = "ctb"), person("Raul", "Marcos", , "raul.marcos@bsc.es", role = "ctb"), + person("Paola", "Marson", , "paola.marson@meteo.fr", role = "ctb"), person("Niti", "Mishra", , "niti.mishra@bsc.es", role = "ctb"), person("Jesus", "Peña", , "jesus.pena@bsc.es", role = "ctb"), person("Francesc", "Roura-Adserias", , "francesc.roura@bsc.es", role = "ctb"), @@ -43,6 +45,7 @@ Description: Exploits dynamical seasonal forecasts in order to provide Terzago et al. (2018) . Torralba et al. (2017) . D'Onofrio et al. (2014) . + Verfaillie et al. (2017) . Van Schaeybroeck et al. (2019) . Yiou et al. (2013) . Depends: diff --git a/R/CST_AdamontAnalog.R b/R/CST_AdamontAnalog.R index a2036878..83f4795b 100644 --- a/R/CST_AdamontAnalog.R +++ b/R/CST_AdamontAnalog.R @@ -1,24 +1,36 @@ -#'CST_AdamontAnalog finds analogous data in the reference dataset to experiment data based on -#'weather types +#'CST_AdamontAnalog finds analogous data in the reference dataset to experiment +#'data based on weather types #' -#'@description This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data. +#'@description This function searches for analogs in a reference dataset for +#'experiment data, based on corresponding weather types. The experiment data is +#'typically a hindcast, observations are typically provided by reanalysis data. #'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version #'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation #' -#'@param exp experiment data an object of class \code{s2dv_cube}, can be output from quantile correction using CST_AdamontQQCorr.R -#'@param wt_exp corresponding weather types (same dimensions as \code{exp$data} but lat/lon) -#'@param obs reference data, also of class \code{s2dv_cube}. Note that lat/lon dimensions need to be the same as \code{exp} -#'@param wt_obs corresponding weather types (same dimensions as \code{obs$data} but lat/lon) -#'@param nanalogs integer defining the number of analog values to return (default: 5) -#'@param method a character string indicating the method used for analog definition +#'@param exp experiment data an object of class \code{s2dv_cube}, can be output +#'from quantile correction using CST_AdamontQQCorr.R +#'@param wt_exp corresponding weather types (same dimensions as \code{exp$data} +#'but lat/lon) +#'@param obs reference data, also of class \code{s2dv_cube}. Note that lat/lon +#'dimensions need to be the same as \code{exp} +#'@param wt_obs corresponding weather types (same dimensions as \code{obs$data} +#'but lat/lon) +#'@param nanalogs integer defining the number of analog values to return +#'(default: 5) +#'@param method a character string indicating the method used for analog +#'definition #' Coded are 'pattcorr': pattern correlation #' 'rain1' (for precip patterns): rain occurrence consistency -#' 'rain01' (for precip patterns): rain occurrence/non occurrence consistency -#'@param thres real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1 -#'@param search_obsdims list of dimensions in \code{obs} along which analogs are searched for +#' 'rain01' (for precip patterns): rain occurrence/non +#' occurrence consistency +#'@param thres real number indicating the threshold to define rain +#'occurrence/non occurrence in rain(0)1 +#'@param search_obsdims list of dimensions in \code{obs} along which analogs are +#'searched for #'@param londim name of longitude dimension #'@param latdim name of latitude dimension -#'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data +#'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs +#'analog values for each value of \code{exp} input data #'@import multiApply #'@import ClimProjDiags #'@examples @@ -29,8 +41,13 @@ #'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) # analog_vals <- CST_AdamontAnalog(exp=lonlat_data$exp, obs=lonlat_data$obs, wt_exp=wt_exp, wt_obs=wt_obs, nanalogs=2) #'} -CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, method = 'pattcorr', thres = NULL, search_obsdims = c('member', 'sdate', 'ftime'), londim = 'lon', latdim = 'lat') { +CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, + method = 'pattcorr', thres = NULL, + search_obsdims = c('member', 'sdate', 'ftime'), + londim = 'lon', latdim = 'lat') { + dimnames <- names(dim(obs$data)) + dimnamesexp <- names(dim(exp$data)) if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ", "as output by CSTools::CST_Load.") @@ -41,11 +58,30 @@ CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, method = 'patt if (is.null(nanalogs)){ nanalogs <- 5 } - dimnames <- names(dim(obs$data)) - if (is.null(which(dimnames == latdim)) || is.null(which(dimnames == londim))) { - stop("'londim' or 'latdim' inputs should correspond to 'obs' dimension names") + if (is.null(which(dimnames == latdim)) || is.null(which(dimnames == londim))){ + stop("'londim' or 'latdim' inputs should correspond to 'obs$data' ", + "dimension names") + } + if (!all(search_obsdims %in% dimnames)){ + stop("Names in parameter 'search_obsdims' should match 'obs$data' ", + "dimension names.") + } + if (!all(dim(wt_exp) %in% dim(exp$data))){ + stop("Dimensions for 'wt_exp' should match 'exp$data' except lat/lon") + } + if (!all(dim(wt_obs) %in% dim(obs$data))){ + stop("Dimensions for 'wt_obs' should match 'obs$data' except lat/lon") + } + plat_exp <- which(dimnamesexp==latdim) + plon_exp <- which(dimnamesexp==londim) + plat_obs <- which(dimnames==latdim) + plon_obs <- which(dimnames==londim) + if ((dim(obs$data)[plon_obs]!=dim(exp$data)[plon_exp]) || + (dim(obs$data)[plat_obs]!=dim(exp$data)[plat_exp])){ + stop("Element 'data' from parameters 'obs' and 'exp' should have", + "same lon / lat dimensions if working with regular grids.") } - # More sanity checks needed? Same dimensions in obs and exp; wt_obs and wt_exp + # End of sanity checks; call AdamontAnalog function analog_vals <- AdamontAnalog(exp = exp$data, obs = obs$data, wt_exp = wt_exp, wt_obs=wt_obs, nanalogs = nanalogs, method = method, thres = thres, @@ -55,27 +91,41 @@ CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, method = 'patt return(analog_vals) } -#'AdamontAnalog finds analogous data in the reference dataset to experiment data based on -#'weather types +#'AdamontAnalog finds analogous data in the reference dataset to experiment data +#'based on weather types #' -#'@description This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data. +#'@description This function searches for analogs in a reference dataset for +#'experiment data, based on corresponding weather types. The experiment data is +#'typically a hindcast, observations are typically provided by reanalysis data. #'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version #'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation #' -#'@param exp experiment data, such as \code{$data} array of an object of class \code{s2dv_cube}, can be output from quantile correction using CST_AdamontQQCorr.R -#'@param wt_exp corresponding weather types (same dimensions as \code{exp} but lat/lon) -#'@param obs reference data, can also be \code{$data} array of class \code{s2dv_cube}. Note that lat/lon dimensions need to be the same as \code{exp} -#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but lat/lon) -#'@param nanalogs integer defining the number of analog values to return (default: 5) -#'@param method a character string indicating the method used for analog definition +#'@param exp experiment data, such as \code{$data} array of an object of class +#'\code{s2dv_cube}, can be output from quantile correction using +#'\code{CST_AdamontQQCorr.R} +#'@param wt_exp corresponding weather types (same dimensions as \code{exp} but +#'lat/lon) +#'@param obs reference data, can also be \code{$data} array of class +#'\code{s2dv_cube}. Note that lat/lon dimensions need to be the same as +#'\code{exp} +#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but * +#'lat/lon) +#'@param nanalogs integer defining the number of analog values to return +#'(default: 5) +#'@param method a character string indicating the method used for analog +#'definition #' Coded are 'pattcorr': pattern correlation #' 'rain1' (for precip patterns): rain occurrence consistency -#' 'rain01' (for precip patterns): rain occurrence/non occurrence consistency -#'@param thres real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1 -#'@param search_obsdims list of dimensions in \code{obs} along which analogs are searched for +#' 'rain01' (for precip patterns): rain occurrence/non +#' occurrence consistency +#'@param thres real number indicating the threshold to define rain occurrence +#'/non occurrence in rain(0)1 +#'@param search_obsdims list of dimensions in \code{obs} along which analogs are +#'searched for #'@param londim name of longitude dimension #'@param latdim name of latitude dimension -#'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data +#'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs +#'analog values for each value of \code{exp} input data #'@import multiApply #'@import ClimProjDiags #'@examples @@ -84,39 +134,70 @@ CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, method = 'patt #'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) #'wt_obs <- sample(1:3, 6*3, replace=T) #'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) -# analog_vals <- AdamontAnalog(exp=lonlat_data$exp$data, obs=lonlat_data$obs$data, wt_exp=wt_exp, wt_obs=wt_obs, nanalogs=2) +# analog_vals <- AdamontAnalog(exp=lonlat_data$exp$data, +#' obs=lonlat_data$obs$data, wt_exp=wt_exp, wt_obs=wt_obs, nanalogs=2) #'} -AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs=5, method = 'pattcorr', thres = NULL,search_obsdims = c('member', 'sdate', 'ftime'), londim = 'lon', latdim = 'lat') { +AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs=5, + method = 'pattcorr', thres = NULL, + search_obsdims = c('member', 'sdate', 'ftime'), + londim = 'lon', latdim = 'lat') { # exp: lat, lon, sdate, ftime, member # obs: lat, lon, dims for searching 'sdate' 'ftime'... # wt_exp: sdate, ftime, member # wt_obs: the dims for searching + dimnames <- names(dim(obs)) + dimnamesexp <- names(dim(exp)) if (method %in% c('rain1','rain01') & is.null(thres)){ stop("Threshold 'thres' must be defined with methods 'rain1' and 'rain01'") } if (method == 'pattcorr' & !is.null(thres)){ warning("Parameter 'thres' is not used with method 'pattcorr'.") } + if (is.null(which(dimnamesexp == latdim)) || is.null(which(dimnamesexp==londim))){ + stop("'londim' or 'latdim' input doesn't match with 'exp' dimension names") + } # Position of lat/lon dimensions in exp data - poslatexp <- which(names(dim(exp)) == latdim) - poslonexp <- which(names(dim(exp)) == londim) - # Reshaping obs: + poslatexp <- which(dimnamesexp == latdim) + poslonexp <- which(dimnamesexp == londim) + poslatobs <- which(dimnames == latdim) + poslonobs <- which(dimnames == londim) + if (!all(search_obsdims %in% dimnames)){ + stop("Names in parameter 'search_obsdims' should match 'obs' ", + "dimension names.") + } + if (!all(dim(wt_exp) %in% dim(exp))){ + stop("Dimensions for 'wt_exp' should match 'exp' except lat/lon") + } + if (!all(dim(wt_obs) %in% dim(obs))){ + stop("Dimensions for 'wt_obs' should match 'obs' except lat/lon") + } + if ((dim(obs)[poslonobs]!=dim(exp)[poslonexp]) || + (dim(obs)[poslatobs]!=dim(exp)[poslatexp])){ + stop("Parameters 'obs' and 'exp' should have same lon / lat dimensions.") + } + + ## Reshaping obs: ## The dimensions where to search in a single dim if (length(search_obsdims) > 1) { for (i in 1:(length(search_obsdims) - 1)) { - obs <- MergeDims(obs, search_obsdims[i:(i + 1)], rename = search_obsdims[i + 1]) - wt_obs <- MergeDims(wt_obs, search_obsdims[i:(i + 1)], rename = search_obsdims[i + 1]) + obs <- MergeDims(obs, search_obsdims[i:(i + 1)], + rename = search_obsdims[i + 1]) + wt_obs <- MergeDims(wt_obs, search_obsdims[i:(i + 1)], + rename = search_obsdims[i + 1]) } } names(dim(obs))[which(names(dim(obs)) == search_obsdims[length(search_obsdims)])] <- 'time' names(dim(wt_obs))[which(names(dim(wt_obs)) == search_obsdims[length(search_obsdims)])] <- 'time' # Split 'time' dim in weather types - obs <- SplitDim(obs, split_dim = 'time', indices = as.vector(wt_obs), new_dim_name='type') + obs <- SplitDim(obs, split_dim = 'time', indices = as.vector(wt_obs), + new_dim_name='type') - analog_vals <- Apply(list(exp, obs, wt_exp), target_dims = list(c(londim, latdim), - c(londim, latdim, 'time', 'type'), NULL), - .analogs, method = method, thres = thres)$output1 + analog_vals <- Apply(list(exp, obs, wt_exp), + target_dims = list(c(londim, latdim), + c(londim, latdim, 'time', 'type'), + NULL), + .analogs, method = method, thres = thres)$output1 # Reshaping output: analog_vals <- Subset(analog_vals,along='type',indices=1,drop='selected') @@ -125,9 +206,11 @@ AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs=5, method = 'pattc postime <- which(names(dim(analog_vals)) == 'time') # Dimension with N analogs pos <- 1:length(dim(analog_vals)) if (poslatexp > poslonexp){ - analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslon,poslat)) + analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)], + postime,poslon,poslat)) } else { - analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],postime,poslat,poslon)) + analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)], + postime,poslat,poslon)) } # Renaming 'time' dim to 'analog' names(dim(analog_vals))[which(names(dim(analog_vals)) == 'time')] <- 'analog' @@ -135,16 +218,17 @@ AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs=5, method = 'pattc } -.analogs <- function(exp, obs, wt_exp, nanalogs = 5, method = 'pattcorr', thres = NULL, - londimexp = 'lon', latdimexp = 'lat', londimobs = 'lon', latdimobs = 'lat') { +.analogs <- function(exp, obs, wt_exp, nanalogs = 5, method = 'pattcorr', + thres = NULL, londimexp = 'lon', latdimexp = 'lat', + londimobs = 'lon', latdimobs = 'lat') { # exp: lon, lat # obs: lon, lat, time, wt # wt_exp: wt single scalar search_analog <- switch(method, 'rain1' = .rain1, 'rain01' = .rain01, - 'pattcorr' = .pattcor, - stop(paste0("Adamont Analog function supports methods ", - "'rain1', 'rain01', 'pattcorr' only"))) + 'pattcorr' = .pattcor, + stop(paste0("Adamont Analog function only supports ", + "methods 'rain1', 'rain01', 'pattcorr'"))) obs <- Subset(obs, along = 'type', indices = wt_exp) accuracy <- Apply(list(exp, obs), target_dims = list(c(londimexp, latdimexp), diff --git a/R/CST_AdamontQQCorr.R b/R/CST_AdamontQQCorr.R index 0c810e4c..c8719a88 100644 --- a/R/CST_AdamontQQCorr.R +++ b/R/CST_AdamontQQCorr.R @@ -1,19 +1,37 @@ -#'CST_AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types +#'CST_AdamontQQCorr computes quantile-quantile correction of seasonal or +#'decadal forecast data using weather types #' -#'@description This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference \code{obs}, typically provided by reanalysis data. -#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version -#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#'@description This function computes a quantile mapping based on weather types +#'for experiment data (typically a hindcast) onto reference \code{obs}, +#'typically provided by reanalysis data. +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} +#'@author Paola Marson, \email{paola.marson@meteo.fr} +#'@author Gildas Dayon, \email{gildas.dayon@meteo.fr} #' #'@param exp experiment data an object of class \code{s2dv_cube} -#'@param wt_exp corresponding weather types (same dimensions as \code{exp$data} but lat/lon) -#'@param obs reference data, also of class \code{s2dv_cube}. lat/lon dimensions can differ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and AdamontNearestNeighbors \code{NN} output should be provided -#'@param corrdims list of dimensions in \code{exp} for which quantile mapping correction is applied -#'@param londim character name of longitude dimension in \code{exp} and \code{obs} -#'@param latdim character name of latitude dimension in \code{exp} and \code{obs} -#'@param regrid boolean indicating whether Nearest Neighbors regridding is needed -#'@param NN list (output from Adamont_NearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o) +#'@param wt_exp corresponding weather types (same dimensions as \code{exp$data} +#' but lat/lon) +#'@param obs reference data, also of class \code{s2dv_cube}. lat/lon dimensions +#' can differ from \code{exp} if non rectilinear latlon grids are used, +#' in which case regrid should be set to TRUE and .NearestNeighbors \code{NN} +#' output should be provided +#'@param corrdims list of dimensions in \code{exp} for which quantile mapping +#' correction is applied +#'@param londim character name of longitude dimension in \code{exp} and +#' \code{obs} +#'@param latdim character name of latitude dimension in \code{exp} and +#' \code{obs} +#'@param regrid (optional) boolean indicating whether .NearestNeighbors +#' regridding is needed +#'@param NN list (optional, if regrid=TRUE) output from .NearestNeighbors +#' mapping (nlat, nlon) onto (nlat_o, nlon_o) +#' +#'@return an object of class \code{s2dv_cube} containing experiment data on the +#' lat/lon grid of \code{obs} input data, corrected by quantile mapping +#' depending on the weather types \code{wt_exp} #' #'@import qmap +#'@import ClimProjDiags #'@import multiApply #'@import abind #'@examples @@ -22,45 +40,104 @@ #'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) #'wt_obs <- sample(1:3, 6*3, replace=T) #'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) -#'exp_corr <- CST_AdamontQQCorr(exp=lonlat_data$exp, wt_exp=wt_exp, obs=lonlat_data$obs, wt_obs=wt_obs, corrdims = c('dataset','member','sdate','ftime')) +#'exp_corr <- CST_AdamontQQCorr(exp=lonlat_data$exp, wt_exp=wt_exp, +#' obs=lonlat_data$obs, wt_obs=wt_obs, +#' corrdims = c('dataset','member','sdate','ftime')) #'} -CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member','sdate','ftime'), londim='lon', latdim='lat') { +CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, + corrdims = c('member','sdate','ftime'), + londim='lon', latdim='lat') { if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')){ stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ", "as output by CSTools::CST_Load.") } dimnames <- names(dim(obs$data)) - if (is.null(which(dimnames == latdim)) || is.null(which(dimnames==londim))){ - stop("'londim' or 'latdim' input doesn't match with 'obs' dimension names") + dimnamesexp <- names(dim(exp$data)) + if (is.null(which(dimnames==latdim)) || is.null(which(dimnames==londim))){ + stop("'londim' or 'latdim' input doesn't match with 'obs$data' + dimension names") + } + if (is.null(which(dimnamesexp==latdim)) || + is.null(which(dimnamesexp==londim))){ + stop("'londim' or 'latdim' input doesn't match with 'exp$data' + dimension names") + } + if (is.null(which(corrdims=='time')) || is.null(which(corrdims=='ftime'))){ + warning("Forecast time should be one of the dimensions for the correction + specified in corrdims input list") + } + if (!all(corrdims %in% dimnamesexp)){ + stop("Names in parameter 'corrdims' should match input dimension names.") + } + if (!all(dim(wt_exp) %in% dim(exp$data))){ + stop("Dimensions for 'wt_exp' should match 'exp$data' except lat/lon") + } + if (!all(dim(wt_obs) %in% dim(obs$data))){ + stop("Dimensions for 'wt_obs' should match 'obs$data' except lat/lon") } if ((length(dim(exp$lon))==2) || (length(dim(obs$lon))==2)){ myNN <- .NearestNeighbors(exp=exp, obs=obs, method='ADA') - exp_corr <- AdamontQQCorr(exp=exp$data, wt_exp=wt_exp, obs=obs$data, wt_obs=wt_obs, corrdims=corrdims, londim=londim, latdim=latdim, regrid=TRUE, NN=myNN) + exp_corr <- AdamontQQCorr(exp=exp$data, wt_exp=wt_exp, obs=obs$data, + wt_obs=wt_obs, corrdims=corrdims, + londim=londim, latdim=latdim, + regrid=TRUE, NN=myNN) } else { - # If not (standard case) - exp_corr <- AdamontQQCorr(exp=exp$data, wt_exp=wt_exp, obs=obs$data, wt_obs=wt_obs, corrdims=corrdims, londim=londim, latdim=latdim, regrid=FALSE) + ## If not (standard case) + ## exp$data lat/lon dimensions should match obs$data + plat_exp <- which(dimnamesexp==latdim) + plon_exp <- which(dimnamesexp==londim) + plat_obs <- which(dimnames==latdim) + plon_obs <- which(dimnames==londim) + if ((dim(obs$data)[plon_obs]!=dim(exp$data)[plon_exp]) || + (dim(obs$data)[plat_obs]!=dim(exp$data)[plat_exp])){ + stop("Element 'data' from parameters 'obs' and 'exp' should have", + "same lon / lat dimensions if working with regular grids.") + } + exp_corr <- AdamontQQCorr(exp=exp$data, wt_exp=wt_exp, obs=obs$data, + wt_obs=wt_obs, corrdims=corrdims, + londim=londim, latdim=latdim, regrid=FALSE) } return(exp_corr) } -#'AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types +#'AdamontQQCorr computes quantile-quantile correction of seasonal or decadal +#'forecast data using weather types #' -#'@description This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference \code{obs}, typically provided by reanalysis data. +#'@description This function computes a quantile mapping based on weather types +#'for experiment data (typically a hindcast) onto reference \code{obs}, +#'typically provided by reanalysis data. #'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version #'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation #' -#'@param exp array with named dimensions (such as \code{$data} array of experiment data from an object of class \code{s2dv_cube}) -#'@param wt_exp corresponding weather types (same dimensions as \code{exp} but lat/lon) -#'@param obs array with named dimensions with reference data (can also be \code{$data} array of class \code{s2dv_cube}). lat/lon dimensions can differ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and AdamontNearestNeighbors \code{NN} output should be provided -#'@param corrdims list of dimensions in \code{exp} for which quantile mapping correction is applied -#'@param londim character name of longitude dimension in \code{exp} and \code{obs} -#'@param latdim character name of latitude dimension in \code{exp} and \code{obs} -#'@param regrid boolean indicating whether Nearest Neighbors regridding is needed -#'@param NN list (output from AdamontNearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o) +#'@param exp array with named dimensions (such as \code{$data} array of +#'experiment data from an object of class \code{s2dv_cube}) +#'@param wt_exp corresponding weather types (same dimensions as \code{exp} but +#'lat/lon) +#'@param obs array with named dimensions with reference data (can also be +#'\code{$data} array of class \code{s2dv_cube}). lat/lon dimensions can differ +#'from \code{exp} if non rectilinear latlon grids are used, in which case +#'regrid should be set to TRUE and .NearestNeighbors \code{NN} output should be +#'provided +#'@param corrdims list of dimensions in \code{exp} for which quantile mapping +#'correction is applied +#'@param londim character name of longitude dimension in \code{exp} and +#'\code{obs} +#'@param latdim character name of latitude dimension in \code{exp} and +#'\code{obs} +#'@param regrid (optional) boolean indicating whether .NearestNeighbors +#'regridding is needed +#'@param NN (optional, if regrid=TRUE) list (output from .NearestNeighbors) +#'maps (nlat, nlon) onto (nlat_o, nlon_o) +#' +#'@return an array (such as \code{$data} array from an object of class +#'\code{s2dv_cube}) with named dimensions, containing experiment data on the +#'lat/lon grid of \code{obs} array, corrected by quantile mapping depending on +#'the weather types \code{wt_exp} #' #'@import qmap +#'@import ClimProjDiags #'@import multiApply #'@import abind #'@examples @@ -69,14 +146,37 @@ CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member','s #'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) #'wt_obs <- sample(1:3, 6*3, replace=T) #'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) -#'exp_corr <- AdamontQQCorr(exp=lonlat_data$exp$data, wt_exp=wt_exp, obs=lonlat_data$obs$data, wt_obs=wt_obs, corrdims = c('dataset','member','sdate','ftime')) +#'exp_corr <- AdamontQQCorr(exp=lonlat_data$exp$data, wt_exp=wt_exp, +#' obs=lonlat_data$obs$data, wt_obs=wt_obs, +#' corrdims = c('dataset','member','sdate','ftime')) #'} -AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member', 'sdate', 'ftime'), londim='lon', latdim='lat', regrid=FALSE, NN=NULL) { +AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, + corrdims = c('member', 'sdate', 'ftime'), + londim='lon', latdim='lat', regrid=FALSE, NN=NULL) { + dimnames <- names(dim(obs)) + dimnamesexp <- names(dim(exp)) + if (is.null(which(dimnames == latdim)) || is.null(which(dimnames==londim))){ + stop("'londim' or 'latdim' input doesn't match with 'obs' dimension names") + } + if (is.null(which(corrdims=='time')) || is.null(which(corrdims=='ftime'))){ + warning("Forecast time should be one of the dimensions for the correction + specified in corrdims input list") + } + if (!all(corrdims %in% dimnamesexp)){ + stop("Names in parameter 'corrdims' should match input dimension names.") + } + if (!all(dim(wt_exp) %in% dim(exp))){ + stop("Dimensions for 'wt_exp' should match 'exp' except lat/lon") + } + if (!all(dim(wt_obs) %in% dim(obs))){ + stop("Dimensions for 'wt_obs' should match 'obs' except lat/lon") + } + if ((regrid == 'TRUE') & is.null(NN)){ + stop("regrid set to TRUE: provide nearest neighbors input NN") + } # The regridding part should only be done if lat/lon dimensions of obs and # exp differ. - # INCLUDE SANITY CHECKS: LATDIM AND LONDIM MUST ALSO MATCH EXP LAT/LON DIM NAMES - # IF REGRID = TRUE AND NN = NULL: error !! if (regrid == 'TRUE'){ obsdims <- names(dim(obs)) poslat <- which(obsdims == latdim) @@ -87,25 +187,37 @@ AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member', 'sdat names(dim(ilat_o))[1] <- latdim ilon_o <- array(c(1:nlon_o)) names(dim(ilon_o))[1] <- londim - # First step if obs data is higher resolution than exp data is to use nearest neighbor - # to compute downscaling of exp data - exp_corr <- Apply(list(exp,ilat_o,ilon_o),target_dims=list(c(latdim,londim),latdim,londim),.getNN,NN=NN)$output1 + ## First step if obs data is higher resolution than exp data is to use + ## nearest neighbor to compute downscaling of exp data + exp_corr <- Apply(list(exp,ilat_o,ilon_o), + target_dims=list(c(latdim,londim),latdim,londim), + .getNN,NN=NN)$output1 - # Reorder exp_corr dimensions to match exp dimensions + ## Reorder exp_corr dimensions to match exp dimensions dexpc <- match(names(dim(exp)), names(dim(exp_corr))) exp_corr <- aperm(exp_corr,dexpc) dimnames(exp_corr) <- dimnames(exp)[dexpc] - # Keep original wt_exp for remapping data + ## Keep original wt_exp for remapping data wt_exp2 <- wt_exp - - # Both exp and obs data are now on the same grid + ## Both exp and obs data are now on the same grid } else { + ## exp lat/lon dimensions should match obs + plat_exp <- which(dimnamesexp==latdim) + plon_exp <- which(dimnamesexp==londim) + plat_obs <- which(dimnames==latdim) + plon_obs <- which(dimnames==londim) + if ((dim(obs)[plon_obs]!=dim(exp)[plon_exp]) || + (dim(obs)[plat_obs]!=dim(exp)[plat_exp])){ + stop("Parameters 'obs' and 'exp' should have same lon / lat", + " dimensions if regrid set to 'FALSE' (regular grid case).") + } exp_corr <- exp - # Keep original wt_exp for remapping data + ## Keep original wt_exp for remapping data wt_exp2 <- wt_exp } - # Use CST_QuantileMapping function for quantile mapping, depending on weather type + ## Use CST_QuantileMapping function for quantile mapping + ## depending on weather type for (i in 1:(length(corrdims) - 1)) { obs <- MergeDims(obs, corrdims[i:(i+1)], rename=corrdims[i+1]) wt_obs <- MergeDims(wt_obs, corrdims[i:(i+1)], rename=corrdims[i+1]) @@ -117,13 +229,15 @@ AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member', 'sdat names(dim(exp_corr))[which(names(dim(exp_corr)) == corrdims[length(corrdims)])] <- 'time' names(dim(wt_exp2))[which(names(dim(wt_exp2)) == corrdims[length(corrdims)])] <- 'time' # Split 'time' dim in weather types - obs <- SplitDim(obs, split_dim='time',indices=as.vector(wt_obs), new_dim_name='type') - exp_corr <- SplitDim(exp_corr, split_dim='time',indices=as.vector(wt_exp2), new_dim_name='type') - # Add NAs to exp_corr if needed to have compatible sample dimensions + obs <- SplitDim(obs, split_dim='time',indices=as.vector(wt_obs), + new_dim_name='type') + exp_corr <- SplitDim(exp_corr, split_dim='time',indices=as.vector(wt_exp2), + new_dim_name='type') + ## Add NAs to exp_corr if needed to have compatible sample dimensions numtobs <- dim(obs)[which(names(dim(obs))=='time')] numtexp <- dim(exp_corr)[which(names(dim(exp_corr))=='time')] if (numtexp%%numtobs > 0){ - # Create extra dimension and include NAs + ## Create extra dimension and include NAs ndimexp <- names(dim(exp_corr)) ndimobs <- names(dim(obs)) postime <- which(ndimexp=='time') @@ -135,23 +249,26 @@ AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member', 'sdat dimobs <- c(dim(obs),1) dim(obs) <- dimobs names(dim(obs)) <- c(ndimobs,'index') - res <- QuantileMapping(exp=exp_corr,obs=obs,sample_dims=c('time','index'),method='RQUANT') + res <- QuantileMapping(exp=exp_corr,obs=obs,sample_dims=c('time','index'), + method='RQUANT') res <- MergeDims(res,c('time','index')) - # Remove the extra NA values added previously + ## Remove the extra NA values added previously res <- Subset(res,along='time',indices=1:numtexp) } else { - # Apply QuantileMapping to exp_corr depending on weather type - res <- QuantileMapping(exp=exp_corr,obs=obs,sample_dims='time',samplemethod='RQUANT') + ## Apply QuantileMapping to exp_corr depending on weather type + res <- QuantileMapping(exp=exp_corr,obs=obs,sample_dims='time', + samplemethod='RQUANT') } rm(exp_corr) # Save space in memory - # Reshape exp_corr data onto time dimension before 'Split' + ## Reshape exp_corr data onto time dimension before 'Split' rep_pos <- array(NA,c(time=length(wt_exp2))) pos_time <- which(names(dim(res)) == 'time') pos_type <- which(names(dim(res)) == 'type') for (x in unique(wt_exp2)){ rep_pos[which(wt_exp2==x)]<-1:length(which(wt_exp2==x)) } - exp_corr <- .unsplit_wtype(exp=res,wt_exp=wt_exp2,rep_pos=rep_pos,pos_time=pos_time) + exp_corr <- .unsplit_wtype(exp=res,wt_exp=wt_exp2,rep_pos=rep_pos, + pos_time=pos_time) # Now reshape exp_corr data onto original dimensions dim(exp_corr) <- c(dim(wt_exp), dim(exp_corr)[-c(pos_time,pos_type)]) return(exp_corr) @@ -161,12 +278,15 @@ AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, corrdims = c('member', 'sdat return(exp[NN$imin_lat[ilat,ilon],NN$imin_lon[ilat,ilon]]) } -.unsplit_wtype <- function(exp=exp,dim_wt='type',wt_exp=wt_exp,dim_time='time',rep_pos=rep_pos,pos_time=1){ +.unsplit_wtype <- function(exp=exp,dim_wt='type',wt_exp=wt_exp, + dim_time='time',rep_pos=rep_pos,pos_time=1){ # Initiate output - new <- Subset(Subset(exp, along=dim_wt, indices=wt_exp[1]), along=dim_time, indices=rep_pos[1]) + new <- Subset(Subset(exp, along=dim_wt, indices=wt_exp[1]), along=dim_time, + indices=rep_pos[1]) dimnames <- names(dim(new)) for (x in 2:length(wt_exp)){ - dat <- Subset(Subset(exp, along=dim_wt, indices=wt_exp[x]), along=dim_time, indices=rep_pos[x]) + dat <- Subset(Subset(exp, along=dim_wt, indices=wt_exp[x]), + along=dim_time, indices=rep_pos[x]) new <- abind::abind(new,dat,along=pos_time) } names(dim(new)) <- dimnames diff --git a/man/AdamontAnalog.Rd b/man/AdamontAnalog.Rd index 576f9ce0..728ed9bc 100644 --- a/man/AdamontAnalog.Rd +++ b/man/AdamontAnalog.Rd @@ -2,8 +2,8 @@ % Please edit documentation in R/CST_AdamontAnalog.R \name{AdamontAnalog} \alias{AdamontAnalog} -\title{AdamontAnalog finds analogous data in the reference dataset to experiment data based on -weather types} +\title{AdamontAnalog finds analogous data in the reference dataset to experiment data +based on weather types} \usage{ AdamontAnalog( exp, @@ -19,34 +19,48 @@ AdamontAnalog( ) } \arguments{ -\item{exp}{experiment data, such as \code{$data} array of an object of class \code{s2dv_cube}, can be output from quantile correction using CST_AdamontQQCorr.R} +\item{exp}{experiment data, such as \code{$data} array of an object of class +\code{s2dv_cube}, can be output from quantile correction using +\code{CST_AdamontQQCorr.R}} -\item{obs}{reference data, can also be \code{$data} array of class \code{s2dv_cube}. Note that lat/lon dimensions need to be the same as \code{exp}} +\item{obs}{reference data, can also be \code{$data} array of class +\code{s2dv_cube}. Note that lat/lon dimensions need to be the same as +\code{exp}} -\item{wt_exp}{corresponding weather types (same dimensions as \code{exp} but lat/lon)} +\item{wt_exp}{corresponding weather types (same dimensions as \code{exp} but +lat/lon)} -\item{wt_obs}{corresponding weather types (same dimensions as \code{obs} but lat/lon)} +\item{wt_obs}{corresponding weather types (same dimensions as \code{obs} but * +lat/lon)} -\item{nanalogs}{integer defining the number of analog values to return (default: 5)} +\item{nanalogs}{integer defining the number of analog values to return +(default: 5)} -\item{method}{a character string indicating the method used for analog definition -Coded are 'pattcorr': pattern correlation - 'rain1' (for precip patterns): rain occurrence consistency - 'rain01' (for precip patterns): rain occurrence/non occurrence consistency} +\item{method}{a character string indicating the method used for analog +definition + Coded are 'pattcorr': pattern correlation + 'rain1' (for precip patterns): rain occurrence consistency + 'rain01' (for precip patterns): rain occurrence/non + occurrence consistency} -\item{thres}{real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1} +\item{thres}{real number indicating the threshold to define rain occurrence +/non occurrence in rain(0)1} -\item{search_obsdims}{list of dimensions in \code{obs} along which analogs are searched for} +\item{search_obsdims}{list of dimensions in \code{obs} along which analogs are +searched for} \item{londim}{name of longitude dimension} \item{latdim}{name of latitude dimension} } \value{ -analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data +analog_vals an object of class \code{s2dv_cube} containing nanalogs +analog values for each value of \code{exp} input data } \description{ -This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data. +This function searches for analogs in a reference dataset for +experiment data, based on corresponding weather types. The experiment data is +typically a hindcast, observations are typically provided by reanalysis data. } \examples{ \dontrun{ @@ -54,6 +68,7 @@ wt_exp <- sample(1:3, 15*6*3, replace=T) dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) wt_obs <- sample(1:3, 6*3, replace=T) dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) + obs=lonlat_data$obs$data, wt_exp=wt_exp, wt_obs=wt_obs, nanalogs=2) } } \author{ diff --git a/man/AdamontQQCorr.Rd b/man/AdamontQQCorr.Rd index d37a7086..2006a7ad 100644 --- a/man/AdamontQQCorr.Rd +++ b/man/AdamontQQCorr.Rd @@ -2,7 +2,8 @@ % Please edit documentation in R/CST_AdamontQQCorr.R \name{AdamontQQCorr} \alias{AdamontQQCorr} -\title{AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types} +\title{AdamontQQCorr computes quantile-quantile correction of seasonal or decadal +forecast data using weather types} \usage{ AdamontQQCorr( exp, @@ -17,24 +18,43 @@ AdamontQQCorr( ) } \arguments{ -\item{exp}{array with named dimensions (such as \code{$data} array of experiment data from an object of class \code{s2dv_cube})} +\item{exp}{array with named dimensions (such as \code{$data} array of +experiment data from an object of class \code{s2dv_cube})} -\item{wt_exp}{corresponding weather types (same dimensions as \code{exp} but lat/lon)} +\item{wt_exp}{corresponding weather types (same dimensions as \code{exp} but +lat/lon)} -\item{obs}{array with named dimensions with reference data (can also be \code{$data} array of class \code{s2dv_cube}). lat/lon dimensions can differ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and AdamontNearestNeighbors \code{NN} output should be provided} +\item{obs}{array with named dimensions with reference data (can also be +\code{$data} array of class \code{s2dv_cube}). lat/lon dimensions can differ +from \code{exp} if non rectilinear latlon grids are used, in which case +regrid should be set to TRUE and .NearestNeighbors \code{NN} output should be +provided} -\item{corrdims}{list of dimensions in \code{exp} for which quantile mapping correction is applied} +\item{corrdims}{list of dimensions in \code{exp} for which quantile mapping +correction is applied} -\item{londim}{character name of longitude dimension in \code{exp} and \code{obs}} +\item{londim}{character name of longitude dimension in \code{exp} and +\code{obs}} -\item{latdim}{character name of latitude dimension in \code{exp} and \code{obs}} +\item{latdim}{character name of latitude dimension in \code{exp} and +\code{obs}} -\item{regrid}{boolean indicating whether Nearest Neighbors regridding is needed} +\item{regrid}{(optional) boolean indicating whether .NearestNeighbors +regridding is needed} -\item{NN}{list (output from AdamontNearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o)} +\item{NN}{(optional, if regrid=TRUE) list (output from .NearestNeighbors) +maps (nlat, nlon) onto (nlat_o, nlon_o)} +} +\value{ +an array (such as \code{$data} array from an object of class +\code{s2dv_cube}) with named dimensions, containing experiment data on the +lat/lon grid of \code{obs} array, corrected by quantile mapping depending on +the weather types \code{wt_exp} } \description{ -This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference \code{obs}, typically provided by reanalysis data. +This function computes a quantile mapping based on weather types +for experiment data (typically a hindcast) onto reference \code{obs}, +typically provided by reanalysis data. } \examples{ \dontrun{ @@ -42,7 +62,9 @@ wt_exp <- sample(1:3, 15*6*3, replace=T) dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) wt_obs <- sample(1:3, 6*3, replace=T) dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) -exp_corr <- AdamontQQCorr(exp=lonlat_data$exp$data, wt_exp=wt_exp, obs=lonlat_data$obs$data, wt_obs=wt_obs, corrdims = c('dataset','member','sdate','ftime')) +exp_corr <- AdamontQQCorr(exp=lonlat_data$exp$data, wt_exp=wt_exp, + obs=lonlat_data$obs$data, wt_obs=wt_obs, + corrdims = c('dataset','member','sdate','ftime')) } } \author{ diff --git a/man/CST_AdamontAnalog.Rd b/man/CST_AdamontAnalog.Rd index 579e8304..7d10ae02 100644 --- a/man/CST_AdamontAnalog.Rd +++ b/man/CST_AdamontAnalog.Rd @@ -2,8 +2,8 @@ % Please edit documentation in R/CST_AdamontAnalog.R \name{CST_AdamontAnalog} \alias{CST_AdamontAnalog} -\title{CST_AdamontAnalog finds analogous data in the reference dataset to experiment data based on -weather types} +\title{CST_AdamontAnalog finds analogous data in the reference dataset to experiment +data based on weather types} \usage{ CST_AdamontAnalog( exp, @@ -19,34 +19,46 @@ CST_AdamontAnalog( ) } \arguments{ -\item{exp}{experiment data an object of class \code{s2dv_cube}, can be output from quantile correction using CST_AdamontQQCorr.R} +\item{exp}{experiment data an object of class \code{s2dv_cube}, can be output +from quantile correction using CST_AdamontQQCorr.R} -\item{obs}{reference data, also of class \code{s2dv_cube}. Note that lat/lon dimensions need to be the same as \code{exp}} +\item{obs}{reference data, also of class \code{s2dv_cube}. Note that lat/lon +dimensions need to be the same as \code{exp}} -\item{wt_exp}{corresponding weather types (same dimensions as \code{exp$data} but lat/lon)} +\item{wt_exp}{corresponding weather types (same dimensions as \code{exp$data} +but lat/lon)} -\item{wt_obs}{corresponding weather types (same dimensions as \code{obs$data} but lat/lon)} +\item{wt_obs}{corresponding weather types (same dimensions as \code{obs$data} +but lat/lon)} -\item{nanalogs}{integer defining the number of analog values to return (default: 5)} +\item{nanalogs}{integer defining the number of analog values to return +(default: 5)} -\item{method}{a character string indicating the method used for analog definition -Coded are 'pattcorr': pattern correlation - 'rain1' (for precip patterns): rain occurrence consistency - 'rain01' (for precip patterns): rain occurrence/non occurrence consistency} +\item{method}{a character string indicating the method used for analog +definition + Coded are 'pattcorr': pattern correlation + 'rain1' (for precip patterns): rain occurrence consistency + 'rain01' (for precip patterns): rain occurrence/non + occurrence consistency} -\item{thres}{real number indicating the threshold to define rain occurrence/non occurrence in rain(0)1} +\item{thres}{real number indicating the threshold to define rain +occurrence/non occurrence in rain(0)1} -\item{search_obsdims}{list of dimensions in \code{obs} along which analogs are searched for} +\item{search_obsdims}{list of dimensions in \code{obs} along which analogs are +searched for} \item{londim}{name of longitude dimension} \item{latdim}{name of latitude dimension} } \value{ -analog_vals an object of class \code{s2dv_cube} containing nanalogs analog values for each value of \code{exp} input data +analog_vals an object of class \code{s2dv_cube} containing nanalogs +analog values for each value of \code{exp} input data } \description{ -This function searches for analogs in a reference dataset for experiment data, based on corresponding weather types. The experiment data is typically a hindcast, observations are typically provided by reanalysis data. +This function searches for analogs in a reference dataset for +experiment data, based on corresponding weather types. The experiment data is +typically a hindcast, observations are typically provided by reanalysis data. } \examples{ \dontrun{ diff --git a/man/CST_AdamontQQCorr.Rd b/man/CST_AdamontQQCorr.Rd index dcd31017..fc4b1932 100644 --- a/man/CST_AdamontQQCorr.Rd +++ b/man/CST_AdamontQQCorr.Rd @@ -2,7 +2,8 @@ % Please edit documentation in R/CST_AdamontQQCorr.R \name{CST_AdamontQQCorr} \alias{CST_AdamontQQCorr} -\title{CST_AdamontQQCorr computes quantile-quantile correction of seasonal or decadal forecast data using weather types} +\title{CST_AdamontQQCorr computes quantile-quantile correction of seasonal or +decadal forecast data using weather types} \usage{ CST_AdamontQQCorr( exp, @@ -17,22 +18,38 @@ CST_AdamontQQCorr( \arguments{ \item{exp}{experiment data an object of class \code{s2dv_cube}} -\item{wt_exp}{corresponding weather types (same dimensions as \code{exp$data} but lat/lon)} +\item{wt_exp}{corresponding weather types (same dimensions as \code{exp$data} +but lat/lon)} -\item{obs}{reference data, also of class \code{s2dv_cube}. lat/lon dimensions can differ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and AdamontNearestNeighbors \code{NN} output should be provided} +\item{obs}{reference data, also of class \code{s2dv_cube}. lat/lon dimensions +can differ from \code{exp} if non rectilinear latlon grids are used, +in which case regrid should be set to TRUE and .NearestNeighbors \code{NN} +output should be provided} -\item{corrdims}{list of dimensions in \code{exp} for which quantile mapping correction is applied} +\item{corrdims}{list of dimensions in \code{exp} for which quantile mapping +correction is applied} -\item{londim}{character name of longitude dimension in \code{exp} and \code{obs}} +\item{londim}{character name of longitude dimension in \code{exp} and +\code{obs}} -\item{latdim}{character name of latitude dimension in \code{exp} and \code{obs}} +\item{latdim}{character name of latitude dimension in \code{exp} and +\code{obs}} -\item{regrid}{boolean indicating whether Nearest Neighbors regridding is needed} +\item{regrid}{(optional) boolean indicating whether .NearestNeighbors +regridding is needed} -\item{NN}{list (output from Adamont_NearestNeighbors) maps (nlat, nlon) onto (nlat_o, nlon_o)} +\item{NN}{list (optional, if regrid=TRUE) output from .NearestNeighbors +mapping (nlat, nlon) onto (nlat_o, nlon_o)} +} +\value{ +an object of class \code{s2dv_cube} containing experiment data on the + lat/lon grid of \code{obs} input data, corrected by quantile mapping + depending on the weather types \code{wt_exp} } \description{ -This function computes a quantile mapping based on weather types for experiment data (typically a hindcast) onto reference \code{obs}, typically provided by reanalysis data. +This function computes a quantile mapping based on weather types +for experiment data (typically a hindcast) onto reference \code{obs}, +typically provided by reanalysis data. } \examples{ \dontrun{ @@ -40,11 +57,15 @@ wt_exp <- sample(1:3, 15*6*3, replace=T) dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3) wt_obs <- sample(1:3, 6*3, replace=T) dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3) -exp_corr <- CST_AdamontQQCorr(exp=lonlat_data$exp, wt_exp=wt_exp, obs=lonlat_data$obs, wt_obs=wt_obs, corrdims = c('dataset','member','sdate','ftime')) +exp_corr <- CST_AdamontQQCorr(exp=lonlat_data$exp, wt_exp=wt_exp, + obs=lonlat_data$obs, wt_obs=wt_obs, + corrdims = c('dataset','member','sdate','ftime')) } } \author{ -Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +Lauriane Batté, \email{lauriane.batte@meteo.fr} + +Paola Marson, \email{paola.marson@meteo.fr} -Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +Gildas Dayon, \email{gildas.dayon@meteo.fr} } -- GitLab From 1890eee3b7cbea10bea2036c4d0163daac6f556e Mon Sep 17 00:00:00 2001 From: lbatteCNRM Date: Fri, 19 Feb 2021 13:59:08 +0100 Subject: [PATCH 22/24] Fixes to sanity checks --- R/CST_AdamontAnalog.R | 12 ++++++++---- R/CST_AdamontQQCorr.R | 23 +++++++++++------------ man/NearestNeighbors.Rd | 30 ------------------------------ 3 files changed, 19 insertions(+), 46 deletions(-) delete mode 100644 man/NearestNeighbors.Rd diff --git a/R/CST_AdamontAnalog.R b/R/CST_AdamontAnalog.R index 83f4795b..18350806 100644 --- a/R/CST_AdamontAnalog.R +++ b/R/CST_AdamontAnalog.R @@ -58,9 +58,13 @@ CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, if (is.null(nanalogs)){ nanalogs <- 5 } - if (is.null(which(dimnames == latdim)) || is.null(which(dimnames == londim))){ - stop("'londim' or 'latdim' inputs should correspond to 'obs$data' ", - "dimension names") + if (!(latdim %in% dimnames) || !(londim %in% dimnames)){ + stop("'londim' or 'latdim' input doesn't match with 'obs$data' dimension", + " names") + } + if (!(latdim %in% dimnamesexp) || !(londim %in% dimnamesexp)){ + stop("'londim' or 'latdim' input doesn't match with 'exp$data' dimension", + " names") } if (!all(search_obsdims %in% dimnames)){ stop("Names in parameter 'search_obsdims' should match 'obs$data' ", @@ -154,7 +158,7 @@ AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs=5, if (method == 'pattcorr' & !is.null(thres)){ warning("Parameter 'thres' is not used with method 'pattcorr'.") } - if (is.null(which(dimnamesexp == latdim)) || is.null(which(dimnamesexp==londim))){ + if (!(latdim %in% dimnamesexp) || !(londim %in% dimnamesexp)){ stop("'londim' or 'latdim' input doesn't match with 'exp' dimension names") } # Position of lat/lon dimensions in exp data diff --git a/R/CST_AdamontQQCorr.R b/R/CST_AdamontQQCorr.R index c8719a88..fe13feb1 100644 --- a/R/CST_AdamontQQCorr.R +++ b/R/CST_AdamontQQCorr.R @@ -54,16 +54,15 @@ CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, } dimnames <- names(dim(obs$data)) dimnamesexp <- names(dim(exp$data)) - if (is.null(which(dimnames==latdim)) || is.null(which(dimnames==londim))){ - stop("'londim' or 'latdim' input doesn't match with 'obs$data' - dimension names") + if (!(latdim %in% dimnames) || !(londim %in% dimnames)){ + stop("'londim' or 'latdim' input doesn't match with 'obs$data' dimension", + " names") } - if (is.null(which(dimnamesexp==latdim)) || - is.null(which(dimnamesexp==londim))){ - stop("'londim' or 'latdim' input doesn't match with 'exp$data' - dimension names") + if (!(latdim %in% dimnamesexp) || !(londim %in% dimnamesexp)){ + stop("'londim' or 'latdim' input doesn't match with 'exp$data' dimension", + " names") } - if (is.null(which(corrdims=='time')) || is.null(which(corrdims=='ftime'))){ + if (!(('time' %in% corrdims) || ('ftime' %in% corrdims))){ warning("Forecast time should be one of the dimensions for the correction specified in corrdims input list") } @@ -156,12 +155,12 @@ AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, dimnames <- names(dim(obs)) dimnamesexp <- names(dim(exp)) - if (is.null(which(dimnames == latdim)) || is.null(which(dimnames==londim))){ + if (!(latdim %in% dimnames) || !(londim %in% dimnames)){ stop("'londim' or 'latdim' input doesn't match with 'obs' dimension names") } - if (is.null(which(corrdims=='time')) || is.null(which(corrdims=='ftime'))){ - warning("Forecast time should be one of the dimensions for the correction - specified in corrdims input list") + if (!(('time' %in% corrdims) || ('ftime' %in% corrdims))){ + warning("Forecast time should be one of the dimensions for the correction", + " specified in corrdims input list") } if (!all(corrdims %in% dimnamesexp)){ stop("Names in parameter 'corrdims' should match input dimension names.") diff --git a/man/NearestNeighbors.Rd b/man/NearestNeighbors.Rd deleted file mode 100644 index 11e6b2ff..00000000 --- a/man/NearestNeighbors.Rd +++ /dev/null @@ -1,30 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AdamontNearestNeighbors.R -\name{NearestNeighbors} -\alias{NearestNeighbors} -\title{ADAMONT Nearest Neighbors computes the distance between reference data grid centroid and SF data grid} -\usage{ -NearestNeighbors(exp, obs, method = "ADA") -} -\arguments{ -\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment longitudes in \code{$lon} and latitudes in \code{$lat}} - -\item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the reference data on a different grid, with longitudes in \code{$lon} and latitudes in \code{$lat}.} - -\item{method}{a string among three options ('ADA': standard ADAMONT distance, 'simple': lon/lat straight euclidian distance, 'radius': distance on the sphere)} -} -\value{ -NN a list, containing the following: - min_lon: array of dimensions \code{obs$lon} giving the longitude of closest gridpoint in exp - min_lat: array of dimensions \code{obs$lat} giving the latitude of closest gridpoint in exp - imin_lon: array of dimensions \code{obs$lon} giving the longitude index of closest gridpoint in exp - imin_lat: array of dimensions \code{obs$lat} giving the latitude index of closest gridpoint in exp -} -\description{ -This function computes the nearest neighbor for each reference data (lon, lat) point in the experiment dataset by computing the distance between the reference dataset grid and the experiment data. This is the first step in the ADAMONT method adapted from Verfaillie et al. (2018). -} -\author{ -Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version - -Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation -} -- GitLab From d1cd32d830f823d29a01bb56902f60d3f4aec31a Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 19 Feb 2021 18:44:52 +0100 Subject: [PATCH 23/24] fixing documentation and code --- DESCRIPTION | 2 +- NAMESPACE | 1 + R/AdamontNearestNeighbors.R | 101 ------------------------------- R/CST_AdamontAnalog.R | 4 +- R/CST_AdamontQQCorr.R | 115 +++++++++++++++++++++++++++++++++--- R/CST_SplitDim.R | 79 +++++++++++++++++++++---- man/AdamontQQCorr.Rd | 3 + man/CST_AdamontQQCorr.Rd | 9 +-- man/CST_SplitDim.Rd | 16 ++++- man/SplitDim.Rd | 12 +++- 10 files changed, 211 insertions(+), 131 deletions(-) delete mode 100644 R/AdamontNearestNeighbors.R diff --git a/DESCRIPTION b/DESCRIPTION index 12fe8f75..a2c5b370 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -81,4 +81,4 @@ VignetteBuilder: knitr License: Apache License 2.0 Encoding: UTF-8 LazyData: true -RoxygenNote: 7.1.1 +RoxygenNote: 7.0.2 diff --git a/NAMESPACE b/NAMESPACE index d8ba17fa..659a4821 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -86,6 +86,7 @@ importFrom(plyr,.) importFrom(plyr,dlply) importFrom(reshape2,melt) importFrom(s2dv,Reorder) +importFrom(s2dverification,Subset) importFrom(utils,glob2rx) importFrom(utils,head) importFrom(utils,tail) diff --git a/R/AdamontNearestNeighbors.R b/R/AdamontNearestNeighbors.R deleted file mode 100644 index 4c6dbe07..00000000 --- a/R/AdamontNearestNeighbors.R +++ /dev/null @@ -1,101 +0,0 @@ -#' ADAMONT Nearest Neighbors computes the distance between reference data grid centroid and SF data grid -#' -#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version -#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation -#'@description This function computes the nearest neighbor for each reference data (lon, lat) point in the experiment dataset by computing the distance between the reference dataset grid and the experiment data. This is the first step in the ADAMONT method adapted from Verfaillie et al. (2018). -#' -#'@param method a string among three options ('ADA': standard ADAMONT distance, 'simple': lon/lat straight euclidian distance, 'radius': distance on the sphere) -#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment longitudes in \code{$lon} and latitudes in \code{$lat} -#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the reference data on a different grid, with longitudes in \code{$lon} and latitudes in \code{$lat}. -#' -#'@return NN a list, containing the following: -#' min_lon: array of dimensions \code{obs$lon} giving the longitude of closest gridpoint in exp -#' min_lat: array of dimensions \code{obs$lat} giving the latitude of closest gridpoint in exp -#' imin_lon: array of dimensions \code{obs$lon} giving the longitude index of closest gridpoint in exp -#' imin_lat: array of dimensions \code{obs$lat} giving the latitude index of closest gridpoint in exp -#' -#'@import s2dverification -#'@import ncdf4 - -.NearestNeighbors <- function (exp, obs, method='ADA') { - - if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { - stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ", - "as output by CSTools::CST_Load.") - } - exp_lon <- exp$lon - exp_lat <- exp$lat - obs_lon <- obs$lon - obs_lat <- obs$lat - dim_exp_lon <- dim(exp_lon) - dim_exp_lat <- dim(exp_lat) - dim_obs_lon <- dim(obs_lon) - dim_obs_lat <- dim(obs_lat) - # Check if one of the grids is non-regular: - if ((length(dim_exp_lon)==2) || (length(dim_obs_lon)==2)){ - # Flatten longitudes and latitudes in case of 2-D longitudes and latitudes (Lambert grids, etc.) - if ((length(dim_exp_lon)==2) & (length(dim_exp_lat)==2)){ - dim(exp_lon) <- c(dim_exp_lon[1]*dim_exp_lon[2]) - dim(exp_lat) <- c(dim_exp_lat[1]*dim_exp_lat[2]) - } - if ((length(dim_obs_lon)==2) & (length(dim_obs_lat)==2)){ - dim(obs_lon) <- c(dim_obs_lon[1]*dim_obs_lon[2]) - dim(obs_lat) <- c(dim_obs_lat[1]*dim_obs_lat[2]) - } - # Now lat and lon arrays have 1 dimension, length npt (= nlat*nlon) - OBS_grid <- cbind(obs_lon,obs_lat) - EXP_grid <- cbind(exp_lon,exp_lat) - dist_min<-min_lon<-min_lat<-imin_lon<-imin_lat<-array(dim=nrow(OBS_grid)) - if (method == 'ADA'){ - C<-cos(OBS_grid[,2]*pi/180)^2 - for (i in 1:nrow(OBS_grid)){ - dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+C[i]*(OBS_grid[i,1]-EXP_grid[,1])^2 - dist_min[i]<-min(dist) - min_lon[i]<-EXP_grid[which.min(dist),1] - min_lat[i]<-EXP_grid[which.min(dist),2] - imin_lon[i]<-which(exp_lon==min_lon[i]) - imin_lat[i]<-which(exp_lat==min_lat[i]) - } - } else if (method == 'simple'){ - for (i in 1:nrow(OBS_grid)){ - dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+(OBS_grid[i,1]-EXP_grid[,1])^2 - dist_min[i]<-min(dist) - min_lon[i]<-EXP_grid[which.min(dist),1] - min_lat[i]<-EXP_grid[which.min(dist),2] - imin_lon[i]<-which(exp_lon==min_lon[i]) - imin_lat[i]<-which(exp_lat==min_lat[i]) - } - } else if (method == 'radius'){ - R <- 6371e3 # metres, Earth radius - EXP_gridr<-EXP_grid*pi/180 - OBS_gridr<-OBS_grid*pi/180 - for (i in 1:nrow(OBS_grid)){ - a<-sin((OBS_gridr[i,2]-EXP_gridr[,2])/2)^2 + cos(OBS_gridr[i,2])*cos(EXP_gridr[,2])*sin((OBS_gridr[i,1]-EXP_gridr[,1])/2)^2 - c<-2*atan2(sqrt(a),sqrt(1-a)) - dist<-R*c - dist_min[i]<-min(dist) - min_lon[i]<-EXP_grid[which.min(dist),1] - min_lat[i]<-EXP_grid[which.min(dist),2] - imin_lon[i]<-which(exp_lon==min_lon[i]) - imin_lat[i]<-which(exp_lat==min_lat[i]) - } - } else { - stop("AdamontNearestNeighbors supports method = 'ADA', 'simple' or 'radius' only.") - } - - # Reshape outputs to original grid - dim(min_lon)=dim_obs_lon - dim(min_lat)=dim_obs_lat - dim(imin_lon)=dim_obs_lon - dim(imin_lat)=dim_obs_lat - - } else { - # Regular lon/lat grid case: has been handled by CST_Load() - stop("AdamontNearestNeighbors is meant for non-regular lat/lon grids; use e.g. CST_Load to interpolate exp onto obs grid") - } - - NN=list(min_lon=min_lon, min_lat=min_lat, imin_lon=imin_lon, imin_lat=imin_lat) - - return(NN) -} - diff --git a/R/CST_AdamontAnalog.R b/R/CST_AdamontAnalog.R index 18350806..01e62d97 100644 --- a/R/CST_AdamontAnalog.R +++ b/R/CST_AdamontAnalog.R @@ -186,9 +186,9 @@ AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs=5, if (length(search_obsdims) > 1) { for (i in 1:(length(search_obsdims) - 1)) { obs <- MergeDims(obs, search_obsdims[i:(i + 1)], - rename = search_obsdims[i + 1]) + rename_dim = search_obsdims[i + 1]) wt_obs <- MergeDims(wt_obs, search_obsdims[i:(i + 1)], - rename = search_obsdims[i + 1]) + rename_dim = search_obsdims[i + 1]) } } names(dim(obs))[which(names(dim(obs)) == search_obsdims[length(search_obsdims)])] <- 'time' diff --git a/R/CST_AdamontQQCorr.R b/R/CST_AdamontQQCorr.R index fe13feb1..be1dbbc8 100644 --- a/R/CST_AdamontQQCorr.R +++ b/R/CST_AdamontQQCorr.R @@ -15,16 +15,14 @@ #' can differ from \code{exp} if non rectilinear latlon grids are used, #' in which case regrid should be set to TRUE and .NearestNeighbors \code{NN} #' output should be provided +#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but +#'lat/lon) #'@param corrdims list of dimensions in \code{exp} for which quantile mapping #' correction is applied #'@param londim character name of longitude dimension in \code{exp} and #' \code{obs} #'@param latdim character name of latitude dimension in \code{exp} and #' \code{obs} -#'@param regrid (optional) boolean indicating whether .NearestNeighbors -#' regridding is needed -#'@param NN list (optional, if regrid=TRUE) output from .NearestNeighbors -#' mapping (nlat, nlon) onto (nlat_o, nlon_o) #' #'@return an object of class \code{s2dv_cube} containing experiment data on the #' lat/lon grid of \code{obs} input data, corrected by quantile mapping @@ -119,6 +117,8 @@ CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, #'from \code{exp} if non rectilinear latlon grids are used, in which case #'regrid should be set to TRUE and .NearestNeighbors \code{NN} output should be #'provided +#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but +#'lat/lon) #'@param corrdims list of dimensions in \code{exp} for which quantile mapping #'correction is applied #'@param londim character name of longitude dimension in \code{exp} and @@ -218,10 +218,10 @@ AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, ## Use CST_QuantileMapping function for quantile mapping ## depending on weather type for (i in 1:(length(corrdims) - 1)) { - obs <- MergeDims(obs, corrdims[i:(i+1)], rename=corrdims[i+1]) - wt_obs <- MergeDims(wt_obs, corrdims[i:(i+1)], rename=corrdims[i+1]) - exp_corr <- MergeDims(exp_corr, corrdims[i:(i+1)], rename=corrdims[i+1]) - wt_exp2 <- MergeDims(wt_exp2, corrdims[i:(i+1)], rename=corrdims[i+1]) + obs <- MergeDims(obs, corrdims[i:(i+1)], rename_dim=corrdims[i+1]) + wt_obs <- MergeDims(wt_obs, corrdims[i:(i+1)], rename_dim=corrdims[i+1]) + exp_corr <- MergeDims(exp_corr, corrdims[i:(i+1)], rename_dim=corrdims[i+1]) + wt_exp2 <- MergeDims(wt_exp2, corrdims[i:(i+1)], rename_dim=corrdims[i+1]) } names(dim(obs))[which(names(dim(obs)) == corrdims[length(corrdims)])] <- 'time' names(dim(wt_obs))[which(names(dim(wt_obs)) == corrdims[length(corrdims)])] <- 'time' @@ -291,4 +291,103 @@ AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, names(dim(new)) <- dimnames return(new) } +#' ADAMONT Nearest Neighbors computes the distance between reference data grid centroid and SF data grid +#' +#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version +#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation +#'@description This function computes the nearest neighbor for each reference data (lon, lat) point in the experiment dataset by computing the distance between the reference dataset grid and the experiment data. This is the first step in the ADAMONT method adapted from Verfaillie et al. (2018). +#' +#'@param method a string among three options ('ADA': standard ADAMONT distance, 'simple': lon/lat straight euclidian distance, 'radius': distance on the sphere) +#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment longitudes in \code{$lon} and latitudes in \code{$lat} +#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the reference data on a different grid, with longitudes in \code{$lon} and latitudes in \code{$lat}. +#' +#'@return NN a list, containing the following: +#' min_lon: array of dimensions \code{obs$lon} giving the longitude of closest gridpoint in exp +#' min_lat: array of dimensions \code{obs$lat} giving the latitude of closest gridpoint in exp +#' imin_lon: array of dimensions \code{obs$lon} giving the longitude index of closest gridpoint in exp +#' imin_lat: array of dimensions \code{obs$lat} giving the latitude index of closest gridpoint in exp +#' +#'@import s2dverification +#'@import ncdf4 +#'@noRd +.NearestNeighbors <- function (exp, obs, method='ADA') { + + if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { + stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + exp_lon <- exp$lon + exp_lat <- exp$lat + obs_lon <- obs$lon + obs_lat <- obs$lat + dim_exp_lon <- dim(exp_lon) + dim_exp_lat <- dim(exp_lat) + dim_obs_lon <- dim(obs_lon) + dim_obs_lat <- dim(obs_lat) + # Check if one of the grids is non-regular: + if ((length(dim_exp_lon)==2) || (length(dim_obs_lon)==2)){ + # Flatten longitudes and latitudes in case of 2-D longitudes and latitudes (Lambert grids, etc.) + if ((length(dim_exp_lon)==2) & (length(dim_exp_lat)==2)){ + dim(exp_lon) <- c(dim_exp_lon[1]*dim_exp_lon[2]) + dim(exp_lat) <- c(dim_exp_lat[1]*dim_exp_lat[2]) + } + if ((length(dim_obs_lon)==2) & (length(dim_obs_lat)==2)){ + dim(obs_lon) <- c(dim_obs_lon[1]*dim_obs_lon[2]) + dim(obs_lat) <- c(dim_obs_lat[1]*dim_obs_lat[2]) + } + # Now lat and lon arrays have 1 dimension, length npt (= nlat*nlon) + OBS_grid <- cbind(obs_lon,obs_lat) + EXP_grid <- cbind(exp_lon,exp_lat) + dist_min<-min_lon<-min_lat<-imin_lon<-imin_lat<-array(dim=nrow(OBS_grid)) + if (method == 'ADA'){ + C<-cos(OBS_grid[,2]*pi/180)^2 + for (i in 1:nrow(OBS_grid)){ + dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+C[i]*(OBS_grid[i,1]-EXP_grid[,1])^2 + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else if (method == 'simple'){ + for (i in 1:nrow(OBS_grid)){ + dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+(OBS_grid[i,1]-EXP_grid[,1])^2 + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else if (method == 'radius'){ + R <- 6371e3 # metres, Earth radius + EXP_gridr<-EXP_grid*pi/180 + OBS_gridr<-OBS_grid*pi/180 + for (i in 1:nrow(OBS_grid)){ + a<-sin((OBS_gridr[i,2]-EXP_gridr[,2])/2)^2 + cos(OBS_gridr[i,2])*cos(EXP_gridr[,2])*sin((OBS_gridr[i,1]-EXP_gridr[,1])/2)^2 + c<-2*atan2(sqrt(a),sqrt(1-a)) + dist<-R*c + dist_min[i]<-min(dist) + min_lon[i]<-EXP_grid[which.min(dist),1] + min_lat[i]<-EXP_grid[which.min(dist),2] + imin_lon[i]<-which(exp_lon==min_lon[i]) + imin_lat[i]<-which(exp_lat==min_lat[i]) + } + } else { + stop("AdamontNearestNeighbors supports method = 'ADA', 'simple' or 'radius' only.") + } + + # Reshape outputs to original grid + dim(min_lon)=dim_obs_lon + dim(min_lat)=dim_obs_lat + dim(imin_lon)=dim_obs_lon + dim(imin_lat)=dim_obs_lat + } else { + # Regular lon/lat grid case: has been handled by CST_Load() + stop("AdamontNearestNeighbors is meant for non-regular lat/lon grids; use e.g. CST_Load to interpolate exp onto obs grid") + } + + NN=list(min_lon=min_lon, min_lat=min_lat, imin_lon=imin_lon, imin_lat=imin_lat) + + return(NN) +} diff --git a/R/CST_SplitDim.R b/R/CST_SplitDim.R index 46cd97cc..5e85182d 100644 --- a/R/CST_SplitDim.R +++ b/R/CST_SplitDim.R @@ -8,9 +8,12 @@ #'@param split_dim a character string indicating the name of the dimension to split #'@param indices a vector of numeric indices or dates. If left at NULL, the dates provided in the s2dv_cube object (element Dates) will be used. #'@param freq a character string indicating the frequency: by 'day', 'month' and 'year' or 'monthly' (by default). 'month' identifies months between 1 and 12 independently of the year they belong to, while 'monthly' differenciates months from different years. +#'@param new_dim_name a character string indicating the name of the new dimension. +#'@param insert_ftime an integer indicating the number of time steps to add at the begining of the time series. #' +#'@details Parameter 'insert_ftime' has been included for the case of using daily data, requiring split the temporal dimensions by months (or similar) and the first lead time doesn't correspondt to the 1st day of the month. In this case, the insert_ftime could be used, to get a final output correctly organized. E.g.: leadtime 1 is the 2nd of November and the input time series extend to the 31st of December. When requiring split by month with \code{inset_ftime = 1}, the 'monthly' dimension of length two will indicate the month (position 1 for November and position 2 for December), dimension 'time' will be length 31. For November, the position 1 and 31 will be NAs, while from positon 2 to 30 will be filled with the data provided. This allows to select correctly days trhough time dimension. #'@import abind -#'@import s2dverification +#'@importFrom s2dverification Subset #'@examples #' #'data <- 1 : 20 @@ -34,12 +37,46 @@ #'dim(new_data$data) #'@export CST_SplitDim <- function(data, split_dim = 'time', indices = NULL, - freq = 'monthly') { + freq = 'monthly', new_dim_name = NULL, insert_ftime = NULL) { if (!inherits(data, 's2dv_cube')) { - stop("Parameter 'data' must be of the class 's2dv_cube', ", - "as output by CSTools::CST_Load.") + stop("Parameter 'data' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") } + if (!is.null(insert_ftime)) { + if (!is.numeric(insert_ftime)) { + stop("Parameter 'insert_ftime' should be an integer.") + } else { + if (length(insert_ftime) > 1) { + warning("Parameter 'insert_ftime' must be of length 1, and only the", + " first element will be used.") + insert_ftime <- insert_ftime[1] + } + # adding NAs at the begining of the data in ftime dim + ftimedim <- which(names(dim(data$data)) == 'ftime') + dims <- dim(data$data) + dims[ftimedim] <- insert_ftime + empty_array <- array(NA, dims) + data$data <- abind(empty_array, data$data, along = ftimedim) + names(dim(data$data)) <- names(dims) + # adding dates to Dates for the new NAs introduced + if ((data$Dates[[1]][2] - data$Dates[[1]][1]) == 1) { + timefreq <- 'days' + } else { + timefreq <- 'months' + warning("Time frequency of forecast time is considered monthly.") + } + start <- data$Dates[[1]] + dim(start) <- c(ftime = length(start)/dims['sdate'], sdate = dims['sdate']) + #new <- array(NA, prod(dim(data$data)[c('ftime', 'sdate')])) + # Pending fix transform to UTC when concatenaiting + data$Dates$start <- do.call(c, lapply(1:dim(start)[2], function(x) { + seq(start[1,x] - as.difftime(insert_ftime, + units = timefreq), + start[dim(start)[1],x], by = timefreq, tz = "UTC")})) + } + } if (is.null(indices)) { + if (any(split_dim %in% c('ftime', 'time', 'sdate'))) { if (is.list(data$Dates)) { indices <- data$Dates[[1]] } else { @@ -53,9 +90,10 @@ CST_SplitDim <- function(data, split_dim = 'time', indices = NULL, indices <- indices[1 : dim(data$data)[which(names(dim(data$data)) == split_dim)]] } + } } data$data <- SplitDim(data$data, split_dim = split_dim, indices = indices, - freq = freq) + freq = freq, new_dim_name = new_dim_name) return(data) } #'Function to Split Dimension @@ -67,10 +105,10 @@ CST_SplitDim <- function(data, split_dim = 'time', indices = NULL, #'@param data an n-dimensional array with named dimensions #'@param split_dim a character string indicating the name of the dimension to split #'@param indices a vector of numeric indices or dates -#'@param freq a character string indicating the frequency: by 'day', 'month' and 'year' or 'monthly' (by default). 'month' identifies months between 1 and 12 independetly of the year they belong to, while 'monthly' differenciates months from different years. Parameter 'freq' can also be numeric indicating the length in which to subset the dimension -#' +#'@param freq a character string indicating the frequency: by 'day', 'month' and 'year' or 'monthly' (by default). 'month' identifies months between 1 and 12 independetly of the year they belong to, while 'monthly' differenciates months from different years. Parameter 'freq' can also be numeric indicating the length in which to subset the dimension. +#'@param new_dim_name a character string indicating the name of the new dimension. #'@import abind -#'@import s2dverification +#'@importFrom s2dverification Subset #'@examples #' #'data <- 1 : 20 @@ -85,7 +123,8 @@ CST_SplitDim <- function(data, split_dim = 'time', indices = NULL, #'new_data <- SplitDim(data, indices = time, freq = 'month') #'new_data <- SplitDim(data, indices = time, freq = 'year') #'@export -SplitDim <- function(data, split_dim = 'time', indices, freq = 'monthly') { +SplitDim <- function(data, split_dim = 'time', indices, freq = 'monthly', + new_dim_name = NULL) { # check data if (is.null(data)) { stop("Parameter 'data' cannot be NULL.") @@ -123,6 +162,7 @@ SplitDim <- function(data, split_dim = 'time', indices, freq = 'monthly') { } indices <- rep(1 : (dims[pos_split] / freq), freq) indices <- sort(indices) + repited <- sort(unique(indices)) } } else if (is.numeric(indices)) { if (!is.null(freq)) { @@ -131,6 +171,7 @@ SplitDim <- function(data, split_dim = 'time', indices, freq = 'monthly') { "parameter 'indices' is numeric.") } } + repited <- sort(unique(indices)) } else { # Indices should be Dates and freq character if (!is.character(freq)) { @@ -161,19 +202,33 @@ SplitDim <- function(data, split_dim = 'time', indices, freq = 'monthly') { if (!is.numeric(indices)) { if (freq == 'day') { indices <- as.numeric(strftime(indices, format = "%d")) + repited <- unique(indices) } else if (freq == 'month') { indices <- as.numeric(strftime(indices, format = "%m")) + repited <- unique(indices) } else if (freq == 'year') { indices <- as.numeric(strftime(indices, format = "%Y")) + repited <- unique(indices) } else if (freq == 'monthly' ) { indices <- as.numeric(strftime(indices, format = "%m%Y")) + repited <- unique(indices) } else { stop("Parameter 'freq' must be numeric or a character: ", "by 'day', 'month', 'year' or 'monthly' (for ", "distinguishable month).") } } - repited <- unique(indices) + # check new_dim_name + if (!is.null(new_dim_name)) { + if (!is.character(new_dim_name)) { + stop("Parameter 'new_dim_name' must be character string") + } + if (length(new_dim_name) > 1) { + new_dim_name <- new_dim_name[1] + warning("Parameter 'new_dim_name' has length greater than 1 ", + "and only the first elemenst is used.") + } + } max_times <- max(unlist(lapply(repited, function(x){sum(indices == x)}))) data <- lapply(repited, function(x) {rebuild(x, data, along = split_dim, @@ -184,6 +239,9 @@ SplitDim <- function(data, split_dim = 'time', indices, freq = 'monthly') { } else { names(dim(data)) <- c(names(dims), 'index') } + if (!is.null(new_dim_name)) { + names(dim(data)) <- c(names(dims), new_dim_name) + } return(data) } @@ -200,3 +258,4 @@ rebuild <- function(x, data, along, indices, max_times) { } return(a) } + diff --git a/man/AdamontQQCorr.Rd b/man/AdamontQQCorr.Rd index 2006a7ad..ec49bad3 100644 --- a/man/AdamontQQCorr.Rd +++ b/man/AdamontQQCorr.Rd @@ -30,6 +30,9 @@ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and .NearestNeighbors \code{NN} output should be provided} +\item{wt_obs}{corresponding weather types (same dimensions as \code{obs} but +lat/lon)} + \item{corrdims}{list of dimensions in \code{exp} for which quantile mapping correction is applied} diff --git a/man/CST_AdamontQQCorr.Rd b/man/CST_AdamontQQCorr.Rd index fc4b1932..967ce171 100644 --- a/man/CST_AdamontQQCorr.Rd +++ b/man/CST_AdamontQQCorr.Rd @@ -26,6 +26,9 @@ can differ from \code{exp} if non rectilinear latlon grids are used, in which case regrid should be set to TRUE and .NearestNeighbors \code{NN} output should be provided} +\item{wt_obs}{corresponding weather types (same dimensions as \code{obs} but +lat/lon)} + \item{corrdims}{list of dimensions in \code{exp} for which quantile mapping correction is applied} @@ -34,12 +37,6 @@ correction is applied} \item{latdim}{character name of latitude dimension in \code{exp} and \code{obs}} - -\item{regrid}{(optional) boolean indicating whether .NearestNeighbors -regridding is needed} - -\item{NN}{list (optional, if regrid=TRUE) output from .NearestNeighbors -mapping (nlat, nlon) onto (nlat_o, nlon_o)} } \value{ an object of class \code{s2dv_cube} containing experiment data on the diff --git a/man/CST_SplitDim.Rd b/man/CST_SplitDim.Rd index ee93aedc..80a94da3 100644 --- a/man/CST_SplitDim.Rd +++ b/man/CST_SplitDim.Rd @@ -4,7 +4,14 @@ \alias{CST_SplitDim} \title{Function to Split Dimension} \usage{ -CST_SplitDim(data, split_dim = "time", indices = NULL, freq = "monthly") +CST_SplitDim( + data, + split_dim = "time", + indices = NULL, + freq = "monthly", + new_dim_name = NULL, + insert_ftime = NULL +) } \arguments{ \item{data}{a 's2dv_cube' object} @@ -14,10 +21,17 @@ CST_SplitDim(data, split_dim = "time", indices = NULL, freq = "monthly") \item{indices}{a vector of numeric indices or dates. If left at NULL, the dates provided in the s2dv_cube object (element Dates) will be used.} \item{freq}{a character string indicating the frequency: by 'day', 'month' and 'year' or 'monthly' (by default). 'month' identifies months between 1 and 12 independently of the year they belong to, while 'monthly' differenciates months from different years.} + +\item{new_dim_name}{a character string indicating the name of the new dimension.} + +\item{insert_ftime}{an integer indicating the number of time steps to add at the begining of the time series.} } \description{ This function split a dimension in two. The user can select the dimension to split and provide indices indicating how to split that dimension or dates and the frequency expected (monthly or by day, month and year). The user can also provide a numeric frequency indicating the length of each division. } +\details{ +Parameter 'insert_ftime' has been included for the case of using daily data, requiring split the temporal dimensions by months (or similar) and the first lead time doesn't correspondt to the 1st day of the month. In this case, the insert_ftime could be used, to get a final output correctly organized. E.g.: leadtime 1 is the 2nd of November and the input time series extend to the 31st of December. When requiring split by month with \code{inset_ftime = 1}, the 'monthly' dimension of length two will indicate the month (position 1 for November and position 2 for December), dimension 'time' will be length 31. For November, the position 1 and 31 will be NAs, while from positon 2 to 30 will be filled with the data provided. This allows to select correctly days trhough time dimension. +} \examples{ data <- 1 : 20 diff --git a/man/SplitDim.Rd b/man/SplitDim.Rd index f07e4756..a4904306 100644 --- a/man/SplitDim.Rd +++ b/man/SplitDim.Rd @@ -4,7 +4,13 @@ \alias{SplitDim} \title{Function to Split Dimension} \usage{ -SplitDim(data, split_dim = "time", indices, freq = "monthly") +SplitDim( + data, + split_dim = "time", + indices, + freq = "monthly", + new_dim_name = NULL +) } \arguments{ \item{data}{an n-dimensional array with named dimensions} @@ -13,7 +19,9 @@ SplitDim(data, split_dim = "time", indices, freq = "monthly") \item{indices}{a vector of numeric indices or dates} -\item{freq}{a character string indicating the frequency: by 'day', 'month' and 'year' or 'monthly' (by default). 'month' identifies months between 1 and 12 independetly of the year they belong to, while 'monthly' differenciates months from different years. Parameter 'freq' can also be numeric indicating the length in which to subset the dimension} +\item{freq}{a character string indicating the frequency: by 'day', 'month' and 'year' or 'monthly' (by default). 'month' identifies months between 1 and 12 independetly of the year they belong to, while 'monthly' differenciates months from different years. Parameter 'freq' can also be numeric indicating the length in which to subset the dimension.} + +\item{new_dim_name}{a character string indicating the name of the new dimension.} } \description{ This function split a dimension in two. The user can select the dimension to split and provide indices indicating how to split that dimension or dates and the frequency expected (monthly or by day, month and year). The user can also provide a numeric frequency indicating the length of each division. -- GitLab From 7437fd5666b263868521c4454c4aaac7f76c1599 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 22 Feb 2021 09:42:48 +0100 Subject: [PATCH 24/24] Fix dependencies --- NAMESPACE | 1 - R/CST_AdamontAnalog.R | 4 ++-- R/CST_AdamontQQCorr.R | 4 ++-- 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 659a4821..efc5c835 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -42,7 +42,6 @@ export(SplitDim) export(WeatherRegime) export(as.s2dv_cube) export(s2dv_cube) -import(ClimProjDiags) import(abind) import(ggplot2) import(multiApply) diff --git a/R/CST_AdamontAnalog.R b/R/CST_AdamontAnalog.R index 01e62d97..104913ac 100644 --- a/R/CST_AdamontAnalog.R +++ b/R/CST_AdamontAnalog.R @@ -32,7 +32,7 @@ #'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs #'analog values for each value of \code{exp} input data #'@import multiApply -#'@import ClimProjDiags +#'@importFrom s2dverification Subset #'@examples #'\dontrun{ #'wt_exp <- sample(1:3, 15*6*3, replace=T) @@ -131,7 +131,7 @@ CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, #'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs #'analog values for each value of \code{exp} input data #'@import multiApply -#'@import ClimProjDiags +#'@importFrom s2dverification Subset #'@examples #'\dontrun{ #'wt_exp <- sample(1:3, 15*6*3, replace=T) diff --git a/R/CST_AdamontQQCorr.R b/R/CST_AdamontQQCorr.R index be1dbbc8..3f9d7920 100644 --- a/R/CST_AdamontQQCorr.R +++ b/R/CST_AdamontQQCorr.R @@ -29,7 +29,7 @@ #' depending on the weather types \code{wt_exp} #' #'@import qmap -#'@import ClimProjDiags +#'@importFrom s2dverification Subset #'@import multiApply #'@import abind #'@examples @@ -136,7 +136,7 @@ CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, #'the weather types \code{wt_exp} #' #'@import qmap -#'@import ClimProjDiags +#'@importFrom s2dverification Subset #'@import multiApply #'@import abind #'@examples -- GitLab