Commit 99efbbbb authored by Franssen, Wietse's avatar Franssen, Wietse
Browse files

Added new biascorrection made by Wouter for precipitation

parent bc1079d9
rm(list=ls())
library(ecomsUDG.Raccess)
library(downscaleR)
library(WFRTools)
library(ggplot2)
source(file = "./infoGeneral.R")
source(file = "./functionsGeneral.R")
source(file = "./functions/infoGeneral.R")
source(file = "./functions/functionsGeneral.R")
submitscript <- FALSE
if (submitscript) {
members <- c(1:15)
initYears <-c(1981:2010)
currMonths <- c(X:X)
targetMonths <- c(X:X)
leadMonths<-c(0:6)
locName <- 'X'
inPathNoBC<-"./testData/reFormatted_noBC"
inPathBC<-"./testData/reFormatted_BC"
outPath<-"./testData/pics/checkBias"
resolution<-"0.50"
inPathNoBC<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_noBC", resolution, locName)
inPathBC<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_BC", resolution, locName)
outPath<-"../Analysis_output/pics/checkBias_rev1.1_otherVars"
} else {
members <-c(1:15)
initYears <-c(1981:2010)
currMonths<-1
targetMonths<-c(1:12)
leadMonths<-c(0:6)
#locName<-"GHA"
locName<-"EU"
inPathNoBC<-"./testData/reFormatted_noBC"
inPathBC<-"./testData/reFormatted_BC"
outPath<-"./testData/pics/checkBias"
resolution<-"0.50"
inPathNoBC<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_noBC", resolution, locName)
inPathBC<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_BC", resolution, locName)
outPath<-"../Analysis_output/pics/checkBias_rev1.1_otherVars"
}
variables<-names(variableInfo)
variables<-c( "pr" , "psl" , "rsds" , "rlds" , "tasmin","tasmax", "huss" , "sfcWind")
variables<-c("tas")
variables<-c("sfcWind")
#variables<-c("pr")
dir.create(outPath, recursive = TRUE, showWarnings = FALSE)
print("start")
for (variableName in variables) {
for (currMonth in currMonths) {
for (targetMonth in targetMonths) {
for (leadMonth in leadMonths) {
targetSYear <- getInitTargetInfo( initYear = initYears[1], targetMonth = targetMonth, leadMonth = leadMonth )$targetYear
targetEYear <- getInitTargetInfo( initYear = initYears[length(initYears)], targetMonth = targetMonth, leadMonth = leadMonth )$targetYear
######## OPEN FILES #####################
noBCFile <- sprintf("%s/%s_forcing_seas15_%s_noBC_E%02d-%02d_%4d-%4d_%02d_LM%d.RData",
noBCFile <- sprintf("%s/%s_forcing_seas15_%s_noBC_E%02d-%02d_TAR%4d-%4d_%02d_LM%d.RData",
inPathNoBC, variableName, locName,
members[1], members[length(members)],
initYears[1], initYears[length(initYears)],
currMonth, leadMonth)
targetSYear, targetEYear,
targetMonth, leadMonth)
print(sprintf("Open: %s", noBCFile))
load(noBCFile)
RData_noBC<-RData
rm(RData)
BCFile <- sprintf("%s/%s_forcing_seas15_%s_BC_E%02d-%02d_%4d-%4d_%02d_LM%d.RData",
BCFile <- sprintf("%s/%s_forcing_seas15_%s_BC_E%02d-%02d_TAR%4d-%4d_%02d_LM%d.RData",
inPathBC, variableName, locName,
members[1], members[length(members)],
initYears[1], initYears[length(initYears)],
currMonth, leadMonth)
targetSYear, targetEYear,
targetMonth, leadMonth)
print(sprintf("Open: %s", BCFile))
load(BCFile)
RData_BC<-RData
rm(RData)
obsFile <- sprintf("/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/DATA_BIAS/DATA/DATA_RAW_R/%s_WFD/%s_wfdei_%4d-%4d_%02d.RData",
if (variableName == "pr" ) {
obsFile <- sprintf("/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/biascorrection/DATA/wfdei_rev3.0/0.50deg/%s/%s_wfdei_GPCC_EU_%4d-%4d_%02d.RData",
locName, variableName,
initYears[1], initYears[length(initYears)],
currMonth)
targetMonth)
} else {
obsFile <- sprintf("/home/WUR/frans004/L_BACKUP/PROJECTS/EUPORIAS/biascorrection/DATA/wfdei_rev3.0/0.50deg/%s/%s_wfdei_EU_%4d-%4d_%02d.RData",
locName, variableName,
initYears[1], initYears[length(initYears)],
targetMonth)
}
print(sprintf("Open: %s", obsFile))
load(obsFile)
obs_new<-obs
newGrid<-getGrid(RData_noBC)
indexesLat<-c(which(obs$xyCoords$y==newGrid$y[1]):which(obs$xyCoords$y==newGrid$y[2]))
indexesLon<-c(which(obs$xyCoords$x==newGrid$x[1]):which(obs$xyCoords$x==newGrid$x[2]))
obs_new$xyCoords<-RData_noBC$xyCoords
obs_new$Data<-obs$Data[,indexesLat,indexesLon]
attr(obs_new$Data,"dimensions") <- c("time","lat","lon")
obs<-obs_new
rm(obs_new)
obs<-RData
rm(RData)
#############################
......@@ -87,23 +89,24 @@ for (variableName in variables) {
barLimitsDiff<-c(colorbarLimits[[variableName]]$diffMin,colorbarLimits[[variableName]]$diffMax)
print("obs_mean...")
obs_mean<-obs;obs_mean$Data<-apply(obs$Data, c(2,3), mean)
p_obs<-ggplotje(obs_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI (", month.name[currMonth], ")"))
# p_obs<-ggplotje(obs_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI (", month.name[currMonth], ")"),barLimits = barLimits)
p_obs<-ggplotje(obs_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI (", month.name[targetMonth], ")"),barLimits = barLimits)
# p_obs<-ggplotje(obs_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI (", month.name[targetMonth], ")"))
print("BC_mean...")
RData_BC_mean<-RData_BC;RData_BC_mean$Data<-apply(RData_BC$Data, c(1,3,4), mean)
p_BC<-ggplotje(RData_BC_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nseas15_BC (", month.name[currMonth], ", leadMonth:", leadMonth, ")"),barLimits = barLimits)
# p_BC<-ggplotje(RData_BC_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nseas15_BC (", month.name[currMonth], ", LM:", leadMonth, "/SM:", iStartMonth, "/S:", iSeason, ")"))
p_BC<-ggplotje(RData_BC_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nseas15_BC (", month.name[targetMonth], ", leadMonth:", leadMonth, ")"),barLimits = barLimits)
# p_BC<-ggplotje(RData_BC_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nseas15_BC (", month.name[targetMonth], ", LM:", leadMonth, "/SM:", iStartMonth, "/S:", iSeason, ")"))
print("noBC_mean...")
RData_noBC_mean<-RData_noBC;RData_noBC_mean$Data<-apply(RData_noBC$Data, c(1,3,4), mean)
p_noBC<-ggplotje(RData_noBC_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nseas15_noBC (", month.name[currMonth], ", leadMonth:", leadMonth, ")"),barLimits = barLimits)
# p_noBC<-ggplotje(RData_noBC_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nseas15_noBC (", month.name[currMonth], ", LM:", leadMonth, "/SM:", iStartMonth, "/S:", iSeason, ")"))
p_noBC<-ggplotje(RData_noBC_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nseas15_noBC (", month.name[targetMonth], ", leadMonth:", leadMonth, ")"),barLimits = barLimits)
# p_noBC<-ggplotje(RData_noBC_mean, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nseas15_noBC (", month.name[targetMonth], ", LM:", leadMonth, "/SM:", iStartMonth, "/S:", iSeason, ")"))
RData_BC_mean2<-RData_BC
RData_BC_mean2$Data<-apply(RData_BC_mean$Data, c(2,3), mean)
obsdiff<-obs_mean
obsdiff$Data<-obs_mean$Data-RData_BC_mean2$Data
# p_obsVsBC<-ggplotje(obsdiff, title = "Observed vs BC")
p_obsVsBC<-ggplotje(obsdiff, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI - seas15_BC"),barLimits = barLimitsDiff,barDiff=TRUE,backgroudColor = "grey85")
# p_obsVsBC<-ggplotje(obsdiff, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI - seas15_BC"),barLimits = barLimitsDiff,barDiff=TRUE,backgroudColor = "grey85")
p_obsVsBC<-ggplotje(obsdiff, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI - seas15_BC"),backgroudColor = "grey85")
# p_obsVsBC<-ggplotje(obsdiff, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI - seas15_BC"),barDiff=TRUE,backgroudColor = "grey85")
RData_noBC_mean2<-RData_BC
......@@ -111,7 +114,8 @@ for (variableName in variables) {
obsdiffnoBC<-obs_mean
obsdiffnoBC$Data<-obs_mean$Data-RData_noBC_mean2$Data
# p_obsVsnoBC<-ggplotje(obsdiffnoBC, title = "Observed vs noBC")
p_obsVsnoBC<-ggplotje(obsdiffnoBC, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI - seas15_noBC"),barLimits = barLimitsDiff,barDiff=TRUE,backgroudColor = "grey85")
# p_obsVsnoBC<-ggplotje(obsdiffnoBC, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI - seas15_noBC"),barLimits = barLimitsDiff,barDiff=TRUE,backgroudColor = "grey85")
p_obsVsnoBC<-ggplotje(obsdiffnoBC, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI - seas15_noBC"),backgroudColor = "grey85")
# p_obsVsnoBC<-ggplotje(obsdiffnoBC, title = paste0(variableInfo[[variableName]]$longName, " ", variableInfo[[variableName]]$units, "\nWFDEI - seas15_noBC"),barDiff=TRUE,backgroudColor = "grey85")
......@@ -157,11 +161,11 @@ for (variableName in variables) {
p6<-densPlot(10.25,46.25)
}
nameFilePicture <- sprintf("%s/%s_%s_E%02d-%02d_%4d-%4d_%02d_LM%d",
nameFilePicture <- sprintf("%s/%s_%s_E%02d-%02d_TAR%4d-%4d_%02d_LM%d",
outPath, variableName, locName,
members[1], members[length(members)],
initYears[1], initYears[length(initYears)],
currMonth, leadMonth)
targetMonth, leadMonth)
print(paste0("Saving picture : ", nameFilePicture))
......
rm(list=ls())
library(ecomsUDG.Raccess)
library(downscaleR)
library(MASS)
source(file = "./functions/functionsGeneral.R")
source(file = "./functions/infoGeneral.R")
source(file = "./functions/functionMakeSelection.R")
source(file = "./functions/functionBiascorrectionWUR.R")
#source(file = "./functions/functionBiasCorrectionECOMS_precip_adapt.R")
#source(file = "./misc.R")
submitscript <- FALSE
if (submitscript) {
members <- c(1:15)
initYears <-c(1981:2010)
currMonths <- c(X:X)
targetMonths <- c(X:X)
leadMonths<-c(0:6)
locName <- 'X'
resolution<-"0.50"
obsPath<-sprintf("../DATA/wfdei_rev3.0/%sdeg/%s/", resolution, locName)
# inPath<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_noBC", resolution, locName)
inPath<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_noBC", resolution, locName)
outPath<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_BC", resolution, locName)
# outPath<-sprintf("../DATA/System4_seasonal_15_rev1.2_1.0mm_adapt/%sdeg/%s_BC", resolution, locName)
# outPath<-sprintf("../DATA/System4_seasonal_15_rev1.1_CRU_0.1mm_new/%sdeg/%s_BC", resolution, locName)
outPath<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg_targetYears/%s_BC", resolution, locName)
} else {
members <-c(1:15)
initYears <-c(1981:2010)
currMonths<-1
targetMonths<-6
leadMonths<-c(0:6)
#locName<-"GHA"
locName<-"EU"
resolution<-"0.50"
obsPath<-sprintf("../DATA/wfdei_rev3.0/%sdeg/%s/", resolution, locName)
inPath<-"./testData/reFormatted_noBC"
outPath<-"./testData/reFormatted_BC"
inPath<-sprintf("../DATA/System4_seasonal_15_rev1.1/%sdeg/%s_noBC", resolution, locName)
outPath<-sprintf("../DATA/System4_seasonal_15_rev1.2_1.0mm_adapt/%sdeg/%s_BC", resolution, locName)
}
variables<-names(variableInfo)
variables<-c( "pr", "sfcWind", "rsds", "rlds", "huss", "tas", "tasmin", "tasmax")
#variables<-c("tas")
variables<-c( "tas", "tasmin", "tasmax", "rsds", "rlds", "pr", "huss", "sfcWind")
variables<-c( "tas", "tasmin", "tasmax", "rsds", "rlds", "huss", "sfcWind")
variables<-c("pr")
dir.create(outPath, recursive = TRUE, showWarnings = FALSE)
print("start")
print("start GPCC")
for (variableName in variables) {
for (currMonth in currMonths) {
for (targetMonth in targetMonths) {
for (leadMonth in leadMonths) {
iFile <- sprintf("%s/%s_forcing_seas15_%s_noBC_E%02d-%02d_%4d-%4d_%02d_LM%d.RData",
targetSYear <- getInitTargetInfo( initYear = initYears[1], targetMonth = targetMonth, leadMonth = leadMonth )$targetYear
targetEYear <- getInitTargetInfo( initYear = initYears[length(initYears)], targetMonth = targetMonth, leadMonth = leadMonth )$targetYear
# if (targetSYear != initYears[1]) {
iFile <- sprintf("%s/%s_forcing_seas15_%s_noBC_E%02d-%02d_TAR%4d-%4d_%02d_LM%d.RData",
inPath, variableName, locName,
members[1], members[length(members)],
initYears[1], initYears[length(initYears)],
currMonth, leadMonth)
targetSYear, targetEYear,
targetMonth, leadMonth)
print(sprintf("Open: %s", iFile))
load(iFile)
RData_noBC<-RData
rm(RData)
obsFile <- sprintf("%s/%s_wfdei_%s_%4d-%4d_%02d.RData",
obsPath, variableName, locName,
initYears[1], initYears[length(initYears)],
currMonth)
if (variableName == "pr") {
obsFile <- sprintf("%s/%s_wfdei_GPCC_%s_%4d-%4d_%02d.RData",
obsPath, variableName, locName,
initYears[1], initYears[length(initYears)],
targetMonth)
# obsFile <- sprintf("%s/%s_wfdei_%s_%4d-%4d_%02d.RData",
# obsPath, variableName, locName,
# initYears[1], initYears[length(initYears)],
# targetMonth)
} else {
obsFile <- sprintf("%s/%s_wfdei_%s_%4d-%4d_%02d.RData",
obsPath, variableName, locName,
initYears[1], initYears[length(initYears)],
targetMonth)
}
print(sprintf("Open: %s", obsFile))
load(obsFile)
......@@ -67,18 +94,36 @@ for (variableName in variables) {
# RData<-RData_new
# rm(RData_new)
## Make uniform and match the selections
dates1<-as.character(as.POSIXct(RData_obs$Dates$start, tz = "GMT"))
dates2<-as.character(as.POSIXct(RData_noBC$Dates$start, tz = "GMT"))
overlapping_dates<-as.Date(intersect(dates1, dates2))
rm(dates2,dates1)
RData_obs <- subsetRData(RData_obs, dates = overlapping_dates)
RData_noBC_pred <- subsetRData(RData_noBC, dates = overlapping_dates)
print(sprintf("Bais correction..."))
RData <- biasCorrection (RData_obs, RData_noBC, RData_noBC, method = "qqmap", multi.member = FALSE, pr.threshold = 0.1)# bias
oFile <- sprintf("%s/%s_forcing_seas15_%s_BC_E%02d-%02d_%4d-%4d_%02d_LM%d.RData",
if (variableName == "pr") {
## Correction: PRECIP
print(sprintf("Bias correction..."))
RData <- biasCorrect(RData_obs, RData_noBC, nbin = 100, epsilon = 0.00001)
} else {
## Correction: other variables
print(sprintf("Bias correction..."))
RData <- biasCorrection (RData_obs, RData_noBC_pred, RData_noBC, method = "qqmap", multi.member = FALSE, pr.threshold = 0.1)# bias
# RData <- biasCorrection (RData_obs, RData_noBC, RData_noBC, method = "qqmap", multi.member = FALSE, pr.threshold = 0.1)# bias
}
## Write output
oFile <- sprintf("%s/%s_forcing_seas15_%s_BC_E%02d-%02d_TAR%4d-%4d_%02d_LM%d.RData",
outPath, variableName, locName,
members[1], members[length(members)],
initYears[1], initYears[length(initYears)],
currMonth, leadMonth)
targetSYear, targetEYear,
targetMonth, leadMonth)
print(sprintf("Saving: %s", oFile))
save(file= oFile, RData)
# }
}
}
}
rm(list=ls())
library(downscaleR)
library(WFRTools)
source(file = "./functions/infoGeneral.R")
source(file = "./functions/functionsGeneral.R")
initYears <-c(1981:2010)
currMonths<-c(1:12)
#locName<-"GHA"
locName<-"EU"
resolution<-"0.75"
inPath<-sprintf("../DATA/wfdei_rev3.0/%sdeg/%s", resolution, locName)
outPath<-sprintf("../DATA/wfdei_rev3.1/%sdeg/%s", resolution, locName)
variables<-c("pr")
dir.create(outPath, recursive = TRUE, showWarnings = FALSE)
print("start")
for (variableName in variables) {
for (currMonth in currMonths) {
inFile <- sprintf("%s/%s_wfdei_GPCC_%s_%4d-%4d_%02d.RData",
inPath, variableName, locName,
initYears[1], initYears[length(initYears)],
currMonth)
print(sprintf("Open: %s", inFile))
load(inFile)
RData$Data<-aperm(RData$Data, c(3,2,1))
attr(RData$Data,"dimensions") <- c("time","lat","lon")
outFile <- sprintf("%s/%s_wfdei_GPCC_%s_%4d-%4d_%02d.RData",
outPath, variableName, locName,
initYears[1], initYears[length(initYears)],
currMonth)
print(sprintf("Saving: %s", outFile))
save(file=outFile, RData)
}
}
This diff is collapsed.
This diff is collapsed.
biasCorrectCell<-function(obscell, predcell, minobs, npred, nbin = 100, epsilon = 0.00001){
# Due to set.seed, the same random numbers are generated each time
# the routine is run
set.seed (423)
rannumdis <- runif(npred, -0.5, 0.5)
set.seed (187)
rannum0 <- runif(npred, 0.0, 1.0)
# To avoid that non-zero values occur more than once
indab0 <- which (predcell > 0.0)
predcell[indab0] <- predcell[indab0] + epsilon * rannumdis[indab0]
# Precipitation cannot be less than 0
indbel0 <- which (predcell < 0.0)
predcell[indbel0] = 0.0
# The array with the corrected hindcasts gets the same structure as the
# hindcasts
simcell <- predcell
# Determine the indices of the valid observations and hindcasts
indvalobs <- which (!is.na (obscell) == T)
indvalpred <- which (!is.na (obscell) == T)
# Go to next cell if there are no valid observations and hindcasts
if (length (indvalobs) == 0) return(simcell)
if (length (indvalpred) == 0) return(simcell)
if (length (indvalobs) < minobs) stop ("Less observations than expected")
# Select valid observations and the corresponding predictions
obsval = obscell[indvalobs]
predval = predcell[indvalobs]
if (sum (is.na (predval)) > 0) stop ("NaN in hindcasts")
# Set bin properties, once per data set
nobs0 = length(obsval)
# Boundaries of the bins (counts)
# bounnsam can be a non-integer
# ibts and iets are indices of that delimit the bins containing the
# sorted observations and hindcasts
# nsambin is the number of samples per bin
bounnsam <- seq(0, nobs0, length.out = nbin + 1)
ibts <- floor(bounnsam[1:nbin] + 1)
iets <- floor(bounnsam[2:(nbin+1)])
nsambin <- iets - ibts + 1
meanobsbin <- vector ("double", nbin)
meanpredbin <- vector ("double", nbin)
thrpredbin <- vector ("double", nbin+1)
thrpredbin[1] = -999.0e8
thrpredbin[nbin+1] = +999.0e8
nobscell = length(obsval)
if (nobscell != nobs0) stop ("number of observations differs between cells")
# Sort the observations and the predictions
obssort <- sort (obsval)
predsort <- sort (predval)
# Compute bin mean values, the correction for each bin and
# the correction for each bin
for (ibin in 1:nbin) {
obsbin = obssort [ibts[ibin]:iets[ibin]]
meanobsbin[ibin] = mean(obsbin)
predbin = predsort [ibts[ibin]:iets[ibin]]
meanpredbin[ibin] = mean(predbin)
if (ibin != nbin) thrpredbin[ibin+1] = 0.5 * (predsort[iets[ibin]] +
predsort[iets[ibin] + 1])
}
diffbin = meanobsbin - meanpredbin
# Number of bins with no precipitation
nbins0obs <- length (which (meanobsbin == 0))
nbins0pred <- length (which (meanpredbin == 0))
morepred0 = F
# If there are more bins with zero precipitation in the hindcasts
# than in the observations (morepred0 = T), there is an issue,
# which solution is prepared here
# allcumsum0bins is the cumulative number of zeros
# over the bins containing zeros (in the hindcasts)
# thrattrrandom is the same as allcumsum0bins but rescaled to {0,1]
if (nbins0pred > nbins0obs) {
morepred0 = T
nsam0bins = nsambin[1:nbins0pred]
cumsumsam0bins <- cumsum (nsam0bins)
n0transbin = 0
if (nbins0pred != nbin) {
itrbin = nbins0pred + 1
n0transbin <- length (which (predsort [ibts[itrbin]:iets[itrbin]] == 0))
}
n0tot <- sum (nsam0bins) + n0transbin
allcumsum0bins <- c (cumsumsam0bins, n0tot)
thrattrrandom <- 1.0 * allcumsum0bins / n0tot
}
irn = 0
# Correct simulated values (simcell)
for (isam in (1:npred)) {
# Take the uncorrected value
predsam = predcell[isam]
# Skip NaN
if ( is.na(predsam)) {
simcell[isam] = NaN
next
}
# Determine bin number for the special cases (morepred0 = T)
# and simulated zero values
if (morepred0 & predsam == 0) {
irn = irn + 1
rnhere = rannum0[irn]
diffthr = thrattrrandom - rnhere
indbinpos <- which(diffthr >= 0.0)
ibinsel <- indbinpos[1]
}
# Determine bin number for all other cases (probably most cases)
else {
diffthr = thrpredbin - predsam
indbinpos <- which(diffthr >= 0.0)
ibinsel <- indbinpos[1] - 1
}
# Apply correction and keep hold of results per bin for controlling
# the routine
simcell[isam] = predcell[isam] + diffbin[ibinsel]
# nsimbin[ibinsel] = nsimbin[ibinsel] + 1
# sumsimbin[ibinsel] = sumsimbin[ibinsel] + simcell[isam]
}
# Precipitation cannot be less than 0
indbel0 <- which (simcell < 0.0)
simcell[indbel0] = 0.0
# Checks
# meansimbin = sumsimbin / nsimbin
diffmean = abs (mean (simcell [indvalobs]) - mean (obscell [indvalobs]))
# diffbinmax = max (abs (meansimbin - meanobsbin), na.rm = T)
if (diffmean > 0.01) {
print (paste ("error: ", diffmean, morepred0, length(indbel0), sep = " "))
}
return(simcell)
}
biasCorrect <- function(obs, pred, nbin = 100, epsilon = 0.00001) {
## Make uniform and match the selections
dates1<-as.character(as.POSIXct(obs$Dates$start, tz = "GMT"))
dates2<-as.character(as.POSIXct(pred$Dates$start, tz = "GMT"))
overlapping_dates<-as.Date(intersect(dates1, dates2))
element<-is.element(dates2,dates1)
rm(dates2,dates1)
obs_new<-obs
ntime<-dim(pred$Data)[2]
nlat<-dim(obs$Data)[2]
nlon<-dim(obs$Data)[3]
print(dim(obs_new$Data))
obs_new$Data <- array(NA, dim = c(ntime, nlat, nlon))
print(dim(obs_new$Data))
obs_new$Data[which(element==TRUE),,]<-obs$Data[]
# obs_new$Data[,10,10]
# obs$Data[,10,10]
obs<-obs_new
#############
obssel <- obs$Data
predsel <- pred$Data
RData_BC <- pred
dimpred <- dim(predsel)
nmem <- dimpred[1]
nsam <- dimpred[2]
nlat <- dimpred[3]
nlon <- dimpred[4]
npred = nmem * nsam
obsselmem <- array(NA, dim = c(npred, nlat, nlon))
for (imem in (1:nmem)) {
ib = (imem - 1) * nsam + 1
ie = ib + nsam - 1
obsselmem[ib:ie, , ] <- obssel
}
########## START
minobs = 28 * 29
minobs = minobs * nmem
for (ilon in (1:nlon)) {
#ilon<-30
print (paste (ilon, Sys.time(), sep = " "))
for (ilat in (1:nlat)) {
#ilat<-40
# Selection of all data for a single cell
obscell <- obsselmem[ ,ilat, ilon]
predcell <- aperm( predsel[ , ,ilat, ilon])
rm(nmem,nsam,ib,ie,imem,obssel)
simcell <- biasCorrectCell(obscell, predcell, npred = npred, minobs = minobs, nbin = nbin, epsilon = epsilon)
RData_BC$Data[,,ilat,ilon]<-simcell
}