QM based on doy.

R/03a_QM_DOY_utils.R

+# A list of functions used to compute the quantile mapping on a 91-day moving window, 
+# using .01 to .99 quantiles for deriving the correction 
+# applying a frequency adjustment for precipitation
+# functions derived from https://github.com/SvenKotlarski/qmCH2018 with small changes
+
+qm.cdf.doy <- function(series, yearbegin, yearend, doywindow, minq, maxq, incq, var, pr.wet){
+  
+  d		<- qm.doywindowgen(doywindow)
+  cdf_doy		<- array(NA,dim=c(365,length(seq(minq,maxq,by=incq))))
+  doy.series	<- qm.doystring(yearbegin,yearend)
+  fwet_doy	<- rep(NA,365)
+  cdfwet_doy	<- array(NA,dim=c(365,101))
+  
+  # *** set 29th February to 1st March ***
+  doy.series[doy.series > 500] <- 60
+  
+  # *** check of input data: length OK? ***
+  if(length(series) != length(doy.series)){
+    print("ERROR: data series and doy series not equally long")
+    cdf_doy <- NULL
+  }else{
+    for(i in 1:365){
+      #			*** match: returns vector of length(doy.series): Position of matches in d$doy.matrix[i,] or NA if no match
+      #			*** => all positions in doy.series that correspond to the window defined in d$doy.matrix[i,] are marked
+      cdf_doy[i,] <- quantile(series[which(is.na(match(doy.series,d$doy.matrix[i,])) == FALSE)], probs=seq(minq,maxq,by=incq),na.rm=TRUE)
+      
+      #			*** possible (shorter) alternative ***
+      #			cdf_doy[i,] <- quantile(series[which(!is.na(match(doy.series,d$doy.matrix[i,])))],probs=c(seq(minq,maxq,by=incq)),na.rm=TRUE)
+      
+      #			*** for precipitation additionally compute wet day frequency and wet day distribution for each DOY ***
+      if (var == 'pr') {
+        # *** all values for current DOY including window ***
+        series.doywin <- series[which(is.na(match(doy.series,d$doy.matrix[i,])) == FALSE)]
+        # *** total number of days w/o NA ***
+        ndays <- length(series.doywin[!is.na(series.doywin)])
+        # *** number of wet days w/o NA ***
+        ndays.wet <- length(which(series.doywin >= pr.wet))
+        # *** wet day fraction for current DOY ***
+        fwet_doy[i] <- ndays.wet / ndays
+        # *** all values for current DOY including window but excluding dry days ***
+        series.doywin.wet <- series.doywin[series.doywin >= pr.wet]
+        # *** wet day CDF for current DOY including window ***
+        cdfwet_doy[i,] <- quantile(series.doywin.wet,probs=seq(0,1,0.01),na.rm=TRUE)
+      }
+    }
+  }
+  return(list(cdf.matrix = cdf_doy, f.wet = fwet_doy, cdf.wet = cdfwet_doy))
+}
+
+
+qm.empqm.doy.fa <- function(series, doy.series, cdf.mod, cdf.obs, cdf.obs.wet, cor.fun, fwet.obs, fwet.mod, pr.wet, method=c("linear","bin"), var, minq, maxq, incq){
+
+  if(length(series) != length(doy.series)){ # *** a little check
+    print(paste("ERROR: length of time series (",length(series),") and DOY series (",length(doy.series),") do not match.", sep=""))
+  }else{
+    series_q		<- rep(NA,length(series))
+    series_c		<- rep(NA,length(series))
+    for(i in 1:length(series)){
+      doy <- doy.series[i]
+      if(doy == 999){ doy <- 60 }
+      if(is.na(series[i]) == "TRUE"){
+        series_q[i] <- NA
+        series_c[i] <- NA
+      }else{
+        # *** OPTION A: if dry day in model and if modelled wet day frequency <= observed wet day frequency apply frenquency adjustment ***
+        if ((series[i]==0) & (fwet.mod[doy] <= fwet.obs[doy])) {
+          # *** generate random number ***
+          rnum <- runif(1)
+          # *** determine fraction of dry days in model that needs to become wet in order to represent observed wet day frequency ***
+          probwet <- 1 - ((1-fwet.obs[doy])/(1-fwet.mod[doy]))
+          # *** if rnum <= probwet: fill with a precipitation day; else assign original value, i.e. zero precip ***
+          if (rnum <= probwet) {
+            # *** randomly draw from observed wet day distribution and assign that value ***
+            fun <- approxfun(seq(0,1,0.01),cdf.obs.wet[doy,])
+            random.number <- runif(1)
+            series_c[i] <- fun(random.number)
+          }else{
+            series_c[i] <- series[i]
+          }
+          # *** in any case, assign NA to series of modelled quantile numbers: means FA has been carried out! ***
+          series_q[i]	<- NA
+        }else{
+          # *** OPTION B: else, don't apply frenquency adjustment; apply regular correction function ***
+          # *** using quantile-"bins" ***
+          if(match.arg(method) == "bin"){
+            # *** Note: this command implicitly selects the first (last) quantile considered for new extremes , i.e.  ***
+            # *** for smaller (larger) values than those occuring in the calbration period (in case of zeros: first   ***
+            # *** occurence). If minq=0.01 and maxq=0.99 nex extremes are hence corrected according to the correction ***
+            # *** function for the 1st and 99th percentile.                                                           ***
+            ix		<- which.min(abs(cdf.mod[doy,]-series[i]))
+            # *** Compute quantile index  ***
+            series_q[i]	<- minq + (ix-1)*incq
+            # *** Compute corrected value ***
+            series_c[i]	<- series[i]-cor.fun[doy,ix]
+          }
+          
+          # *** using linear interpolation between quantiles ***
+          if(match.arg(method) == "linear"){
+            # *** number of interpolated sampling points ***
+            npoints <- 1000
+            cor		<- approx(cdf.mod[doy,],n=npoints,method="linear")$y
+            # *** Note: this command implicitly selects the first (last) quantile considered for new extremes , i.e.  ***
+            # *** for smaller (larger) values than those occuring in the calbration period (in case of zeros: first   ***
+            # *** occurence). If minq=0.01 and maxq=0.99 nex extremes are hence corrected according to the correction ***
+            # *** function for the 1st and 99th percentile.                                                           ***
+            ix		<- which.min(abs(cor-series[i]))
+            # *** Compute quantile index ***
+            series_q[i]	<- minq + (ix-1) * ((maxq-minq)/(npoints-1))
+            # *** Compute correction function ***
+            cf		<- approx(cor.fun[doy,],n=npoints,method="linear")$y
+            # *** Compute corrected value ***
+            series_c[i]	<- series[i]-cf[ix]
+          }
+        }
+        
+      }
+    }
+    
+    # *** for huss (specific humidity), sfcWind (wind speed), hurs (relative humidity), rsds (global radiation): set negative values to 0 ***
+    #if ((var == 'huss') || (var == 'sfcWind') || (var == 'hurs') || (var == 'rsds')) series_c[which(series_c<0)] <- 0
+    
+    return(list(qm.input.series = series, qm.corrected.series = series_c, quantile.index  = series_q))
+  }
+}
+
+
+
+qm.doywindowgen	<- function(window.width){
+  
+  doys 		<- c(1:365)
+  ndays 		<- length(doys)
+  half.width	<- (window.width/2) - 0.5
+  
+  doy_ix <- array(data=NA,dim=c(ndays,window.width))
+  temp.vec <- rep(1:ndays,3)
+  
+  for (day in doys) { doy_ix[day,] <- temp.vec[((ndays+day)-half.width):((ndays+day)+half.width)] }
+  
+  return(list(half.width = half.width,window.width = window.width, doy.matrix = doy_ix))
+}
+
+qm.doystring <- function(yearbegin,yearend){
+  
+  doyreg 		<- c(1:365)
+  doyleap		<- c(1:59,999,60:365)
+  doyscenario     <- integer()
+  
+  for (countyear in yearbegin:yearend) {
+    if (leap_year(countyear)) doyscenario <- c(doyscenario,doyleap) else doyscenario <- c(doyscenario,doyreg)
+  }
+  return(doyscenario)
+}

R/03b_QM_DOY.R

+library(raster)
+library(rasterVis)
+library(rgdal)
+library(ncdf4)
+library(ncdf.tools)
+library(ncdf4.helpers)
+library(data.table)
+library(magrittr)
+#library(doParallel)
+#library(foreach)
+library(PCICt)
+
+# the code performes the QM correction on daily data by using a 91-day moving window
+# the correction is performed for quantiles from 0.01 to 0.99 in order to better account for extremes
+# a frequency adjustment is also performed for precipitation only
+# NOTE: quite slow - optimization and parallelization is needed
+
+#registerDoParallel(parallel::detectCores() - 1)
+
+path_in_mod<-"/../CORDEX/" #already remapped and cropped files 
+variable<-"tasmin"
+inventory<-readRDS(paste0(path_in_mod,"/",variable,"_inventory.RDS")) #or a list for the selected models
+CORDEX.models<-unique(paste0(inventory$gcm,"+",inventory$institute_rcm)) # total number of used models
+hist_years <- 1950:2005 # years in historical period
+proj_years <- 2006:2100 # years in projections
+cal_years  <- 1981:2010 # years for the calibration
+cor_years  <- 1950:2100 # years to correct
+wday.width<-91 # window centred on the DOY to correct
+minq<-0.01; maxq<-0.99; incq<-0.01 # definition of quantiles min, max and step
+FA <- TRUE # frequency adjustment for precipitation?
+source("/../03a_QM_DOY_utils.R") # load the functions required for QM
+# step 1
+# load daily observations for the variable under evaluation and for the calibration period 
+# observations are split by month and stored in separate folders for each year
+
+path_in_obs<-"/../OBS_DATA/"
+list<-list.files(paste0(path_in_obs,variable,"/",cal.years[1]),pattern=".nc",full.names = T)
+obs.stack<-stack(lapply(list,stack))
+for(i in cal.years[-1]){
+  list<-list.files(paste0(path_in_obs,variable,"/",i),pattern=".nc",full.names = T)
+  suppressWarnings(tmp<-stack(lapply(list,stack)))
+  obs.stack<-stack(obs.stack,tmp)
+}
+
+idx.TSTcell<-which(!is.na(getValues(obs.stack[[1]]))) # indices of cells which will be processed
+scenario<-"rcp85"
+# convert stack into matrix with only points of interest 
+# row=time col=points
+obs.dat<-as.data.frame(obs.stack)
+obs.dat<-obs.dat[idx.TSTcell,]; obs.dat<-t(data.matrix(obs.dat))
+
+# step 2 
+# load CORDEX for the specific simulation and run the QM correction
+dates.ref<-seq(as.Date(paste0(min(cal_years),"/01/01")),as.Date(paste0(max(cal_years),"/12/31")),"day")
+dates.hist.full<-seq(as.Date(paste0(min(hist_years),"/01/01")),as.Date(paste0(max(hist_years),"/12/31"),"day")) # CORDEX historical period
+dates.proj.full<-seq(as.Date(paste0(min(proj_years),"/01/01")),as.Date(paste0(max(proj_years),"/01/01")),"day") # CORDEX scenario period
+dates.full<-c(dates.hist.full,dates.proj.full) # full analysed period
+
+# loop over each gcm+rcm simulation --> this loop should be parallelized
+for(i.mod in seq_along(CORDEX.models)){ 
+  
+  i.gcm<-strsplit(CORDEX.models[i.mod],"+",fixed=T)[1]
+  i.rcm<-strsplit(CORDEX.models[i.mod],"+",fixed=T)[2]
+  
+  #load full data series for the gcm+rcm simulation
+  dat.mod<-load_CORDEX(gcm=i.gcm,rcm=i.rcm,
+                       path_in_mod=path_in_mod, variable=variable,
+                       cells=idx.TSTcell,
+                       dates.historical=dates.hist.full,
+                       dates.proections=dates.proj.full)
+  
+  # matrix for storing QM corrected simulations
+  # row=time, col=points
+  qm.mat<-matrix(ncol=length(idx.TSTcell),nrow=length(dates.full))
+  
+  for(i.cell in 1:length(idx.TSTcell)){
+    # load the series for the test point
+    series.mod.cal<-dat.mod[which(year(dates.full)%in%cal_years),i.cell]
+    series.obs.cal<-as.numeric(obs.dat.conv[which(year(dates.ref)%in%cal_years),i.cell])
+    series.mod.all<-dat.mod[which(year(dates.full)%in%cor_years),i.cell]
+    
+    # set to zero pr correspomding to dry days
+    if ((variiable=='pr') & FA) {
+      series.mod.cal[series.mod.cal < pr.wet] <- 0
+      series.obs.cal[series.obs.cal < pr.wet] <- 0
+      series.mod.all[series.mod.all < pr.wet] <- 0
+    }
+    
+    # extract the DOY index for the full model series 
+    # the index for 29th Feb is set to 999 in all leap years
+    # 29th Feb is then handled as 1st Mar
+    doy.full <- qm.doystring(min(cor_years), max(cor_years))
+    
+    # the CDF for each DOY for both model and observation data is computed
+    # the CDF for wet days and the wet day frequency for each DOY are also reported. They will be used for the frequency adjustment
+    c_mod <- qm.cdf.doy(series.mod.cal, min(cal_years), max(cal_years), wday.width, minq, maxq, incq, variable, pr.wet)
+    c_obs <- qm.cdf.doy(series.obs.cal, min(cal_years), max(cal_years), wday.width, minq, maxq, incq, variable, pr.wet)
+    
+    
+    # here the correction for each quantile and DOY (additive) is computed
+    if((length(c_mod$cdf.matrix[1,]) != length(c_obs$cdf.matrix[1,])) & (length(c_mod$cdf.matrix[,1]) != length(c_obs$cdf.matrix[,1]))){
+      print("ERROR: doy-based CDF of model and observations do not have the same dimension")
+    }else{
+      dif.cdf	<- c_mod$cdf.matrix-c_obs$cdf.matrix
+    }
+    
+    if(variable == "pr"){
+      # apply the quantile correction by adding an adjustment of precipitation frequency in model data
+      cc <- qm.empqm.doy.fa(series.mod.all, doy.full, c_mod$cdf.matrix, c_obs$cdf.matrix, c_obs$cdf.wet,
+                            dif.cdf, c_obs$f.wet, c_mod$f.wet, pr.wet, method=qm.method, variable, minq, maxq, incq)
+      # check if negative corrected precipitation occur and set to 0
+      cc$qm.corrected.series[cc$qm.corrected.series < pr.wet] <- 0
+    } else {
+      cc <- qm.empqm.doy(series.mod.all.in, doy.full, c_mod$cdf.matrix, dif.cdf, 
+                         method=qm.method, variable, minq, maxq, incq)
+    }
+    
+    # store the corrected series in the output matrix
+    qm.mat[,i.cell]<-cc$qm.corrected.series
+  } # close loop on grid cell
+} # close loop on models
\ No newline at end of file