WSL/SLF GitLab Repository

Commit ffcee683 authored by Adrien Michel's avatar Adrien Michel
Browse files

Majore coode reorganisation, removal or R implementation, update of assessment functions

parent 4a6e2e87
......@@ -8,7 +8,7 @@ export(applyQMCpp)
export(applyQMCppDoy)
export(assessment.plot)
export(assessment.plot.calib)
export(compute.data.season.diff.calib)
export(compute.data.season.mean)
export(compute.quantiles)
export(compute.spatial.pools)
export(compute.temporal.pools)
......@@ -19,16 +19,16 @@ export(get.filter.mask.vectdata)
export(get.mask.na)
export(get.ovelrapp.mask)
export(get.overlapp.netcdf)
export(ma)
export(mask.vectdata)
export(plt.season.diff.calib)
export(plt.seasonnal.summary)
export(plt.seasonnal.summary.calib)
export(plt.seasonnal.summary.details)
export(plt.seasonnal.summary.details.calib)
export(plt.seasonnal.summary.maps.calib)
export(plt.snow.days.comparison)
export(plt.snow.days.comparison.calib)
export(plt.snow.days.maps)
export(plt.snow.days.maps.calib)
export(plt.snow.season.length)
export(plt.snow.season.length.calib)
export(read.asc)
......@@ -37,14 +37,12 @@ export(read.mask.tif)
export(read.swe)
export(read.tif)
export(snow.days.comparison)
export(snow.days.comparison.calib.seasonnal)
export(snow.days.comparison.do)
export(snow.days.comparison.seasonnal)
export(snow.days.comparison.seasonnal.calib)
export(snow.season.length)
export(snow.season.length.do)
export(snowDaysComparisonCpp)
export(snowSeasonLengthCpp)
export(train.model)
export(train.run.assess.model)
export(vectorize.list)
importFrom(Rcpp,evalCpp)
useDynLib(SnowQM, .registration=TRUE)
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
################################################################################
############################# INNER FUNCTIONS ##################################
################################################################################
################################################################################
############################# SINGLE MODEL RUN #################################
################################################################################
################################################################################
####################### CALIBRATION AND VALIDATION #############################
################################################################################
################################################################################
#' XXXXX
#'
#' \code{plt.snow.days.comparison} XXX
#' \code{\link{barplot()}}) XXXX
#'
#' @param data.target The target dataset (a vectdata object, see
#' \code{\link{vectorize.list()}})
#' @param data.training The training dataset (a vectdata object, see
#' \code{\link{vectorize.list()}})
#' @param data.corrected the corrected dataset (a vectdata object, see
#' \code{\link{vectorize.list()}})
#' @param out.plot.dir Path in which the RDS files should be saved. Three files are
#' created: snow_season_target, snow_season_training, snow_season_corrected
#' @param ncores Number of cores to use for parallel computation, using ncores =
#' 1 defaults to serial computing
#' @param swe.threshold The threshold SWE value to define snow days
#'
#' @author Adrien Michel, 2021 (WSL/SLF)
plt.spatial.correlation <- function(data.target.calib, data.training.calib,
data.corrected.calib, data.target.valid,
data.training.valid, data.corrected.valid,
out.plot.dir, ncores = 1, dist.threshold = 5) {
if(ncores > 1){
`%dopar%` <- foreach::`%dopar%`
doParallel::registerDoParallel(ncores)
} else {
`%dopar%` <- foreach::`%do%`
}
time <- Sys.time()
x.coords.vect <- data.target.calib$x.coords[data.target.calib$x.indices]
y.coords.vect <- data.target.calib$y.coords[data.target.calib$y.indices]
x <- matrix(0,ncol=length(data.target.calib$x.indices),nrow=length(data.target.calib$x.indices))
y <- matrix(0,ncol=length(data.target.calib$y.indices),nrow=length(data.target.calib$y.indices))
for(i in 1:length(data.target.calib$x.indices)){
x[,i]=((x.coords.vect-x.coords.vect[i])/1000)^2
y[,i]=((y.coords.vect-y.coords.vect[i])/1000)^2
}
indices=list()
for (i in 1:length(data.target.calib$x.indices)){
indices[[i]] <- which(x[,pixel]+y[,pixel] >0 & x[,pixel]+y[,pixel]<=dist.threshold^2)
}
print(Sys.time()-time)
get.correlation <- function(data,indices) {
data.tmp <- matrix(apply(data$data, 1, mean, na.rm=T))
diff=rep(0,times=length(data$x.indices))
for (pixel in 1:length(data$x.indices)){
diff[pixel] = ((data.tmp[pixel]))/mean(data.tmp[indices[[pixel]]])
}
return(matrix(diff))
}
data.target.calib$correlation=get.correlation(data.target.calib,indices)
plot(data.target.calib,var = "correlation")
data.target.valid$correlation=get.correlation(data.target.valid,indices)
plot(data.target.valid,var = "correlation")
data.tmp <- matrix(apply(data.training.calib$data, 1, mean, na.rm=T))
correlation = foreach::foreach (pixel=1:length(data.training.calib$data[,1])) %dopar% {
dist <- as.matrix(sqrt((data.training.calib$x.coords[data.training.calib$x.indices[pixel]]-
data.training.calib$x.coords[data.training.calib$x.indices])^2+
(data.training.calib$y.coords[data.training.calib$y.indices[pixel]]-
data.training.calib$y.coords[data.training.calib$y.indices])^2))/1000
indices <- which(dist>0 & dist<=dist.threshold)
diff=(data.tmp[pixel]-data.tmp[indices])/(data.tmp[pixel]+data.tmp[indices])
diff[(data.tmp[pixel]+data.tmp[indices]) == 0] <- 0
mean(diff)
}
data.training.calib$correlation=matrix(unlist(correlation))
data.tmp <- matrix(apply(data.corrected.calib$data, 1, mean, na.rm=T))
correlation = foreach::foreach (pixel=1:length(data.corrected.calib$data[,1])) %dopar% {
dist <- as.matrix(sqrt((data.corrected.calib$x.coords[data.corrected.calib$x.indices[pixel]]-
data.corrected.calib$x.coords[data.corrected.calib$x.indices])^2+
(data.corrected.calib$y.coords[data.corrected.calib$y.indices[pixel]]-
data.corrected.calib$y.coords[data.corrected.calib$y.indices])^2))/1000
indices <- which(dist>0 & dist<=dist.threshold)
diff=(data.tmp[pixel]-data.tmp[indices])/(data.tmp[pixel]+data.tmp[indices])
diff[(data.tmp[pixel]+data.tmp[indices]) == 0] <- 0
mean(diff)
}
data.corrected.calib$correlation=matrix(unlist(correlation))
data.tmp <- matrix(apply(data.target.valid$data, 1, mean, na.rm=T))
correlation = foreach::foreach (pixel=1:length(data.target.valid$data[,1])) %dopar% {
dist <- as.matrix(sqrt((data.target.valid$x.coords[data.target.valid$x.indices[pixel]]-
data.target.valid$x.coords[data.target.valid$x.indices])^2+
(data.target.valid$y.coords[data.target.valid$y.indices[pixel]]-
data.target.valid$y.coords[data.target.valid$y.indices])^2))/1000
indices <- which(dist>0 & dist<=dist.threshold)
diff=(data.tmp[pixel]-data.tmp[indices])/(data.tmp[pixel]+data.tmp[indices])
diff[(data.tmp[pixel]+data.tmp[indices]) == 0] <- 0
mean(diff)
}
data.target.valid$correlation=matrix(unlist(correlation))
data.tmp <- matrix(apply(data.training.valid$data, 1, mean, na.rm=T))
correlation = foreach::foreach (pixel=1:length(data.training.valid$data[,1])) %dopar% {
dist <- as.matrix(sqrt((data.training.valid$x.coords[data.training.valid$x.indices[pixel]]-
data.training.valid$x.coords[data.training.valid$x.indices])^2+
(data.training.valid$y.coords[data.training.valid$y.indices[pixel]]-
data.training.valid$y.coords[data.training.valid$y.indices])^2))/1000
indices <- which(dist>0 & dist<=dist.threshold)
diff=(data.tmp[pixel]-data.tmp[indices])/(data.tmp[pixel]+data.tmp[indices])
diff[(data.tmp[pixel]+data.tmp[indices]) == 0] <- 0
mean(diff)
}
data.training.valid$correlation=matrix(unlist(correlation))
data.tmp <- matrix(apply(data.corrected.valid$data, 1, mean, na.rm=T))
correlation = foreach::foreach (pixel=1:length(data.corrected.valid$data[,1])) %dopar% {
dist <- as.matrix(sqrt((data.corrected.valid$x.coords[data.corrected.valid$x.indices[pixel]]-
data.corrected.valid$x.coords[data.corrected.valid$x.indices])^2+
(data.corrected.valid$y.coords[data.corrected.valid$y.indices[pixel]]-
data.corrected.valid$y.coords[data.corrected.valid$y.indices])^2))/1000
indices <- which(dist>0 & dist<=dist.threshold)
diff=(data.tmp[pixel]-data.tmp[indices])/(data.tmp[pixel]+data.tmp[indices])
diff[(data.tmp[pixel]+data.tmp[indices]) == 0] <- 0
mean(diff)
}
data.corrected.valid$correlation=matrix(unlist(correlation))
pdf(paste0(out.plot.dir,"/snow_height_spatial_correlation.pdf"),width=12,height=8)
par(mar=c(3,3,2,5), mgp=c(1.7,0.5,0),mfrow=c(2,3))
plot(data.target.calib,var="correlation", main=paste("Target data calib"),zlim=c(-1,1))
plot(data.training.calib,var="correlation", main=paste("Traning data calib"),zlim=c(-1,1))
plot(data.corrected.calib,var="correlation", main=paste("Corrected data calib"),zlim=c(-1,1))
plot(data.target.valid,var="correlation", main=paste("Target data valid"),zlim=c(-1,1))
plot(data.training.valid,var="correlation", main=paste("Training data valid"),zlim=c(-1,1))
plot(data.corrected.valid,var="correlation", main=paste("Corrected data valid"),zlim=c(-1,1))
data.training.calib$correlation.rel <- data.training.calib$correlation-data.target.calib$correlation
data.corrected.calib$correlation.rel <- data.corrected.calib$correlation-data.target.calib$correlation
data.training.valid$correlation.rel <- data.training.valid$correlation-data.target.valid$correlation
data.corrected.valid$correlation.rel <- data.corrected.valid$correlation-data.target.valid$correlation
par(mfrow=c(2,2))
me <- round(mean(data.training.calib$correlation.rel),4)
mae <- round(mean(abs(data.training.calib$correlation.rel)),3)
plot(data.training.calib,var="correlation.rel", main=paste("Traning data calib", me, mae), zlim=c(-1,1))
me <- round(mean(data.corrected.calib$correlation.rel),4)
mae <- round(mean(abs(data.corrected.calib$correlation.rel)),3)
plot(data.corrected.calib,var="correlation.rel", main=paste("Corrected data calib", me, mae),zlim=c(-1,1))
me <- round(mean(data.training.valid$correlation.rel),4)
mae <- round(mean(abs(data.training.valid$correlation.rel)),3)
plot(data.training.valid,var="correlation.rel", main=paste("Target data valid", me, mae),zlim=c(-1,1))
me <- round(mean(data.corrected.valid$correlation.rel),4)
mae <- round(mean(abs(data.corrected.valid$correlation.rel)),3)
plot(data.corrected.valid,var="correlation.rel", main=paste("Corrected data valid", me, mae),zlim=c(-1,1))
dev.off()
if(ncores > 1){
doParallel::stopImplicitCluster()
}
}
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
......@@ -5,7 +5,8 @@
#' @export
read.swe <- function (data.path, start=NULL, end=NULL, count=NULL,
data.dim="SWE", x.coords.dim="chx", y.coords.dim="chy", mask=NULL,
dem.path, slope.path, azi.path, x.shift=0, y.shift=0, silent=FALSE) {
dem.path, slope.path, azi.path, curvature.path, x.shift=0, y.shift=0,
silent=FALSE) {
ncdata <- read.nc.data(data.path, start, end, count, silent, data.dim,
x.coords.dim, y.coords.dim)
......@@ -18,6 +19,7 @@ read.swe <- function (data.path, start=NULL, end=NULL, count=NULL,
dem <- vectorize.list(read.asc(dem.path, x.shift, y.shift))
slope <- vectorize.list(read.asc(slope.path, x.shift, y.shift))
azi <- vectorize.list(read.asc(azi.path, x.shift, y.shift))
curvature <- vectorize.list(read.asc(curvature.path, x.shift, y.shift))
if(!identical(dem$x.coords,vectdata$x.coords) |
!identical(dem$y.coords,vectdata$y.coords)){
......@@ -43,6 +45,9 @@ read.swe <- function (data.path, start=NULL, end=NULL, count=NULL,
coords.2 <- paste(slope$x.coords[slope$x.indices], slope$y.coords[slope$y.indices],
sep="_")
vectdata$slope <- as.matrix(slope$data[match(coords.1, coords.2)])
coords.2 <- paste(slope$x.coords[slope$x.indices], slope$y.coords[slope$y.indices],
sep="_")
vectdata$curvature <- as.matrix(curvature$data[match(coords.1, coords.2)])
if(!is.null(mask)){
print("[i] Start masking the data")
......@@ -73,9 +78,6 @@ read.nc.data <- function(file.path, start=NULL, end=NULL, count=NULL,
ncdata$time.date <- round(as.POSIXct(ncdf4.helpers::nc.get.time.series(nc,
time.dim.name="time")),"days")
ncdata$lat <- ncdf4::ncvar_get(nc, "lat")
ncdata$lon <- ncdf4::ncvar_get(nc, "lon")
if ((!is.null(start) || !is.null(end)) && is.null(count)) {
subset <- cut.time.start.end(ncdata$time.date, start,end)
ncdata$time <- ncdata$time[subset]
......
......@@ -128,8 +128,8 @@ set.na.vectdata <- function(vectdata, val.to.remove){
#' @param filename A
#' @param ... A
#' @export
read.mask.asc <- function(filename, val.to.remove=NULL, x.shift=NULL,
y.shift=NULL){
read.mask.asc <- function(filename, val.to.remove=NULL, x.shift=0,
y.shift=0){
asc.list <- read.asc(filename, x.shift, y.shift)
mask.vectdata <- vectorize.list(asc.list, "data")
if(!is.null(val.to.remove)){
......@@ -144,8 +144,8 @@ read.mask.asc <- function(filename, val.to.remove=NULL, x.shift=NULL,
#' @param filename A
#' @param ... A
#' @export
read.mask.tif <- function(filename, val.to.remove=NULL, x.shift=NULL,
y.shift=NULL){
read.mask.tif <- function(filename, val.to.remove=NULL, x.shift=0,
y.shift=0){
tif.list <- read.tif(filename, x.shift, y.shift)
mask.vectdata <- vectorize.list(tif.list, "data")
if(!is.null(val.to.remove)){
......
#' Train, apply, and produce assessment plot of the QM model
#' Apply the trained QM model to a given dataset
#'
#' \code{train.run.assess.model} applies a trained QM model (see
#' \code{apply.model} applies a trained QM model (see
#' \code{\link{compute.quantiles()}}) to the dataset given as input.
#'
#' @param data.training.file Dataset on which the QM model should be applied. Dataset
#' @param in.data.file Dataset on which the QM model should be applied. Dataset
#' should be a RDS file in a vectdata format (see \code{\link{vectorize.list()}})
#' @param out.data.file Path of the file to write the corrected dataset. Dataset
#' will be a RDS file in a vectdata format (see \code{\link{vectorize.list()}})
#' @param source.path Directory containing the directories with quantile distributions of the
#' training and target datasets (see \code{\link{vectorize.list()}})
#' @param type R or cpp to choose between the R or the C++ implementation
#' (see \code{\link{apply.qm.internal()}} and \code{\link{applyQMCpp()}}).
#' C++ implementation is faster and recommended, given the package has been
#' installed with OpenMP support.
#' @param ncores number of cores to be used for the quantile mapping part
#' @param ncores.pp number of cores to be used for the post processing part, this
#' part in only implemented in R and memory issues might arise
#' @param swe.threshold The threshold value to define snow day, in mSWE, needed
#' if fun = num.snow.days
#' @param zlim Limit for the map plots in \code{\link{plt.seasonnal.summary.details()}}
#' @param ncores number of cores to be used
#'
#' @author Adrien Michel, 2021 (WSL/SLF)
#' @export
apply.model <- function(in.data.file, out.data.file, quantile.path, ncores)
{
print("[I] Creating output directory")
dir.create(dirname(out.data.file), showWarnings = FALSE)
train.run.assess.model <- function(data.training.file, data.target.file, params.qm,
out.dir, type, ncores, ncores.pp, swe.threshold,
zlim = "auto") {
start.time <- Sys.time()
out.data.file <- sub('\\.RDS$', '', basename(data.training.file))
out.data.file <- paste0(out.dir,"/",out.data.file,"_corrected.RDS")
train.model(data.training.file, data.target.file, params.qm = params.qm,
out.dir = out.dir, type = type, ncores = ncores)
print("[I] Reading data")
data <- readRDS(in.data.file)
apply.model(in.data.file = data.training.file,
out.data.file = out.data.file, quantile.path = out.dir,
type = type, ncores = ncores)
print("[I] Computing days of the year")
doys=as.integer(strftime(data$time.date, format = "%j"))
doys.unique=sort(unique(doys))
`%do%` <- foreach::`%do%`
doy.selection = foreach::foreach(doy=doys.unique) %do% {
which(doys == doy)
}
names(doy.selection)=doys.unique
assessment.plot(data.training.file, data.target.file,
data.corrected.file = out.data.file,
out.path = out.dir, swe.threshold = swe.threshold,
ncores = ncores.pp, zlim)
print("[I] Correcting data")
start.time <- Sys.time()
corrected.swe <- applyQMCpp(data$data,
doyArray = doys,
sourcePath = paste0(quantile.path,"/training_quantiles"),
targetPath = paste0(quantile.path,"/target_quantiles"),
ncores = ncores)
data.corrected <- data
data.corrected$data <- matrix(corrected.swe, ncol=dim(data$data)[2])
print(paste("[I] Total run in",
print(paste("[I] Correction computed in",
as.numeric(round(difftime(Sys.time(),start.time, units = "mins"),2)),
"minutes"))
saveRDS(data.corrected,out.data.file)
print("[I] Cleaning memory")
rm(list=ls())
invisible(gc())
}
\ No newline at end of file
}
......@@ -234,24 +234,13 @@ compute.spatial.pools <- function(data, params.qm, ncores = 1,
}
compute.quantiles.internal <- function(data, spatial.pools, probs)
{
pixels <- seq(dim(data)[1])
out <- matrix(ncol=101, nrow=length(pixels))
for(p in pixels){
out[p,] <- quantile(data[spatial.pools[[p]],], probs = probs, na.rm = TRUE)
}
return(out)
}
#' Title
#'
#' @param x A
#' @param ... A
#' @export
compute.quantiles <- function(data, temporal.pools, spatial.pools, output.path,
type, ncores=1, probs=seq(0, 1, 0.01))
ncores=1, probs=seq(0, 1, 0.01))
{
if(!dir.exists(output.path)) {
......@@ -260,33 +249,14 @@ compute.quantiles <- function(data, temporal.pools, spatial.pools, output.path,
stop(paste0("Output directory '", output.path, "' seems to be not writable"))
}
if(type == "cpp") {
computeQuantilesCpp(swe=data$data, timePools=temporal.pools,
spatialPools=spatial.pools, probs=probs,
outputPath=output.path, ncores=ncores)
} else if (type == "R"){
if(ncores>1) {
`%dopar%` <- foreach::`%dopar%`
doParallel::registerDoParallel(ncores)
dummy <- foreach::foreach(doy=c(1:366), .errorhandling="stop") %dopar% {
PDF <- compute.quantiles.internal(data$data[,temporal.pools[[doy]]],
spatial.pools, probs)
write.table(round(PDF, digits=3), file = paste0(output.path,"/quantiles_", doy, ".txt"),
row.names = FALSE, col.names = FALSE)
}
doParallel::stopImplicitCluster()
} else {
`%do%` <- foreach::`%do%`
dummy <- foreach::foreach(doy=c(1:366), .errorhandling="stop") %do% {
PDF <- compute.quantiles.internal(data$data[,temporal.pools[[doy]]],
spatial.pools,probs)
write.table(round(PDF, digits=3), file = paste0(output.path,"/quantiles_", doy, ".txt"),
row.names = FALSE, col.names = FALSE)
}
}
} else {
stop("Type sould be either cpp or R")
}
computeQuantilesCpp(swe=data$data,
timePools=temporal.pools,
spatialPools=spatial.pools,
xCoords=round(data$x.coords[data$x.indices]),
yCoords=round(data$x.yoords[data$xy.indices]),
probs=probs,
outputPath=output.path, ncores=ncores)
}
#' Title
......@@ -294,9 +264,8 @@ compute.quantiles <- function(data, temporal.pools, spatial.pools, output.path,
#' @param x A
#' @param ... A
#' @export
train.model <- function (data.training.file , data.target.file, params.qm,
out.quantiles.dir, out.plot.dir,
type = "cpp", ncores = 1){
train.model <- function(data.training.file , data.target.file, params.qm,
out.quantiles.dir, out.plot.dir, ncores = 1){
print("[I] Reading training data")
data.training <- readRDS(data.training.file)
......@@ -360,7 +329,7 @@ train.model <- function (data.training.file , data.target.file, params.qm,
temporal.pools = temporal.pools,
spatial.pools = spatial.pools,
output.path = paste0(out.quantiles.dir,"/training_quantiles/"),
type = type, ncores = ncores)
ncores = ncores)
print(paste("[I] Training quantiles computed in",
as.numeric(round(difftime(Sys.time(),start.time, units = "mins"),2)),
"minutes"))
......@@ -371,7 +340,7 @@ train.model <- function (data.training.file , data.target.file, params.qm,
temporal.pools = temporal.pools,
spatial.pools = spatial.pools,
output.path = paste0(out.quantiles.dir,"/target_quantiles/"),
type = type,ncores = ncores)
ncores = ncores)
print(paste("[I] Target quantiles computed in",
as.numeric(round(difftime(Sys.time(),start.time, units = "mins"),2)),
"minutes"))
......
apply.qm.internal <- function(data, doy.selection, source.path, target.path,
type, ncores)
{
if(ncores > 1){
doParallel::registerDoParallel(ncores)
`%do%` <- foreach::`%dopar%`
} else{
`%do%` <- foreach::`%do%`
}
data.corrected <- data
pixels <- c(1:length(data$data[,1]))
for(doy in names(doy.selection)) {
doy.subset <- doy.selection[[doy]]
PDF.source <- as.matrix(data.table::fread(paste0(source.path,
"/quantiles_",doy,".txt"),
header=F))
PDF.dest <- as.matrix(data.table::fread(paste0(target.path,
"/quantiles_",doy,".txt"),
header=F))
if(length(PDF.source[,1]) > length(pixels)) {
stop(paste0("File ",source.path, "/quantiles_",
doy, ".txt has the wrong number of pixels"))
}
if(length(PDF.dest[,1]) > length(pixels)) {
stop(paste0("File ",target.path, "/quantiles_",
doy, ".txt has the wrong number of pixels"))
}
correction = foreach::foreach(p = pixels) %do% {
if(length(PDF.source[p,]) != length(PDF.dest[p,])){
stop("Source and target quantile files have different number of qunatiles")
}
corrected <- rep(NA,length(doy.subset))
i <- 1
for(d in doy.subset) {
quantile <- which(PDF.source[p,] > data$data[p,d])[1]-1
quantile <- bound.quantile(quantile, length(PDF.source[p,]))
corrected[i] <- data$data[p,d] - PDF.source[p,quantile] + PDF.dest[p,quantile]
i <- i+1
}
as.vector(corrected)
}
for(p in pixels) {
i=1
for(s in sel){
data.corrected$data[p,s]=correction[[p]][i]
i=i+1
}
}
}
if(ncores > 1){
doParallel::stopImplicitCluster()
}
return(data.corrected)
}
#' Apply the trained QM model to a given dataset
#'
#' \code{apply.model} applies a trained QM model (see
#' \code{\link{compute.quantiles()}}) to the dataset given as input.
#'
#' @param in.data.file Dataset on which the QM model should be applied. Dataset
#' should be a RDS file in a vectdata format (see \code{\link{vectorize.list()}})
#' @param out.data.file Path of the file to write the corrected dataset. Dataset
#' will be a RDS file in a vectdata format (see \code{\link{vectorize.list()}})
#' @param source.path Directory containing the directories with quantile distributions of the
#' training and target datasets (see \code{\link{vectorize.list()}})
#' @param type R or cpp to choose between the R or the C++ implementation
#' (see \code{\link{apply.qm.internal()}} and \code{\link{applyQMCpp()}}).
#' C++ implementation is faster and recommended, given the package has been
#' installed with OpenMP support.
#' @param ncores number of cores to be used
#'
#' @author Adrien Michel, 2021 (WSL/SLF)
#' @export
apply.model <- function(in.data.file, out.data.file, quantile.path, type, ncores)
{
print("[I] Creating output directory")
dir.create(dirname(out.data.file), showWarnings = FALSE)
print("[I] Reading data")
data <- readRDS(in.data.file)
print("[I] Computing days of the year")
doys=as.integer(strftime(data$time.date, format = "%j"))
doys.unique=sort(unique(doys))
`%do%` <- foreach::`%do%`
doy.selection = foreach::foreach(doy=doys.unique) %do% {
which(doys == doy)
}
names(doy.selection)=doys.unique
print("[I] Correcting data")
start.time <- Sys.time()
if(type == "cpp") {
corrected.swe <- applyQMCpp(data$data,
doyArray = doys,
sourcePath = paste0(quantile.path,"/training_quantiles"),
targetPath = paste0(quantile.path,"/target_quantiles"),
ncores = ncores)
data.corrected <- data