From a881e95c8815b5de6f036b2c4f062575475ad938 Mon Sep 17 00:00:00 2001 From: Hollos Roland Date: Wed, 23 Sep 2020 08:13:12 +0200 Subject: [PATCH 1/7] parallel muso --- RBBGCMuso/NAMESPACE | 3 + RBBGCMuso/R/calibrateMuso.R | 390 +++++++++++++++++++++++++++++++++ RBBGCMuso/R/calibration.R | 84 ++++--- RBBGCMuso/R/changeMuso.R | 11 +- RBBGCMuso/R/setupMuso.R | 42 ++-- RBBGCMuso/R/thread-1 | 1 + RBBGCMuso/R/thread-2 | 1 + RBBGCMuso/R/thread-3 | 1 + RBBGCMuso/R/thread-4 | 1 + RBBGCMuso/man/calibrateMuso.Rd | 38 ++++ RBBGCMuso/man/changemulline.Rd | 12 +- RBBGCMuso/man/musoGlue.Rd | 9 +- RBBGCMuso/man/optiMuso.Rd | 5 +- calibrateMuso.R | 295 +++++++++++++++++++++++++ parallelRun.R | 90 ++++++++ 15 files changed, 914 insertions(+), 69 deletions(-) create mode 100644 RBBGCMuso/R/calibrateMuso.R create mode 100644 RBBGCMuso/R/thread-1 create mode 100644 RBBGCMuso/R/thread-2 create mode 100644 RBBGCMuso/R/thread-3 create mode 100644 RBBGCMuso/R/thread-4 create mode 100644 RBBGCMuso/man/calibrateMuso.Rd create mode 100644 calibrateMuso.R create mode 100644 parallelRun.R diff --git a/RBBGCMuso/NAMESPACE b/RBBGCMuso/NAMESPACE index 4eb00d3..76126a4 100644 --- a/RBBGCMuso/NAMESPACE +++ b/RBBGCMuso/NAMESPACE @@ -1,7 +1,9 @@ # Generated by roxygen2: do not edit by hand export(calibMuso) +export(calibrateMuso) export(changeMusoC) +export(changemulline) export(checkMeteoBGC) export(cleanupMuso) export(compareMuso) @@ -54,6 +56,7 @@ importFrom(dplyr,summarize) importFrom(dplyr,tbl_df) importFrom(ecmwfr,wf_request) importFrom(ecmwfr,wf_set_key) +importFrom(future,future) importFrom(ggplot2,aes) importFrom(ggplot2,aes_string) importFrom(ggplot2,element_blank) diff --git a/RBBGCMuso/R/calibrateMuso.R b/RBBGCMuso/R/calibrateMuso.R new file mode 100644 index 0000000..0deefa5 --- /dev/null +++ b/RBBGCMuso/R/calibrateMuso.R @@ -0,0 +1,390 @@ +#' calibrateMuso +#' +#' This funtion uses the Monte Carlo technique to uniformly sample the parameter space from user defined parameters of the Biome-BGCMuSo model. The sampling algorithm ensures that the parameters are constrained by the model logic which means that parameter dependencies are fully taken into account (parameter dependency means that e.g leaf C:N ratio must be smaller than C:N ratio of litter; more complicated rules apply to the allocation parameters where the allocation fractions to different plant compartments must sum up 1). This function implements a mathematically correct solution to provide uniform distriution for all selected parameters. +#' @author Roland HOLLOS +#' @importFrom future future +#' @export +calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, + endDate = NULL, formatString = "%Y-%m-%d", + dataVar, outLoc = "./calib", + preTag = "cal-", settings = setupMuso(), + outVars = NULL, iterations = 100, + skipSpinup = TRUE, plotName = "calib.jpg", + modifyOriginal=TRUE, likelihood, uncertainity = NULL, + naVal = NULL, postProcString = NULL, + thread_prefix="thread", numCores = (parallel::detectCores()-1), pb = txtProgressBar(min=0, max=iterations, style=3), + maxLikelihoodEpc=TRUE, + pbUpdate = setTxtProgressBar, method="GLUE",lg = FALSE, w=NULL, ...){ + + + file.remove(list.files(path = settings$inputLoc, pattern="progress.txt", recursive = TRUE)) + file.remove(list.files(path = settings$inputLoc, pattern="preservedCalib.csv", recursive = TRUE)) + unlink("thread",recursive=TRUE) + + # ____ _ _ _ _ + # / ___|_ __ ___ __ _| |_ ___ | |_| |__ _ __ ___ __ _ __| |___ + # | | | '__/ _ \/ _` | __/ _ \ | __| '_ \| '__/ _ \/ _` |/ _` / __| + # | |___| | | __/ (_| | || __/ | |_| | | | | | __/ (_| | (_| \__ \ + # \____|_| \___|\__,_|\__\___| \__|_| |_|_| \___|\__,_|\__,_|___/ + + + + copyToThreadDirs(thread_prefix, numcores = numCores, runDir = settings$inputLoc) + + # ____ _ _ _ + # | _ \ _ _ _ __ | |_| |__ _ __ ___ __ _ __| |___ + # | |_) | | | | '_ \ | __| '_ \| '__/ _ \/ _` |/ _` / __| + # | _ <| |_| | | | | | |_| | | | | | __/ (_| | (_| \__ \ + # |_| \_\\__,_|_| |_| \__|_| |_|_| \___|\__,_|\__,_|___/ + + threadCount <- distributeCores(iterations, numCores) + + fut <- lapply(1:numCores, function(i) { + # browser() + future({ + tryCatch(musoSingleThread(measuredData, parameters, startDate, + endDate, formatString, + dataVar, outLoc, + preTag, settings, + outVars, iterations = threadCount[i], + skipSpinup, plotName, + modifyOriginal, likelihood, uncertainity, + naVal, postProcString, i), error = function(e){ + writeLines(as.character(iterations),"progress.txt") + }) + + # musoSingleThread(measuredData, parameters, startDate, + # endDate, formatString, + # dataVar, outLoc, + # preTag, settings, + # outVars, iterations = threadCount[i], + # skipSpinup, plotName, + # modifyOriginal, likelihood, uncertainity, + # naVal, postProcString, i) + }) + }) + + # __ ___ _ _ + # \ \ / / |__ __ _| |_ ___| |__ _ __ _ __ ___ ___ ___ ___ ___ + # \ \ /\ / /| '_ \ / _` | __/ __| '_ \ | '_ \| '__/ _ \ / __/ _ \/ __/ __| + # \ V V / | | | | (_| | || (__| | | | | |_) | | | (_) | (_| __/\__ \__ \ + # \_/\_/ |_| |_|\__,_|\__\___|_| |_| | .__/|_| \___/ \___\___||___/___/ + # |_| + + getProgress <- function(){ + # threadfiles <- list.files(settings$inputLoc, pattern="progress.txt", recursive = TRUE) + threadfiles <- list.files(pattern="progress.txt", recursive = TRUE) + if(length(threadfiles)==0){ + return(0) + } else { + sum(sapply(threadfiles, function(x){ + partRes <- readLines(x) + if(length(partRes)==0){ + return(0) + } else { + return(as.numeric(partRes)) + } + + })) + + } + } + + progress <- 0 + while(progress < iterations){ + Sys.sleep(1) + progress <- tryCatch(getProgress(), error=function(e){progress}) + pbUpdate(pb,as.numeric(progress)) + } + close(pb) + + # ____ _ _ + # / ___|___ _ __ ___ | |__ (_)_ __ ___ + # | | / _ \| '_ ` _ \| '_ \| | '_ \ / _ \ + # | |__| (_) | | | | | | |_) | | | | | __/ + # \____\___/|_| |_| |_|_.__/|_|_| |_|\___| + + resultFiles <- list.files(pattern="preservedCalib.*csv$",recursive=TRUE) + res0 <- read.csv(grep("thread_1/",resultFiles, value=TRUE),stringsAsFactors=FALSE) + resultFilesSans0 <- grep("thread_1/", resultFiles, value=TRUE, invert=TRUE) + # results <- do.call(rbind,lapply(resultFilesSans0, function(f){read.csv(f, stringsAsFactors=FALSE)})) + resultsSans0 <- lapply(resultFilesSans0, function(f){read.csv(f, stringsAsFactors=FALSE, header=FALSE)}) + resultsSans0 <- do.call(rbind,resultsSans0) + colnames(resultsSans0) <- colnames(res0) + results <- (rbind(res0,resultsSans0)) + + switch(method, + "GLUE"={ + musoGlue(results, parameters=parameters,settings=settings, w=w, lg=lg) + }, + "agromo"={ + liks <- results[,sprintf("%s_likelihood",names(likelihood))] + epcIndexes <- value(fut[[1]], stdout = FALSE, signal=FALSE) + epcVals <- results[which.max(liks),1:length(epcIndexes)] + changemulline(filePaths= settings$epcInput[2], epcIndexes, epcVals, src = settings$epcInput[2], outFiles = "maxLikelihood_epc.epc") + names(epcVals) <- epcIndexes + xdate <- as.Date(measuredData$date) + meanM <- measuredData[,sprintf("mean.%s", names(likelihood))] + minsd <- meanM - measuredData[,sprintf("sd.%s", names(likelihood)[1])] + maxsd <- meanM + measuredData[,sprintf("sd.%s", names(likelihood)[1])] + minM <- measuredData[,sprintf("min.%s", names(likelihood)[1])] + maxM <- measuredData[,sprintf("max.%s", names(likelihood)[1])] + plot(xdate, minM, type="l", xlab=NA, ylim=c(min(minM)*0.8, max(maxM)*1.1), ylab = names(likelihood)[1]) + lines(xdate, maxM) + polygon(c(xdate,rev(xdate)),c(minM,rev(maxM)), col="gray",border=NA) + lines(xdate, minsd) + lines(xdate, maxsd) + polygon(c(xdate,rev(xdate)),c(minsd,rev(maxsd)), col="gray30",border=NA) + points(xdate,meanM) + + varIndex <- match(as.character(dataVar),settings$dailyVarCodes) + apriori <- calibMuso(settings) + modDates <- as.Date(row.names(apriori), format="%d.%m.%Y") + lines(modDates, apriori[,varIndex],col="brown") + calibrated <- calibMuso(settings, calibrationPar = as.numeric(names(epcVals)), parameters=epcVals) + lines(modDates, calibrated[,varIndex],col="blue") + + }, + stop(sprintf("method: %s not found, please choose from {GLUE, agromo}. See more about this in the documentation of the function!", method)) + ) + # Here starts maxLikelihoodAgroMo: parameters + + + # Here ends maxLikelihoodAgromo + + # return(epcVals) + # ____ _ _ _ _____ + # / ___| | | | | | ____| + # | | _| | | | | | _| + # | |_| | |__| |_| | |___ + # \____|_____\___/|_____| + + # musoGlue("preservedCalib.csv",w=w, lg = lg) +} + +copyToThreadDirs <- function(prefix="thread", numcores=parallel::detectCores()-1, runDir="."){ + dir.create(file.path(runDir,prefix), showWarnings=TRUE) + fileNames <- grep("^thread.*", list.files(runDir), value=TRUE, invert=TRUE) + invisible(sapply(1:numcores,function(corenum){ + threadDir <- file.path(runDir,prefix,paste0(prefix,"_",corenum)) + dir.create(threadDir, showWarnings=FALSE) + file.copy(from=fileNames,to=threadDir, overwrite=FALSE) + })) +} + +musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, + endDate = NULL, formatString = "%Y-%m-%d", + dataVar, outLoc = "./calib", + preTag = "cal-", settings = setupMuso(), + outVars = NULL, iterations = 300, + skipSpinup = TRUE, plotName = "calib.jpg", + modifyOriginal=TRUE, likelihood, uncertainity = NULL, + naVal = NULL, postProcString = NULL, threadNumber) { + + setwd(paste0(settings$inputLoc, "/thread/thread_", threadNumber)) + iniFiles <- list.files(pattern=".*ini") + if(length(iniFiles)==1){ + iniFiles <- rep(iniFiles, 2) + } + settings <- setupMuso(iniInput = iniFiles) + # Exanding likelihood + likelihoodFull <- as.list(rep(NA,length(dataVar))) + names(likelihoodFull) <- names(dataVar) + if(!missing(likelihood)) { + lapply(names(likelihood),function(x){ + likelihoodFull[[x]] <<- likelihood[[x]] + }) + } + defaultLikelihood <- which(is.na(likelihood)) + if(length(defaultLikelihood)>0){ + likelihoodFull[[defaultLikelihood]] <- (function(x, y){ + exp(-sqrt(mean((x-y)^2))) + }) + } + + mdata <- measuredData + if(is.null(parameters)){ + parameters <- tryCatch(read.csv("parameters.csv", stringsAsFactor=FALSE), error = function (e) { + stop("You need to specify a path for the parameters.csv, or a matrix.") + }) + } else { + if((!is.list(parameters)) & (!is.matrix(parameters))){ + parameters <- tryCatch(read.csv(parameters, stringsAsFactor=FALSE), error = function (e){ + stop("Cannot find neither parameters file neither the parameters matrix") + }) + }} + + outLoc <- normalizePath(outLoc) + outLocPlain <- basename(outLoc) + currDir <- getwd() + + if(!dir.exists(outLoc)){ + dir.create(outLoc) + warning(paste(outLoc," is not exists, so it was created")) + } + + outLoc <- normalizePath(outLoc) + parameterNames <- parameters[,1] + pretag <- file.path(outLoc,preTag) + ##reading the original epc file at the specified + ## row numbers + print("optiMuso is randomizing the epc parameters now...",quote = FALSE) + if(iterations < 3000){ + randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = 3000) + randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),] + } else { + randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations) + } + + origEpc <- readValuesFromFile(settings$epc[2],randVals[[1]]) + partialResult <- matrix(ncol=length(randVals[[1]])+2*length(dataVar)) + colN <- randVals[[1]] + colN[match(parameters[,2],randVals[[1]])] <- parameters[,1] + colN[match(parameters[,2], randVals[[1]])[!is.na(match(parameters[,2],randVals[[1]]))]] <- parameters[,1] + colnames(partialResult) <- c(colN,sprintf("%s_likelihood",names(dataVar)), + sprintf("%s_rmse",names(dataVar))) + numParameters <- length(colN) + partialResult[1:numParameters] <- origEpc + ## Prepare the preservedCalib matrix for the faster + ## run. + + pretag <- file.path(outLoc,preTag) + + musoCodeToIndex <- sapply(dataVar,function(musoCode){ + settings$dailyOutputTable[settings$dailyOutputTable$code == musoCode,"index"] + }) + resultRange <- (numParameters + 1):(ncol(partialResult)) + ## Creating function for generating separate + ## csv files for each run + + settings$iniInput[2] %>% + (function(x) paste0(dirname(x),"/",tools::file_path_sans_ext(basename(x)),"-tmp.",tools::file_ext(x))) %>% + unlink + randValues <- randVals[[2]] + + settings$calibrationPar <- randVals[[1]] + + if(!is.null(naVal)){ + measuredData <- as.data.frame(measuredData) + measuredData[measuredData == naVal] <- NA + } + + alignIndexes <- alignMuso(settings,measuredData) + if(!is.null(uncertainity)){ + uncert <- measuredData[alignIndexes$meas,uncertainity] + } else { + uncert <- NULL + } + # browser() + if(threadNumber == 1){ + origModellOut <- calibMuso(settings=settings, silent=TRUE, skipSpinup = skipSpinup, postProcString=postProcString, modifyOriginal=modifyOriginal) + partialResult[,resultRange] <- calcLikelihoodsAndRMSE(dataVar=dataVar, + mod=origModellOut, + mes=measuredData, + likelihoods=likelihood, + alignIndexes=alignIndexes, + musoCodeToIndex = musoCodeToIndex,uncert=uncert) + write.csv(x=origModellOut, file=paste0(pretag, 1, ".csv")) + write.csv(x=partialResult, file="preservedCalib.csv",row.names=FALSE) + } + print("Running the model with the random epc values...", quote = FALSE) + + # if(!is.null(postProcString)){ + # colNumb <- length(settings$dailyVarCodes) + 1 + # } + + + for(i in 2:(iterations+1)){ + + + tmp <- tryCatch(calibMuso(settings = settings, + parameters = randValues[(i-1),], + silent= TRUE, + skipSpinup = skipSpinup, modifyOriginal=modifyOriginal, postProcString = postProcString), error = function (e) NULL) + if(is.null(tmp)){ + partialResult[,resultRange] <- NA + } else { + partialResult[,resultRange] <- calcLikelihoodsAndRMSE(dataVar=dataVar, + mod=tmp, + mes=measuredData, + likelihoods=likelihood, + alignIndexes=alignIndexes, + musoCodeToIndex = musoCodeToIndex, uncert = uncert) + } + + + partialResult[1:numParameters] <- randValues[(i-1),] + write.table(x=partialResult, file="preservedCalib.csv", append=TRUE, row.names=FALSE, + sep=",", col.names=FALSE) + write.csv(x=tmp, file=paste0(pretag, (i+1),".csv")) + writeLines(as.character(i-1),"progress.txt") + } + + if(threadNumber == 1){ + return(randVals[[1]]) + } + +} + +distributeCores <- function(iterations, numCores){ + perProcess<- iterations %/% numCores + numSimu <- rep(perProcess,numCores) + gainers <- sample(1:numCores, iterations %% numCores) + numSimu[gainers] <- numSimu[gainers] + 1 + numSimu +} + +prepareFromAgroMo <- function(fName){ + obs <- read.table(fName, stringsAsFactors=FALSE, sep = ";", header=T) + obs <- reshape(obs, timevar="var_id", idvar = "date", direction = "wide") + dateCols <- apply(do.call(rbind,(strsplit(obs$date, split = "-"))),2,as.numeric) + colnames(dateCols) <- c("year", "month", "day") + cbind.data.frame(dateCols, obs) +} + + +calcLikelihoodsAndRMSE <- function(dataVar, mod, mes, likelihoods, alignIndexes, musoCodeToIndex, uncert){ + + # NOT COMPATIBLE WITH OLD MEASUREMENT DATA, mes have to be a matrix + likelihoodRMSE <- sapply(names(dataVar),function(key){ + # browser() + modelled <- mod[alignIndexes$mod,musoCodeToIndex[key]] + selected <- grep(sprintf("%s$", key), colnames(mes)) + # browser() + measured <- mes[alignIndexes$meas,selected] + notNA <- sapply(1:nrow(measured), function(x){!any(is.na(measured[x,]))}) + modelled <- modelled[notNA] + measured <- measured[notNA,] + # uncert <- uncert[!is.na(measured)] + + # measured <- measured[!is.na(measured)] + apply(measured, 1, function(x){!any(is.na(x))}) + measured <- t(apply(measured, 1, function(x){if(!any(is.na(x))){x}} )) + if(ncol(measured)!=1){ + m <- measured[,grep("^mean", colnames(measured))] + } + res <- c(likelihoods[[key]](modelled, measured), + sqrt(mean((modelled-m)^2)) + ) + # browser() + res + }) + names(likelihoodRMSE) <- c(sprintf("%s_likelihood",dataVar), sprintf("%s_rmse",dataVar)) + + return(c(likelihoodRMSE[1,],likelihoodRMSE[2,])) +} + +agroLikelihood <- function(modVector,measured){ + mu <- measured[,grep("mean", colnames(measured))] + stdev <- measured[,grep("^sd", colnames(measured))] + ndata <- nrow(measured) + sum(sapply(1:ndata, function(x){ + dnorm(modVector, mu[x], stdev[x], log = TRUE) + }), na.rm=TRUE) +} + + + +maxLikelihoodAgromo <- function (results, imgPath, varName, ...) { + +} diff --git a/RBBGCMuso/R/calibration.R b/RBBGCMuso/R/calibration.R index ce246ef..5c687c3 100644 --- a/RBBGCMuso/R/calibration.R +++ b/RBBGCMuso/R/calibration.R @@ -33,7 +33,7 @@ optiMuso <- function(measuredData, parameters = NULL, startDate = NULL, outVars = NULL, iterations = 30, skipSpinup = TRUE, plotName = "calib.jpg", modifyOriginal=TRUE, likelihood, uncertainity = NULL, - naVal = NULL, postProcString = NULL, w=NULL, lg=FALSE) { + naVal = NULL, postProcString = NULL, w=NULL, lg=FALSE, parallel = TRUE) { # Exanding likelihood likelihoodFull <- as.list(rep(NA,length(dataVar))) names(likelihoodFull) <- names(dataVar) @@ -187,24 +187,24 @@ alignMuso <- function (settings,measuredData) { cbind.data.frame(model=modIndex,meas=measIndex) } -calcLikelihoodsAndRMSE <- function(dataVar, mod, mes, likelihoods, alignIndexes, musoCodeToIndex, uncert){ - - likelihoodRMSE <- sapply(names(dataVar),function(key){ - # browser() - modelled <- mod[alignIndexes$mod,musoCodeToIndex[key]] - measured <- mes[alignIndexes$meas,key] - modelled <- modelled[!is.na(measured)] - # uncert <- uncert[!is.na(measured)] - measured <- measured[!is.na(measured)] - res <- c(likelihoods[[key]](modelled, measured, uncert), - sqrt(mean((modelled-measured)^2)) - ) - res - }) - names(likelihoodRMSE) <- c(sprintf("%s_likelihood",dataVar), sprintf("%s_rmse",dataVar)) - - return(c(likelihoodRMSE[1,],likelihoodRMSE[2,])) -} +# calcLikelihoodsAndRMSE <- function(dataVar, mod, mes, likelihoods, alignIndexes, musoCodeToIndex, uncert){ +# +# likelihoodRMSE <- sapply(names(dataVar),function(key){ +# # browser() +# modelled <- mod[alignIndexes$mod,musoCodeToIndex[key]] +# measured <- mes[alignIndexes$meas,key] +# modelled <- modelled[!is.na(measured)] +# # uncert <- uncert[!is.na(measured)] +# measured <- measured[!is.na(measured)] +# res <- c(likelihoods[[key]](modelled, measured, uncert), +# sqrt(mean((modelled-measured)^2)) +# ) +# res +# }) +# names(likelihoodRMSE) <- c(sprintf("%s_likelihood",dataVar), sprintf("%s_rmse",dataVar)) +# +# return(c(likelihoodRMSE[1,],likelihoodRMSE[2,])) +# } #' musoGlue #' @@ -215,7 +215,12 @@ calcLikelihoodsAndRMSE <- function(dataVar, mod, mes, likelihoods, alignIndexes, #' @export musoGlue <- function(presCalFile, w, delta = 0.17, settings=setupMuso(), parameters=read.csv("parameters.csv", stringsAsFactors=FALSE), lg=FALSE){ - preservedCalib<- read.csv(presCalFile) + if(is.data.frame(presCalFile)){ + preservedCalib <- presCalFile + } else { + preservedCalib <- read.csv(presCalFile) + } + paramIndex <- parameters[(match(colnames(preservedCalib),parameters[,1])),2] paramIndex <- paramIndex[!is.na(paramIndex)] paramIndex <- c(paramIndex, @@ -242,7 +247,8 @@ musoGlue <- function(presCalFile, w, delta = 0.17, settings=setupMuso(), paramet parameterIndexes <- 1:(min(likeIndexes)-1) preservedCalib <- preservedCalib[!is.na(preservedCalib$combined),] unfilteredLikelihood <- preservedCalib$combined - preservedCalibtop5 <- preservedCalib[preservedCalib$combined>quantile(preservedCalib$combined,0.95),] + top5points <- preservedCalib$combined>quantile(preservedCalib$combined,0.95) + preservedCalibtop5 <- preservedCalib[,] optRanges <-t(apply(preservedCalibtop5,2,function(x) quantile(x,c(0.05,0.5,0.95)))) pdf("dotplot.pdf") if(lg){ @@ -261,20 +267,40 @@ musoGlue <- function(presCalFile, w, delta = 0.17, settings=setupMuso(), paramet abline(v=optRanges[i,3],col="red") } + par(pari) dev.off() maxParValues <- preservedCalibtop5[which.max(preservedCalibtop5$combined),] maxParIndexes <- paramIndex write.csv(cbind.data.frame(calibrationPar=maxParValues,parameters=maxParIndexes),"maxLikelihood.csv") write.csv(optRanges,"optRanges.csv") - optInterval <-t(apply(preservedCalibtop5,2,function(x) quantile(x,c(0.5-delta,0.5+delta)))) - optParamRange <- cbind.data.frame(rownames(optInterval)[parameterIndexes],as.numeric(paramIndex),optInterval[parameterIndexes,]) - optimalEpc <- musoRand(optParamRange,iterations = 2) - optimalEpc[[2]] <- optimalEpc[[2]][1,] - write.csv(as.data.frame(optimalEpc),"epcOptim.csv") - print(head(optRanges,n=-2)) - calibMuso(calibrationPar=optimalEpc[[1]],parameters=optimalEpc[[2]]) - file.copy(settings$epcInput[2],"epcOptim.epc") + # browser() + # There are some serious problems with this implementation. The uncertainity bouns are not for the parameters, but for the output values. The median is pointwise median for all simulation. + # And the 95 and 5 percentile also. + # dataVec <- preservedCalibtop5$combined + # closestToMedian <- function (dataVec) { + # match(sort(dataVec)[min(which(sort(dataVec)>=median(dataVec)))], dataVec) + # } + # + # while(is.null(optimalEpc)){ + # match(quantile(preservedCalibtop5$combined,0.5), preservedCalibtop5$combined) + # optInterval <-t(apply(preservedCalibtop5,2,function(x) quantile(x,c(0.5-delta,0.5+delta)))) + # optParamRange <- cbind.data.frame(rownames(optInterval)[parameterIndexes],as.numeric(paramIndex),optInterval[parameterIndexes,]) + # optimalEpc <- tryCatch(musoRand(optParamRange,iterations = 2), error=function(e){NULL}) + # delta <- delta*1.05 + # if(delta > 0.5){ + # delta <- 0.5 + # } + # if((delta == 0.5) && is.null(optimalEpc)){ + # stop("cannot find optimal value in the given range") + # } + # } + # print("getOptim") + # optimalEpc[[2]] <- optimalEpc[[2]][1,] + # write.csv(as.data.frame(optimalEpc),"epcOptim.csv") + # print(head(optRanges,n=-2)) + # calibMuso(calibrationPar=optimalEpc[[1]],parameters=optimalEpc[[2]]) + # file.copy(settings$epcInput[2],"epcOptim.epc") } generateOptEpc <- function(optRanges,delta, maxLikelihood=FALSE){ diff --git a/RBBGCMuso/R/changeMuso.R b/RBBGCMuso/R/changeMuso.R index 97b70e7..9a34a0c 100644 --- a/RBBGCMuso/R/changeMuso.R +++ b/RBBGCMuso/R/changeMuso.R @@ -1,14 +1,11 @@ -#' This is the function which is capable to change multiple specific lines to others using their row numbers. +#' changemulline #' #' The function uses the previous changspecline function to operate. - ##From now changespecline is in the forarcheologist file, because itis no longer needed #' #' @author Roland Hollos -#' @keywords internal -#' +#' @export - -changemulline <- function(filePaths, calibrationPar, contents, src){ +changemulline <- function(filePaths, calibrationPar, contents, src, outFiles=filePaths){ # browser() if(is.null(src)){ src <- filePaths @@ -19,7 +16,7 @@ changemulline <- function(filePaths, calibrationPar, contents, src){ fileStringVector <<- changeByIndex(index, content, fileStringVector) }, calibrationPar, contents) - writeLines(fileStringVector, filePaths) + writeLines(fileStringVector, outFiles) } diff --git a/RBBGCMuso/R/setupMuso.R b/RBBGCMuso/R/setupMuso.R index 8873478..ad7d2f5 100644 --- a/RBBGCMuso/R/setupMuso.R +++ b/RBBGCMuso/R/setupMuso.R @@ -68,7 +68,7 @@ setupMuso <- function(executable=NULL, # inputParser <- function(string,fileName,counter,value=TRUE){ - unlist(strsplit(grep(string,fileName,value=TRUE),"[\ \t]"))[counter] + unlist(strsplit(grep(string,fileName,value=TRUE, perl = TRUE),"[\ \t]", useBytes = TRUE))[counter] } # outMaker <- function(inputVar,grepString,filep){ @@ -104,8 +104,8 @@ setupMuso <- function(executable=NULL, #iniChangedp <- FALSE if(is.null(iniInput)){ - spinups<-grep("s.ini$",list.files(inputLoc),value=TRUE) - normals<-grep("n.ini$",list.files(inputLoc),value=TRUE) + spinups<-grep("s.ini$",list.files(inputLoc),value=TRUE, perl = TRUE) + normals<-grep("n.ini$",list.files(inputLoc),value=TRUE, perl = TRUE) if(length(spinups)==1){ iniInput[1] <- file.path(inputLoc,spinups) @@ -124,8 +124,8 @@ setupMuso <- function(executable=NULL, ##read the ini files for the further changes iniFiles<-lapply(iniInput, function (x) readLines(x,-1)) - iniFiles[[1]] <- gsub("\\","/", iniFiles[[1]],fixed=TRUE) #replacing \ to / - iniFiles[[2]] <- gsub("\\","/", iniFiles[[2]],fixed=TRUE) #replacing \ to / + iniFiles[[1]] <- gsub("\\\\","/", iniFiles[[1]], perl = TRUE) #replacing \ to / + iniFiles[[2]] <- gsub("\\\\","/", iniFiles[[2]], perl = TRUE) #replacing \ to / names(iniFiles) <- c("spinup","normal") @@ -156,21 +156,21 @@ setupMuso <- function(executable=NULL, # if(is.null(mapData)){ # - outIndex<-grep("DAILY_OUTPUT",iniFiles[[2]])+1 - numVar<-as.numeric(unlist(strsplit(iniFiles[[2]][outIndex],"[\ \t]"))[1]) + outIndex<-grep("DAILY_OUTPUT",iniFiles[[2]], perl = TRUE)+1 + numVar<-as.numeric(unlist(strsplit(iniFiles[[2]][outIndex],"[\ \t]", useBytes = TRUE))[1]) dailyVarCodes<-tryCatch(iniFiles[[2]][(outIndex+1):(outIndex+numVar)], error = function(e){ stop("Cannot read indexes of output variables from the normal ini file, please make sure you have not skiped a line after the flag: \"DAILY_OUTPUT\"") }) - dailyVarCodes<-unlist(lapply(dailyVarCodes, function(x) unlist(strsplit(x,"[\ \t]"))[1])) - dailyVarnames<-unlist(lapply(dailyVarCodes, function(x) musoMapping(unlist(strsplit(x,"[\ \t]"))[1]))) + dailyVarCodes<-unlist(lapply(dailyVarCodes, function(x) unlist(strsplit(x,"[\ \t]", useBytes = TRUE))[1])) + dailyVarnames<-unlist(lapply(dailyVarCodes, function(x) musoMapping(unlist(strsplit(x,"[\ \t]", useBytes = TRUE))[1]))) - outIndex<-grep("ANNUAL_OUTPUT",iniFiles[[2]])+1 - numVar<-as.numeric(unlist(strsplit(iniFiles[[2]][outIndex],"[\ \t]"))[1]) + outIndex<-grep("ANNUAL_OUTPUT",iniFiles[[2]], perl = TRUE)+1 + numVar<-as.numeric(unlist(strsplit(iniFiles[[2]][outIndex],"[\ \t]", useBytes = TRUE))[1]) annualVarCodes<-iniFiles[[2]][(outIndex+1):(outIndex+numVar)] - annualVarCodes<-unlist(lapply(annualVarCodes, function(x) unlist(strsplit(x,"[\ \t]"))[1])) - annualVarnames<-unlist(lapply(annualVarCodes, function(x) musoMapping(unlist(strsplit(x,"[\ \t]"))[1]))) + annualVarCodes<-unlist(lapply(annualVarCodes, function(x) unlist(strsplit(x,"[\ \t]", useBytes = TRUE))[1])) + annualVarnames<-unlist(lapply(annualVarCodes, function(x) musoMapping(unlist(strsplit(x,"[\ \t]", useBytes = TRUE))[1]))) outputVars<-list(dailyVarnames,annualVarnames) # browser() # } else { @@ -206,8 +206,8 @@ setupMuso <- function(executable=NULL, } outputName <- character(2) - outputName[1] <- basename(unlist(strsplit(iniFiles[[1]][grep("OUTPUT_CONTROL",iniFiles[[1]])+1],"[\ \t]"))[1]) - outputName[2] <- basename(unlist(strsplit(iniFiles[[2]][grep("OUTPUT_CONTROL",iniFiles[[2]])+1],"[\ \t]"))[1]) + outputName[1] <- basename(unlist(strsplit(iniFiles[[1]][grep("OUTPUT_CONTROL",iniFiles[[1]], perl = TRUE)+1],"[\ \t]", useBytes = TRUE))[1]) + outputName[2] <- basename(unlist(strsplit(iniFiles[[2]][grep("OUTPUT_CONTROL",iniFiles[[2]], perl = TRUE)+1],"[\ \t]", useBytes = TRUE))[1]) ## outputName <- unlist(strsplit(grep("output",grep("prefix",iniFiles[[2]],value=TRUE),value=TRUE),"[\ \t]"))[1] ##THIS IS AN UGLY SOLUTION, WHICH NEEDS AN UPGRADE!!! FiXED (2017.09.11) ## outputName <- unlist(strsplit(grep("prefix for output files",iniFiles[[2]],value=TRUE),"[\ \t]"))[1] @@ -220,7 +220,7 @@ setupMuso <- function(executable=NULL, if(is.null(outputLoc)){ ## outputLoc<-paste((rev(rev(unlist(strsplit(outputName,"/")))[-1])),collapse="/") - outputLoc <- dirname(unlist(strsplit(iniFiles[[2]][grep("OUTPUT_CONTROL",iniFiles[[2]])+1],"[\ \t]"))[1]) + outputLoc <- dirname(unlist(strsplit(iniFiles[[2]][grep("OUTPUT_CONTROL",iniFiles[[2]], perl = TRUE)+1],"[\ \t]", useBytes = TRUE))[1]) if(substr(outputLoc,start = 1,stop = 1)!="/"){ ##if the outputName is not absolute path make it absolute outputLoc <- file.path(inputLoc,outputLoc) @@ -233,12 +233,12 @@ setupMuso <- function(executable=NULL, inputFiles<-c(iniInput,epcInput,metInput) numData<-rep(NA,3) - numYears <- as.numeric(unlist(strsplit(iniFiles[[2]][grep("TIME_DEFINE",iniFiles[[2]])+1],"[\ \t]"))[1]) + numYears <- as.numeric(unlist(strsplit(iniFiles[[2]][grep("TIME_DEFINE",iniFiles[[2]], perl = TRUE)+1],"[\ \t]", useBytes = TRUE))[1]) ## numYears<-unlist(read.table(iniInput[2],skip = 14,nrows = 1)[1]) - numValues <- as.numeric(unlist(strsplit(iniFiles[[2]][grep("DAILY_OUTPUT",iniFiles[[2]])+1],"[\ \t]"))[1]) + numValues <- as.numeric(unlist(strsplit(iniFiles[[2]][grep("DAILY_OUTPUT",iniFiles[[2]], perl = TRUE)+1],"[\ \t]", useBytes = TRUE))[1]) ## numValues will be replaced to numVar ## numValues<-unlist(read.table(iniInput[2],skip=102,nrows = 1)[1]) - startYear <- as.numeric(unlist(strsplit(iniFiles[[2]][grep("TIME_DEFINE",iniFiles[[2]])+2],"[\ \t]"))[1]) + startYear <- as.numeric(unlist(strsplit(iniFiles[[2]][grep("TIME_DEFINE",iniFiles[[2]], perl = TRUE)+2],"[\ \t]", useBytes = TRUE))[1]) numData[1] <- numValues * numYears * 365 # Have to corrigate leapyears numData[2] <- numYears * numValues*12 @@ -263,9 +263,9 @@ setupMuso <- function(executable=NULL, searchBellow <- function(inFile, key, stringP = TRUE, n=1, management = FALSE){ if(stringP){ - unlist(strsplit(inFile[grep(key,inFile,perl=TRUE)+n],split = "\\s+"))[1] + unlist(strsplit(inFile[grep(key,inFile, perl=TRUE)+n],split = "\\s+", useBytes = TRUE))[1] } else { - as.numeric(unlist(strsplit(inFile[grep(key,inFile,perl=TRUE)+n],split = "\\s+"))[1]) + as.numeric(unlist(strsplit(inFile[grep(key,inFile,perl=TRUE)+n],split = "\\s+", useBytes = TRUE))[1]) } } soilFile <- NULL diff --git a/RBBGCMuso/R/thread-1 b/RBBGCMuso/R/thread-1 new file mode 100644 index 0000000..29d6383 --- /dev/null +++ b/RBBGCMuso/R/thread-1 @@ -0,0 +1 @@ +100 diff --git a/RBBGCMuso/R/thread-2 b/RBBGCMuso/R/thread-2 new file mode 100644 index 0000000..29d6383 --- /dev/null +++ b/RBBGCMuso/R/thread-2 @@ -0,0 +1 @@ +100 diff --git a/RBBGCMuso/R/thread-3 b/RBBGCMuso/R/thread-3 new file mode 100644 index 0000000..29d6383 --- /dev/null +++ b/RBBGCMuso/R/thread-3 @@ -0,0 +1 @@ +100 diff --git a/RBBGCMuso/R/thread-4 b/RBBGCMuso/R/thread-4 new file mode 100644 index 0000000..29d6383 --- /dev/null +++ b/RBBGCMuso/R/thread-4 @@ -0,0 +1 @@ +100 diff --git a/RBBGCMuso/man/calibrateMuso.Rd b/RBBGCMuso/man/calibrateMuso.Rd new file mode 100644 index 0000000..5710638 --- /dev/null +++ b/RBBGCMuso/man/calibrateMuso.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calibrateMuso.R +\name{calibrateMuso} +\alias{calibrateMuso} +\title{calibrateMuso} +\usage{ +calibrateMuso( + measuredData, + parameters = NULL, + startDate = NULL, + endDate = NULL, + formatString = "\%Y-\%m-\%d", + dataVar, + outLoc = "./calib", + preTag = "cal-", + settings = setupMuso(), + outVars = NULL, + iterations = 100, + skipSpinup = TRUE, + plotName = "calib.jpg", + modifyOriginal = TRUE, + likelihood, + uncertainity = NULL, + naVal = NULL, + postProcString = NULL, + thread_prefix = "thread", + numCores = (parallel::detectCores() - 1), + pb = txtProgressBar(min = 0, max = iterations, style = 3), + maxLikelihoodEpc = TRUE, + pbUpdate = setTxtProgressBar +) +} +\description{ +This funtion uses the Monte Carlo technique to uniformly sample the parameter space from user defined parameters of the Biome-BGCMuSo model. The sampling algorithm ensures that the parameters are constrained by the model logic which means that parameter dependencies are fully taken into account (parameter dependency means that e.g leaf C:N ratio must be smaller than C:N ratio of litter; more complicated rules apply to the allocation parameters where the allocation fractions to different plant compartments must sum up 1). This function implements a mathematically correct solution to provide uniform distriution for all selected parameters. +} +\author{ +Roland HOLLOS +} diff --git a/RBBGCMuso/man/changemulline.Rd b/RBBGCMuso/man/changemulline.Rd index e7ec966..848e443 100644 --- a/RBBGCMuso/man/changemulline.Rd +++ b/RBBGCMuso/man/changemulline.Rd @@ -2,16 +2,9 @@ % Please edit documentation in R/changeMuso.R \name{changemulline} \alias{changemulline} -\title{This is the function which is capable to change multiple specific lines to others using their row numbers.} +\title{changemulline} \usage{ -changemulline( - filePaths, - calibrationPar, - contents, - fileOut, - fileToChange, - modifyOriginal = FALSE -) +changemulline(filePaths, calibrationPar, contents, src) } \description{ The function uses the previous changspecline function to operate. @@ -19,4 +12,3 @@ The function uses the previous changspecline function to operate. \author{ Roland Hollos } -\keyword{internal} diff --git a/RBBGCMuso/man/musoGlue.Rd b/RBBGCMuso/man/musoGlue.Rd index b4960e7..c3b878c 100644 --- a/RBBGCMuso/man/musoGlue.Rd +++ b/RBBGCMuso/man/musoGlue.Rd @@ -4,7 +4,14 @@ \alias{musoGlue} \title{musoGlue} \usage{ -musoGlue(preservedCalib, w) +musoGlue( + presCalFile, + w, + delta = 0.17, + settings = setupMuso(), + parameters = read.csv("parameters.csv", stringsAsFactors = FALSE), + lg = FALSE +) } \arguments{ \item{plotName}{u} diff --git a/RBBGCMuso/man/optiMuso.Rd b/RBBGCMuso/man/optiMuso.Rd index d338cbc..6543dc9 100644 --- a/RBBGCMuso/man/optiMuso.Rd +++ b/RBBGCMuso/man/optiMuso.Rd @@ -20,9 +20,12 @@ optiMuso( plotName = "calib.jpg", modifyOriginal = TRUE, likelihood, + uncertainity = NULL, naVal = NULL, postProcString = NULL, - w = NULL + w = NULL, + lg = FALSE, + parallel = TRUE ) } \arguments{ diff --git a/calibrateMuso.R b/calibrateMuso.R new file mode 100644 index 0000000..46c46e8 --- /dev/null +++ b/calibrateMuso.R @@ -0,0 +1,295 @@ +#' calibrateMuso +#' +#' This funtion uses the Monte Carlo technique to uniformly sample the parameter space from user defined parameters of the Biome-BGCMuSo model. The sampling algorithm ensures that the parameters are constrained by the model logic which means that parameter dependencies are fully taken into account (parameter dependency means that e.g leaf C:N ratio must be smaller than C:N ratio of litter; more complicated rules apply to the allocation parameters where the allocation fractions to different plant compartments must sum up 1). This function implements a mathematically correct solution to provide uniform distriution for all selected parameters. +#' @author Roland HOLLOS +#' @export +calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, + endDate = NULL, formatString = "%Y-%m-%d", + dataVar, outLoc = "./calib", + preTag = "cal-", settings = setupMuso(), + outVars = NULL, iterations = 30, + skipSpinup = TRUE, plotName = "calib.jpg", + modifyOriginal=TRUE, likelihood, uncertainity = NULL, + naVal = NULL, postProcString = NULL, + tread_prefix="thread", numcores = (parallel::detectCores()-1), pb = txtProgressBar(min=0, max=iterations, style=3), + pbUpdate = setTxtProgressBar){ + + # ____ _ _ _ _ + # / ___|_ __ ___ __ _| |_ ___ | |_| |__ _ __ ___ __ _ __| |___ + # | | | '__/ _ \/ _` | __/ _ \ | __| '_ \| '__/ _ \/ _` |/ _` / __| + # | |___| | | __/ (_| | || __/ | |_| | | | | | __/ (_| | (_| \__ \ + # \____|_| \___|\__,_|\__\___| \__|_| |_|_| \___|\__,_|\__,_|___/ + + + + copyToThreadDirs(thread_prefix, numcores = numcores, runDir = settings$inputLoc) + + # ____ _ _ _ + # | _ \ _ _ _ __ | |_| |__ _ __ ___ __ _ __| |___ + # | |_) | | | | '_ \ | __| '_ \| '__/ _ \/ _` |/ _` / __| + # | _ <| |_| | | | | | |_| | | | | | __/ (_| | (_| \__ \ + # |_| \_\\__,_|_| |_| \__|_| |_|_| \___|\__,_|\__,_|___/ + + threadCount <- distributeCores(iterations, numCores) + + fut <- lapply(1:numCores, function(i) { + # future({ + musoSingleThread(measuredData, parameters, startDate, + endDate, formatString, + dataVar, outLoc, + preTag, settings, + outVars, iterations = threadCount[i], + skipSpinup, plotName, + modifyOriginal, likelihood, uncertainity, + naVal, postProcString, i) + # }) + }) + + # __ ___ _ _ + # \ \ / / |__ __ _| |_ ___| |__ _ __ _ __ ___ ___ ___ ___ ___ + # \ \ /\ / /| '_ \ / _` | __/ __| '_ \ | '_ \| '__/ _ \ / __/ _ \/ __/ __| + # \ V V / | | | | (_| | || (__| | | | | |_) | | | (_) | (_| __/\__ \__ \ + # \_/\_/ |_| |_|\__,_|\__\___|_| |_| | .__/|_| \___/ \___\___||___/___/ + # |_| + + getProgress <- function(){ + threadfiles <- list.files(settings$inputLoc, pattern="progress.txt", recursive = TRUE) + if(length(threadfiles)==0){ + return(0) + } else { + sum(sapply(threadfiles, function(x){ + partRes <- readLines(x) + if(length(partRes)==0){ + return(0) + } else { + return(as.numeric(partRes)) + } + + })) + + } + } + + progress <- 0 + while(progress < 400){ + Sys.sleep(1) + progress <- getProgress() + pbUpdate(pb,as.numeric(progress)) + } + close(pb) + + # ____ _ _ + # / ___|___ _ __ ___ | |__ (_)_ __ ___ + # | | / _ \| '_ ` _ \| '_ \| | '_ \ / _ \ + # | |__| (_) | | | | | | |_) | | | | | __/ + # \____\___/|_| |_| |_|_.__/|_|_| |_|\___| + + + # ____ _ _ _ _____ + # / ___| | | | | | ____| + # | | _| | | | | | _| + # | |_| | |__| |_| | |___ + # \____|_____\___/|_____| + + + # musoGlue("preservedCalib.csv",w=w, lg = lg) +} + +copyToThreadDirs <- function(prefix="thread", numcores=parallel::detectCores()-1, runDir="."){ + dir.create(file.path(runDir,prefix), showWarnings=TRUE) + fileNames <- grep("^thread.*", list.files(runDir), value=TRUE, invert=TRUE) + invisible(sapply(1:numcores,function(corenum){ + threadDir <- file.path(runDir,prefix,paste0(prefix,"_",corenum)) + dir.create(threadDir, showWarnings=FALSE) + file.copy(from=fileNames,to=threadDir, overwrite=FALSE) + })) +} + +musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, + endDate = NULL, formatString = "%Y-%m-%d", + dataVar, outLoc = "./calib", + preTag = "cal-", settings = setupMuso(), + outVars = NULL, iterations = 30, + skipSpinup = TRUE, plotName = "calib.jpg", + modifyOriginal=TRUE, likelihood, uncertainity = NULL, + naVal = NULL, postProcString = NULL, threadNumber) { + + setwd(paste0("thread/thread-",threadNumber)) + # Exanding likelihood + likelihoodFull <- as.list(rep(NA,length(dataVar))) + names(likelihoodFull) <- names(dataVar) + if(!missing(likelihood)) { + lapply(names(likelihood),function(x){ + likelihoodFull[[x]] <<- likelihood[[x]] + }) + } + defaultLikelihood <- which(is.na(likelihood)) + if(length(defaultLikelihood)>0){ + likelihoodFull[[defaultLikelihood]] <- (function(x, y){ + exp(-sqrt(mean((x-y)^2))) + }) + } + + mdata <- measuredData + if(is.null(parameters)){ + parameters <- tryCatch(read.csv("parameters.csv", stringsAsFactor=FALSE), error = function (e) { + stop("You need to specify a path for the parameters.csv, or a matrix.") + }) + } else { + if((!is.list(parameters)) & (!is.matrix(parameters))){ + parameters <- tryCatch(read.csv(parameters, stringsAsFactor=FALSE), error = function (e){ + stop("Cannot find neither parameters file neither the parameters matrix") + }) + }} + + outLoc <- normalizePath(outLoc) + outLocPlain <- basename(outLoc) + currDir <- getwd() + + if(!dir.exists(outLoc)){ + dir.create(outLoc) + warning(paste(outLoc," is not exists, so it was created")) + } + + outLoc <- normalizePath(outLoc) + parameterNames <- parameters[,1] + pretag <- file.path(outLoc,preTag) + ##reading the original epc file at the specified + ## row numbers + print("optiMuso is randomizing the epc parameters now...",quote = FALSE) + if(iterations < 3000){ + randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = 3000) + randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),] + } else { + randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations) + } + + origEpc <- readValuesFromFile(settings$epc[2],randVals[[1]]) + partialResult <- matrix(ncol=length(randVals[[1]])+2*length(dataVar)) + colN <- randVals[[1]] + colN[match(parameters[,2],randVals[[1]])] <- parameters[,1] + colnames(partialResult) <- c(colN,sprintf("%s_likelihood",names(dataVar)), + sprintf("%s_rmse",names(dataVar))) + numParameters <- length(colN) + partialResult[1:numParameters] <- origEpc + ## Prepare the preservedCalib matrix for the faster + ## run. + + pretag <- file.path(outLoc,preTag) + + musoCodeToIndex <- sapply(dataVar,function(musoCode){ + settings$dailyOutputTable[settings$dailyOutputTable$code == musoCode,"index"] + }) + resultRange <- (numParameters + 1):(ncol(partialResult)) + ## Creating function for generating separate + ## csv files for each run + + progBar <- txtProgressBar(1,iterations,style=3) + settings$iniInput[2] %>% + (function(x) paste0(dirname(x),"/",tools::file_path_sans_ext(basename(x)),"-tmp.",tools::file_ext(x))) %>% + unlink + randValues <- randVals[[2]] + + settings$calibrationPar <- randVals[[1]] + + if(!is.null(naVal)){ + measuredData <- as.data.frame(measuredData) + measuredData[measuredData == naVal] <- NA + } + + alignIndexes <- alignMuso(settings,measuredData) + if(!is.null(uncertainity)){ + uncert <- measuredData[alignIndexes$meas,uncertainity] + } else { + uncert <- NULL + } + # browser() + origModellOut <- calibMuso(settings=settings, silent=TRUE, skipSpinup = skipSpinup, postProcString=postProcString, modifyOriginal=modifyOriginal) + partialResult[,resultRange] <- calcLikelihoodsAndRMSE(dataVar=dataVar, + mod=origModellOut, + mes=measuredData, + likelihoods=likelihood, + alignIndexes=alignIndexes, + musoCodeToIndex = musoCodeToIndex,uncert=uncert) + write.csv(x=origModellOut, file=paste0(pretag, 1, ".csv")) + print("Running the model with the random epc values...", quote = FALSE) + + # if(!is.null(postProcString)){ + # colNumb <- length(settings$dailyVarCodes) + 1 + # } + + write.csv(x=partialResult, file="preservedCalib.csv",row.names=FALSE) + for(i in 2:(iterations+1)){ + # browser() + tmp <- tryCatch(calibMuso(settings = settings, + parameters = randValues[(i-1),], + silent= TRUE, + skipSpinup = skipSpinup, modifyOriginal=modifyOriginal, postProcString = postProcString), error = function (e) NULL) + if(is.null(tmp)){ + partialResult[,resultRange] <- NA + } else { + partialResult[,resultRange] <- calcLikelihoodsAndRMSE(dataVar=dataVar, + mod=tmp, + mes=measuredData, + likelihoods=likelihood, + alignIndexes=alignIndexes, + musoCodeToIndex = musoCodeToIndex, uncert = uncert) + } + + + partialResult[1:numParameters] <- randValues[(i-1),] + write.table(x=partialResult, file="preservedCalib.csv", append=TRUE, row.names=FALSE, + sep=",", col.names=FALSE) + write.csv(x=tmp, file=paste0(pretag, (i+1),".csv")) + writeLines(as.character(i),"progress.txt") + } + +} + +distributeCores <- function(iterations, numCores){ + perProcess<- iterations %/% numCores + numSimu <- rep(perProcess,numCores) + gainers <- sample(1:numCores, iterations %% numCores) + numSimu[gainers] <- numSimu[gainers] + 1 + numSimu +} + +prepareFromAgroMo <- function(fName){ + obs <- read.table(fName, stringsAsFactors=FALSE, sep = ";", header=T) + obs <- reshape(obs, timevar="var_id", idvar = "date", direction = "wide") + dateCols <- apply(do.call(rbind,(strsplit(obs$date, split = "-"))),2,as.numeric) + colnames(dateCols) <- c("year", "month", "day") + cbind.data.frame(dateCols, obs) +} + + +calcLikelihoodsAndRMSE <- function(dataVar, mod, mes, likelihoods, alignIndexes, musoCodeToIndex, uncert){ + + likelihoodRMSE <- sapply(names(dataVar),function(key){ + # browser() + modelled <- mod[alignIndexes$mod,musoCodeToIndex[key]] + selected <- grep(sprintf("%s$", key), colnames(mes)) + measured <- mes[alignIndexes$meas,selected] + modelled <- modelled[!is.na(measured)] + # uncert <- uncert[!is.na(measured)] + measured <- measured[!is.na(measured)] + res <- c(likelihoods[[key]](modelled, measured, uncert), + sqrt(mean((modelled-measured)^2)) + ) + res + }) + names(likelihoodRMSE) <- c(sprintf("%s_likelihood",dataVar), sprintf("%s_rmse",dataVar)) + + return(c(likelihoodRMSE[1,],likelihoodRMSE[2,])) +} + +agroLikelihood <- function(modVector,measured){ + mu <- measured[,grep("^mean", colnames(measured))] + stdev <- measured[,grep("^sd", colnames(measured))] + ndata <- nrow(measured) + sum(sapply(1:ndata, function(x){ + dnorm(modVector, mu[x], stdev[x], log = TRUE) + })) +} + +# prepareFromAgroMo("/home/hollorol/agromo/calibration/martonvasar/MV_highN.obs") diff --git a/parallelRun.R b/parallelRun.R new file mode 100644 index 0000000..05da565 --- /dev/null +++ b/parallelRun.R @@ -0,0 +1,90 @@ +library(parallel) +library('future') +plan(multiprocess) +library('RBBGCMuso') + +a <- tempdir() +setwd(a) +file.copy(from="~/R/x86_64-pc-linux-gnu-library/3.5/RBBGCMuso/examples/hhs",to="./", recursive = TRUE) +setwd("hhs") +list.files() + +settings <- setupMuso() +setupMuso() +settings$outputLoc + + +copyToThreadDirs <- function(prefix="thread", numcores=parallel::detectCores()-1, runDir="."){ + dir.create(file.path(runDir,prefix), showWarnings=TRUE) + fileNames <- grep("^thread.*", list.files(runDir), value=TRUE, invert=TRUE) + invisible(sapply(1:numcores,function(corenum){ + threadDir <- file.path(runDir,prefix,paste0(prefix,"_",corenum)) + dir.create(threadDir, showWarnings=FALSE) + file.copy(from=fileNames,to=threadDir, overwrite=FALSE) + })) +} + +copyToThreadDirs() +unlink("thread", recursive=TRUE) + +procFun <- function(index){ + progressState <- tempfile(pattern=paste("thread",index,sep="-", tmpdir="./")) + for(i in 1:100){ + Sys.sleep(1) + writeLines(as.character(i),paste("thread",index,sep="-")) + } +} + +futu <- vector(mode="list", length=4) +names(futu) <- 1:4 +futu + + +getProgress <- function(){ + threadfiles <- list.files(pattern="thread*") + if(length(threadfiles)==0){ + return(0) + } else { + sum(sapply(threadfiles, function(x){ + partRes <- readLines(x) + if(length(partRes)==0){ + return(0) + } else { + return(as.numeric(partRes)) + } + + })) + + } +} + +getProgress() +futu +list.files() +readLines("threadi-1") +procFun(8) +file.remove(pattern="thread*") +file.remove((list.files(pattern="thread*"))) + + +wachProgress <- function(){ + progress <- 0 + while(progress < 400){ + Sys.sleep(1) + progress <- getProgress() + print(paste(as.numeric(progress)/400*100,"%")) + } +} + +for(i in 1:4){ + futu[[as.character(i)]] <- future({procFun(i)}) +} +lapply(1:4,function(i) future({procFun(i)})) +pb <- txtProgressBar(min=0,max=400,style=3) +progress <- 0 +while(progress < 400){ + Sys.sleep(1) + progress <- getProgress() + setTxtProgressBar(pb,as.numeric(progress)) +} +close(pb) From 3a973d0172ee41e165271fd3f4fbd763c3680c6b Mon Sep 17 00:00:00 2001 From: Hollos Roland Date: Thu, 24 Sep 2020 17:36:06 +0200 Subject: [PATCH 2/7] some bugfixing --- RBBGCMuso/R/calibrateMuso.R | 20 +++++++++++++------- calibrateMuso.R | 6 +++--- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/RBBGCMuso/R/calibrateMuso.R b/RBBGCMuso/R/calibrateMuso.R index 0deefa5..2ac1992 100644 --- a/RBBGCMuso/R/calibrateMuso.R +++ b/RBBGCMuso/R/calibrateMuso.R @@ -16,10 +16,10 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, maxLikelihoodEpc=TRUE, pbUpdate = setTxtProgressBar, method="GLUE",lg = FALSE, w=NULL, ...){ - - file.remove(list.files(path = settings$inputLoc, pattern="progress.txt", recursive = TRUE)) - file.remove(list.files(path = settings$inputLoc, pattern="preservedCalib.csv", recursive = TRUE)) - unlink("thread",recursive=TRUE) + future::plan(future::multisession) + file.remove(list.files(path = settings$inputLoc, pattern="progress.txt", recursive = TRUE, full.names=TRUE)) + file.remove(list.files(path = settings$inputLoc, pattern="preservedCalib.csv", recursive = TRUE, full.names=TRUE)) + unlink(file.path(settings$inputLoc,"thread"),recursive=TRUE) # ____ _ _ _ _ # / ___|_ __ ___ __ _| |_ ___ | |_| |__ _ __ ___ __ _ __| |___ @@ -94,9 +94,15 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, while(progress < iterations){ Sys.sleep(1) progress <- tryCatch(getProgress(), error=function(e){progress}) - pbUpdate(pb,as.numeric(progress)) + if(is.null(pb)){ + pbUpdate(as.numeric(progress)) + } else { + pbUpdate(pb,as.numeric(progress)) + } + } + if(!is.null(pb)){ + close(pb) } - close(pb) # ____ _ _ # / ___|___ _ __ ___ | |__ (_)_ __ ___ @@ -164,7 +170,7 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, copyToThreadDirs <- function(prefix="thread", numcores=parallel::detectCores()-1, runDir="."){ dir.create(file.path(runDir,prefix), showWarnings=TRUE) - fileNames <- grep("^thread.*", list.files(runDir), value=TRUE, invert=TRUE) + fileNames <- grep("^thread.*", list.files(runDir,full.names=TRUE), value=TRUE, invert=TRUE) invisible(sapply(1:numcores,function(corenum){ threadDir <- file.path(runDir,prefix,paste0(prefix,"_",corenum)) dir.create(threadDir, showWarnings=FALSE) diff --git a/calibrateMuso.R b/calibrateMuso.R index 46c46e8..8df08d0 100644 --- a/calibrateMuso.R +++ b/calibrateMuso.R @@ -21,7 +21,6 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, # \____|_| \___|\__,_|\__\___| \__|_| |_|_| \___|\__,_|\__,_|___/ - copyToThreadDirs(thread_prefix, numcores = numcores, runDir = settings$inputLoc) # ____ _ _ _ @@ -33,7 +32,7 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, threadCount <- distributeCores(iterations, numCores) fut <- lapply(1:numCores, function(i) { - # future({ + future({ musoSingleThread(measuredData, parameters, startDate, endDate, formatString, dataVar, outLoc, @@ -42,7 +41,7 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, skipSpinup, plotName, modifyOriginal, likelihood, uncertainity, naVal, postProcString, i) - # }) + }) }) # __ ___ _ _ @@ -96,6 +95,7 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, } copyToThreadDirs <- function(prefix="thread", numcores=parallel::detectCores()-1, runDir="."){ + browser() dir.create(file.path(runDir,prefix), showWarnings=TRUE) fileNames <- grep("^thread.*", list.files(runDir), value=TRUE, invert=TRUE) invisible(sapply(1:numcores,function(corenum){ From 2d3e9a1660b40df3edb00fa6bfbd48ccfa50fda6 Mon Sep 17 00:00:00 2001 From: Hollos Roland Date: Thu, 24 Sep 2020 17:41:23 +0200 Subject: [PATCH 3/7] temporarly disable ecmwf --- RBBGCMuso/DESCRIPTION | 1 - 1 file changed, 1 deletion(-) diff --git a/RBBGCMuso/DESCRIPTION b/RBBGCMuso/DESCRIPTION index c2d0d5b..c267c32 100644 --- a/RBBGCMuso/DESCRIPTION +++ b/RBBGCMuso/DESCRIPTION @@ -28,7 +28,6 @@ Imports: data.table, gridExtra, lubridate, - ecmwfr, openxlsx, ncdf4, tcltk From ac118f6e6572781b971a7f60551645b68f82468e Mon Sep 17 00:00:00 2001 From: Hollos Roland Date: Thu, 24 Sep 2020 17:59:39 +0200 Subject: [PATCH 4/7] delete getMeteoData1.R temporally --- RBBGCMuso/R/getMeteoData1.R | 533 ------------------------------------ RBBGCMuso/R/thread-1 | 1 - RBBGCMuso/R/thread-2 | 1 - RBBGCMuso/R/thread-3 | 1 - RBBGCMuso/R/thread-4 | 1 - 5 files changed, 537 deletions(-) delete mode 100644 RBBGCMuso/R/getMeteoData1.R delete mode 100644 RBBGCMuso/R/thread-1 delete mode 100644 RBBGCMuso/R/thread-2 delete mode 100644 RBBGCMuso/R/thread-3 delete mode 100644 RBBGCMuso/R/thread-4 diff --git a/RBBGCMuso/R/getMeteoData1.R b/RBBGCMuso/R/getMeteoData1.R deleted file mode 100644 index 37fa55e..0000000 --- a/RBBGCMuso/R/getMeteoData1.R +++ /dev/null @@ -1,533 +0,0 @@ -#' getMeteoData1BGC -#' -#' This function downloads hourly 2m air temperature and total precipitation values in NetCDF file format -#' at one grid point to create MTClim files which contain daily data. -#' Please note that, to download ERA5, you need Copernicus registration and API key. -#' According to default settings, data will be downloaded in the gridpoint of Martonvasar (Hungary), in 2017. -#' -#' @author Erzsebet Kristof, Roland Hollos -#' @param startYear Start year of the downloading data. It shall be greater than 1978. -#' @param endYear End year of the downloading data. It shall be smaller than 2020. -#' @param lon Geographical longitude of the selected grid point (negative values: W, positive values: E). It shall be between -180 and 180. The value is rounded to two decimal places. -#' @param lat Geographical latitude of the selected grid point (negative values: S, positive values: N). It shall be between -90 and 90. The value is rounded to two decimal places. -#' @param timeOut Time in seconds to wait to download from Copernicus. -#' @param monthList Month selector (e.g. to download data for April then write "01", to download data for September and October then write c("04","05") -#' @param dayList Day selector (e.g. to download data for the 1st day of the month(s) then write "01", to download data for the 2nd and 4th day of the month then write c("02","04") -#' @param hourList Hour selector (e.g. to download data for 00 UTC then write "00:00", to download data for 01 UTC and 15 UTC then write c("01","15") -#' @param fileDir Directory where the .ini and .mtcin files will be created -#' @param destDir Directory where the wth file will be created. If it is NULL then destDir is the same as fileDir. -#' @param apiFile Directory where the cdsapirc file is located -#' @importFrom ecmwfr wf_set_key wf_request -#' @importFrom ncdf4 nc_open ncvar_get nc_close ncvar_put ncdim_def ncvar_def nc_create -#' @importFrom tcltk tk_choose.dir tk_choose.files -#' @export - -getMeteoData1BGC <- function(startYear=2017, endYear=2017, lon=18.8, lat=47.3, timeOut=7200, - monthList=sprintf("%02d",1:12), - dayList=sprintf("%02d",1:31), - hourList=sprintf("%02d:00",0:23), - destDir=NULL, apiFile=NULL, fileDir=NULL) { - - if ((startYear<1979)|(endYear>2019)) { - stop("Error, please choose a year between 1979 and 2019.") - } - - if (startYear>endYear) { - stop("Error, start year shall be larger than end year.") - } - - if(is.null(destDir)){ - destDir <- "mtcFiles" - } - - # print("Please choose the file in the pop-up window which contains the CDS API key.") - - # With tk_choose.files, it is not working in R Server. - if(is.null(apiFile)){ - apiFile <- tk_choose.files(caption = "Please choose the file which contains the CDS API key.") - print(sprintf("Apifile: %s is selected",apiFile)) - } - - if(is.null(fileDir)){ - fileDir <- tk_choose.dir(caption = "Please choose the fileDir.") - } - - apiCodes <- suppressWarnings(readLines(apiFile)[2]) - apiCodes <- unlist(strsplit(apiCodes,split=":")) - userID <- trimws(apiCodes[2]) - APIkey <- trimws(apiCodes[3]) - - # print("Please choose an empty folder in the pop-up window in which MTClim file will be generated.") - currDir <- getwd() - # fileDir <- tk_choose.dir(caption = "Please choose an empty folder in which MTClim file will be generated.") - setwd(fileDir) - dir.create(destDir) - - lonList <- expand.grid(lon,lat)[,1] - latList <- expand.grid(lon,lat)[,2] - - for (grid in 1:length(lonList)) { - - if ((lonList[grid]>=180.01)|(lonList[grid]<=-180.01)) { - stop("Error, please choose a geographical longitude between -180 and 180.") - } - - if ((latList[grid]>=90.01)|(latList[grid]<=-90.01)) { - stop("Error, please choose a geographical latitude between -90 and 90.") - } - - lon <- round(lonList[grid], digits=2) - lat <- round(latList[grid], digits=2) - - lonSub <- gsub("[.]", "_", as.character(lonList[grid])) - latSub <- gsub("[.]", "_", as.character(latList[grid])) - - print("-----------------------------------------------------------------") - print("1/6. WAITING TO COPERNICUS TO VALIDATE USER ID...") - # Download 2m temperature data from ERA5: - userID <- wf_set_key(user = userID, - key = APIkey, - service = "cds") - - print("-----------------------------------------------------------------") - print("2/6. DOWNLOADING TEMPERATURE DATASET FROM ERA5 DATABASE. PLEASE WAIT...") - - request <- list( - "dataset" = "reanalysis-era5-single-levels", - "class" = "ea", - "expver" = "1", - "stream" = "oper", - "product_type" = "reanalysis", - "variable" = "167.128", - "year" = as.character(startYear:endYear), - "month" = monthList, - "day" = dayList, - "time" = hourList, - "area" = paste0(lat,"/",lon,"/",lat,"/",lon), # N, W, S, E - "grid" = "0.1/0.1", - "format" = "netcdf", - "target" = paste0("ERA5_hourly_t2m_",startYear,"-", - endYear,"_lat_",lat,"_lon_",lon,".nc")) - - file <- wf_request(user = userID, # user ID (for authentification) - request = request, # the request - transfer = TRUE, # download the file - path = fileDir, - time_out = timeOut) - - print("-----------------------------------------------------------------") - print("3/6. DOWNLOADING PRECIPITATION DATASET FROM ERA5 DATABASE. PLEASE WAIT...") - - # Download total precipitation data from ERA5: - request <- list( - "dataset" = "reanalysis-era5-single-levels", - "class" = "ea", - "expver" = "1", - "stream" = "oper", - "product_type" = "reanalysis", - "variable" = "228.128", - "year" = as.character(startYear:endYear), - "month" = monthList, - "day" = dayList, - "time" = hourList, - "area" = paste0(lat,"/",lon,"/",lat,"/",lon), # N, W, S, E - "grid" = "0.1/0.1", - "format" = "netcdf", - "target" = paste0("ERA5_hourly_precip_",startYear,"-",endYear,"_lat_",lat,"_lon_",lon,".nc")) - - file <- wf_request(user = userID, # user ID (for authentification) - request = request, # the request - transfer = TRUE, # download the file - path = fileDir, - time_out = timeOut) - - print("-----------------------------------------------------------------") - print("4/6. DOWNLOADING GEOPOTENTIAL DATASET FROM ERA5 DATABASE. PLEASE WAIT...") - - # Download geopotential height data from ERA5: - request <- list( - "dataset" = "reanalysis-era5-single-levels", - "class" = "ea", - "expver" = "1", - "stream" = "oper", - "product_type" = "reanalysis", - "variable" = "129.128", - "year" = as.character(startYear:endYear), - "month" = monthList, - "day" = dayList, - "time" = hourList, - "area" = paste0(lat,"/",lon,"/",lat,"/",lon), # N, W, S, E - "grid" = "0.1/0.1", - "format" = "netcdf", - "target" = paste0("ERA5_hourly_geopot_",startYear,"-",endYear,"_lat_",lat,"_lon_",lon,".nc")) - - file <- wf_request(user = userID, # user ID (for authentification) - request = request, # the request - transfer = TRUE, # download the file - path = fileDir, - time_out = timeOut) - - print("-----------------------------------------------------------------") - print("5/6. CREATING MTCLIM FILES...") - - ### Creating MTClim files: - - ## Daily minimum and maximum temperature: - - file.nc <- nc_open(paste0("ERA5_hourly_t2m_",startYear,"-",endYear,"_lat_",lat,"_lon_",lon,".nc")) - longitude <- ncvar_get(file.nc, "longitude") - latitude <- ncvar_get(file.nc, "latitude") - timesteps <- ncvar_get(file.nc, "time") - date <- as.POSIXct(timesteps*3600, origin="1900-01-01 00:00:0.0",tz="UTC") - var <- ncvar_get(file.nc, "t2m")-273.15 # in Celsius - nc_close(file.nc) - - # Handling leap years: - leapYear <- c("1952-12-31","1956-12-31","1960-12-31","1964-12-31","1968-12-31", - "1972-12-31","1976-12-31","1980-12-31","1984-12-31","1988-12-31", - "1992-12-31","1996-12-31","2000-12-31","2004-12-31","2008-12-31", - "2012-12-31","2016-12-31","2020-12-31","2024-12-31") - - leapYearIndex <- list() - for (index in 1:length(leapYear)){ - if((length(which(grepl(leapYear[index], date)))>0)==TRUE){ - leapYearIndex[[index]] <- which(grepl(leapYear[index], date)) - } - } - - leapYearIndex <- unlist(leapYearIndex) - - leapYearIndex <- na.omit(leapYearIndex) - - if(length(leapYearIndex)>0){ - timesteps <- timesteps[-leapYearIndex] - date <- date[-leapYearIndex] - var <- var[-leapYearIndex] - } - - dailyT2mMin <- c(0) - dailyT2mMax <- c(0) - dailyT2mAvg <- c(0) - - if(dim(var)%%24==0) { - for (i in seq(1,dim(date),24)) { - dailyT2mMin[i] <- min(var[i:(i+23)]) - - dailyT2mMax[i] <- max(var[i:(i+23)]) - - dailyT2mAvg[i] <- mean(var[i:(i+23)]) - - } - } - - dailyT2mMin <- na.omit(dailyT2mMin) - dailyT2mMax <- na.omit(dailyT2mMax) - - ## Estimating elevation: - file.nc <- nc_open(paste0("ERA5_hourly_geopot_",startYear,"-",endYear,"_lat_",lat,"_lon_",lon,".nc")) - var2 <- ncvar_get(file.nc, "z") # in m**2/s**2 - nc_close(file.nc) - - DB <- data.frame(0) - DB[1,1] <- 1 - DB[1,2] <- 0 - DB[1,3] <- lat - DB[1,4] <- lon - DB[1,5] <- mean(var2)/9.80665 - DB[1,6] <- 0 - DB[1,7] <- 0 - DB[1,8] <- mean(var) - DB[1,9] <- mean(dailyT2mAvg, na.rm=TRUE) - - colnames(DB) <- c("new_id", "100m_rep_id", "lat (deg)", "long (deg)", - "elev (m)", "fboritas", "nuts3", "Tavg", "dTavg") - - # Save lon, lat values into csv: - write.table(lon, "Lon.csv", sep=";", row.names=FALSE, col.names=FALSE) - write.table(lat, "Lat.csv", sep=";", row.names=FALSE, col.names=FALSE) - write.table(DB, "DB_GRID_TABLE.csv", sep=";", row.names=FALSE) - - # Save the data into NetCDF files: tmin and tmax - - lon_nc <- longitude - lat_nc <- latitude - time_nc <- timesteps[seq(1,dim(timesteps),24)] - - longitude <- ncdim_def(name="longitude", units="degree", vals=lon_nc, longname="geographical longitude") - latitude <- ncdim_def(name="latitude", units="degree", vals=lat_nc, longname="geographical latitude") - time <- ncdim_def(name="time", units="day", vals=time_nc, longname="timesteps") - - mv <- -999 - - variable_Tmin <- ncvar_def(name="Tmin", units="Celsius", dim=list(longitude,latitude,time), - missval=mv, longname="Daily minimum temperature", prec="double") - - variable_Tmax <- ncvar_def(name="Tmax", units="Celsius", dim=list(longitude,latitude,time), - missval=mv, longname="Daily maximum temperature", prec="double") - - nc_min <- nc_create(filename=paste0("ERA5_daily_t2m_min_",startYear,"-",endYear,"_lat_",lat,"_lon_",lon,".nc"), vars=list(variable_Tmin)) - nc_max <- nc_create(filename=paste0("ERA5_daily_t2m_max_",startYear,"-",endYear,"_lat_",lat,"_lon_",lon,".nc"), vars=list(variable_Tmax)) - - ncvar_put(nc_min, varid=variable_Tmin, vals=dailyT2mMin, start=c(1,1,1), - count=c(dim(lon_nc),dim(lat_nc),length(timesteps[seq(1,dim(timesteps),24)]))) - - ncvar_put(nc_max, varid=variable_Tmax, vals=dailyT2mMax, start=c(1,1,1), - count=c(dim(lon_nc),dim(lat_nc),length(timesteps[seq(1,dim(timesteps),24)]))) - - nc_close(nc_min) - nc_close(nc_max) - - ## Daily precipitation sums: - - file.nc <- nc_open(paste0("ERA5_hourly_precip_",startYear,"-",endYear,"_lat_",lat,"_lon_",lon,".nc")) - longitude <- ncvar_get(file.nc, "longitude") - latitude <- ncvar_get(file.nc, "latitude") - timesteps <- ncvar_get(file.nc, "time") - date <- as.POSIXct(timesteps*3600, origin="1900-01-01 00:00:0.0",tz="UTC") - var <- ncvar_get(file.nc, "tp")*1000 # in mm - nc_close(file.nc) - - # Handling leap years: - leapYearIndex <- list() - for (index in 1:length(leapYear)){ - if((length(which(grepl(leapYear[index], date)))>0)==TRUE){ - leapYearIndex[[index]] <- which(grepl(leapYear[index], date)) - } - } - - leapYearIndex <- unlist(leapYearIndex) - - leapYearIndex <- na.omit(leapYearIndex) - - if(length(leapYearIndex)>0){ - timesteps <- timesteps[-leapYearIndex] - date <- date[-leapYearIndex] - var <- var[-leapYearIndex] - } - - dailyPrecipSum <- c(0) - - if(dim(var)%%24==0) { - for (i in seq(1,dim(date),24)) { - dailyPrecipSum[i] <- sum(var[i:(i+23)]) - } - } - - dailyPrecipSum <- na.omit(dailyPrecipSum) - summary(dailyPrecipSum) - length(dailyPrecipSum) - length(dailyPrecipSum) - - # Save the data into NetCDF files: precip - - # Dimensions: - lon_nc <- longitude - lat_nc <- latitude - time_nc <- timesteps[seq(1,dim(timesteps),24)] - - longitude <- ncdim_def(name="longitude", units="degree", vals=lon_nc, longname="geographical longitude") - latitude <- ncdim_def(name="latitude", units="degree", vals=lat_nc, longname="geographical latitude") - time <- ncdim_def(name="time", units="day", vals=time_nc, longname="timesteps") - - mv <- -999 - - variable <- ncvar_def(name="prcp", units="mm", dim=list(longitude,latitude,time), - missval=mv, longname="Daily precipitation sums", prec="double") - - nc <- nc_create(filename=paste0("ERA5_daily_precip_sum_",startYear,"-",endYear,"_lat_",lat,"_lon_",lon,".nc"), vars=list(variable)) - - ncvar_put(nc, varid=variable, vals=dailyPrecipSum, start=c(1,1,1), - count=c(dim(lon_nc),dim(lat_nc),length(timesteps[seq(1,dim(timesteps),24)]))) - - nc_close(nc) - - listifyNCDF <- function(ncFile, dataLab, lat, lon, gridCoords){ - dataCon <- nc_open(ncFile) - dataVal <- ncvar_get(dataCon, dataLab) - lonVal <- lon - latVal <- lat - dataCoords <- expand.grid(latVal[1],lonVal[1]) - numberOfGridCells <- nrow(gridCoords) - - dataOut <- list() - for(i in 1: numberOfGridCells){ - dataOut[[i]] <- dataVal - } - return(dataOut) - } - - daily_precip.nc <- nc_open(paste0("ERA5_daily_precip_sum_", startYear, - "-", endYear, "_lat_", lat, "_lon_", lon, ".nc")) - daily_tmin.nc <- nc_open(paste0("ERA5_daily_t2m_min_", startYear, - "-", endYear, "_lat_", lat, "_lon_", lon, ".nc")) - daily_tmax.nc <- nc_open(paste0("ERA5_daily_t2m_max_", startYear, - "-", endYear, "_lat_", lat, "_lon_", lon, ".nc")) - - tmin <- ncvar_get(daily_tmin.nc, "Tmin") - precip <- ncvar_get(daily_precip.nc, "prcp") - tmax <- ncvar_get(daily_tmax.nc, "Tmax") - - nc_close(daily_precip.nc) - nc_close(daily_tmin.nc) - nc_close(daily_tmax.nc) - - lon <- read.csv2("Lon.csv", stringsAsFactors = FALSE, header = FALSE) - lat <- read.csv2("Lat.csv", stringsAsFactors = FALSE, header = FALSE) - gridCoords <- read.csv2("DB_GRID_TABLE.csv", stringsAsFactors = FALSE) - gridSkeleton <- gridCoords[,c(1,3,4)] - - precip <- listifyNCDF(ncFile = paste0("ERA5_daily_precip_sum_",startYear,"-", - endYear,"_lat_",lat,"_lon_",lon,".nc"), - dataLab = "prcp", # tp - lon = lon, - lat = lat, - gridCoords = gridSkeleton) - precip <- lapply(precip, function(x) x/10) - - tmin <- listifyNCDF(ncFile = paste0("ERA5_daily_t2m_min_",startYear,"-", - endYear,"_lat_",lat,"_lon_",lon,".nc"), - dataLab = "Tmin", - lon = lon, - lat = lat, - gridCoords = gridSkeleton) - - tmax <- listifyNCDF(ncFile = paste0("ERA5_daily_t2m_max_",startYear,"-", - endYear,"_lat_",lat,"_lon_",lon,".nc"), - dataLab = "Tmax", - lon = lon, - lat = lat, - gridCoords = gridSkeleton) - numberOfSites <- length(tmax) - tmp <- tmin - - mtcMaker <- function(simulationName = "", executable, minYear, startYear, numYears, - maxTemperatures, minTemperatures, precipitations, sites, destDir){ - # executable <- normalizePath(executable) # make absolute from relative path - file.copy(executable,getwd()) - executable <- paste0("./",basename(executable)) - iniTemp <- c(paste0("Simulation is based on: ",simulationName),"","IOFILES", - "1.mtcin", - "1", "CONTROL", "0", ((endYear-startYear)+1)*365, #"365", - "0", "0", "1", "PARAMETERS", - "1000.0", "60.0", "45.0", "1000.0", - "0.0", "0.0", "60.0", "0.0", "0.0", - "-6.0", "-3.0", "END") - endYear <- startYear + numYears - 1 - - mtDate <- function(startYear, numYears, endYear){ - ## generate the first two column for mtcin - year <- rep(startYear:endYear, each = 365) - yearday <- rep(1:365,numYears) - return(cbind(year,yearday)) - } - - sliceData <- function(minYear, startYear, numYears, data){ - ## get the data from the needed perion. - baseDiff <- startYear-minYear - startDate <- baseDiff*365+1 - endDate <- (baseDiff+numYears)*365 - return(data[startDate:endDate]) - } - - nSites <- nrow(sites) - for(i in 1:nSites) { - print(i) - iniTemp[4] <- sprintf("%d.mtcin",i) - iniTemp[8] <- numYears*365 - iniTemp[5] <- file.path(destDir,i) - iniTemp[13] <- iniTemp[16] <- sites[i,5] - iniTemp[15] <- sites[i,3] #LAT is the first col. - - write.table(iniTemp,paste(i, ".ini", sep=""),col.names=FALSE,row.names=FALSE, - sep="\n",quote=FALSE) - - itmax <- sliceData(minYear,startYear,numYears,maxTemperatures[[i]]) - itmin <- sliceData(minYear,startYear,numYears,minTemperatures[[i]]) - iprecip <- sliceData(minYear,startYear,numYears,precipitations[[i]]) - dateData <- mtDate(startYear,numYears,endYear) - - write.table(cbind(dateData,itmax,itmin,iprecip),paste(i, ".mtcin", sep="") - ,row.names=FALSE,col.names=FALSE) - - system2(executable, paste0(i, ".ini"), stdout=NULL) - } - } - - executable <- if(Sys.info()["sysname"]=="Windows"){ - paste0(system.file("",package="MTClimMaker"),"/mtclim43.exe") - } else { - paste0(system.file("",package="MTClimMaker"),"/mtclim43") - } - - mtcMaker(simulationName = paste0("lat. ", gsub("_",".",latSub), ", lon. ", gsub("_",".",lonSub), ", elev. ", - round(as.numeric(gridCoords[1,5]), digits=1), " m, ", "Note that VPD and srad are computed by MTClim."), - executable = executable, - # executable = paste0(system.file("",package="MTClimMaker"),"/mtclim43"), - minYear = startYear, - startYear = startYear, - numYears = (endYear-startYear)+1, - maxTemperatures = tmax, - minTemperatures = tmin, - precipitations = precip, - sites = gridCoords, - destDir = destDir) - setwd(currDir) - - print("-----------------------------------------------------------------") - # print("6/6. MOVING .mtcin, .ini AND .wth FILES INTO THE CORRECT FOLDER.") - print("6/6. RENAMING .mtcin, .ini AND .wth FILES...") - - - ### Move the ini, mtcin, mtc43 files into the correct folder. - # Moving files to the correct folder: - # dir.create(paste0("lat", latSub, "_lon", lonSub)) - # file.copy("1.ini", paste0("lat", latSub, "_lon", lonSub)) - # file.copy("1.mtcin", paste0("lat", latSub, "_lon", lonSub)) - # file.copy(paste0(destDir,"/1.mtc43"), paste0("lat", latSub, "_lon", lonSub)) - - # file.rename(paste0("lat", latSub, "_lon", lonSub, "/1.ini"), - # paste0("lat", latSub, "_lon", lonSub, "/", grid, - # "_lat", latSub, "_lon", lonSub, ".ini")) - # - # file.rename(paste0("lat", latSub, "_lon", lonSub, "/1.mtcin"), - # paste0("lat", latSub, "_lon", lonSub, "/", grid, - # "_lat", latSub, "_lon", lonSub, ".mtcin")) - # - # file.rename(paste0("lat", latSub, "_lon", lonSub, "/1.wth"), - # paste0("lat", latSub, "_lon", lonSub, "/", grid, - # "_lat", latSub, "_lon", lonSub, ".wth")) - ### - - #file.rename("1.ini", - # paste0(grid, - # "_lat", latSub, "_lon", lonSub, ".ini")) # ID number (1,2,...) can be added to the title by using object grid befor "_lat" - - #file.rename("1.mtcin", - # paste0(grid, - # "_lat", latSub, "_lon", lonSub, ".mtcin")) # ID number (1,2,...) can be added to the title by using object grid befor "_lat" - - file.rename(paste0(fileDir, "/", destDir, "/1.mtc43"), - paste0(fileDir, "/", destDir, "/", - "_lat", latSub, "_lon", lonSub, ".wth")) # ID number (1,2,...) can be added to the title by using object grid befor "_lat" - - file.copy(paste0(fileDir, "/", destDir, "/", "_lat", latSub, "_lon", lonSub, ".wth"), fileDir) - # file.rename("1.ini", paste0(sprintf('%0.4d', gridID),".ini")) - # - # file.rename("1.mtcin", paste0(sprintf('%0.4d', gridID),".mtcin")) - # - # file.rename("1.wth", paste0(sprintf('%0.4d', gridID),".wth")) - - print(paste0("MTCLIM FILES ARE CREATED FOR THE GRID POINT: LAT ", - latList[grid], " & LON ", lonList[grid], ".")) - - print(paste0("SUMMARY: ", - "elevation: ", round(as.numeric(gridCoords[1,5]), digits=1), " m, ", - "yearly average temperature: ", round(as.numeric(gridCoords[1,8]), digits=1), " C")) - - # Removing unnecessary files: - file.remove("DB_GRID_TABLE.csv", "Lat.csv", "Lon.csv", "1.ini", "1.mtcin", - "mtclim43.exe", list.files(pattern="ERA5_hourly*"), list.files(pattern="ERA5_daily*")) - } - - unlink(destDir, recursive=TRUE) -} diff --git a/RBBGCMuso/R/thread-1 b/RBBGCMuso/R/thread-1 deleted file mode 100644 index 29d6383..0000000 --- a/RBBGCMuso/R/thread-1 +++ /dev/null @@ -1 +0,0 @@ -100 diff --git a/RBBGCMuso/R/thread-2 b/RBBGCMuso/R/thread-2 deleted file mode 100644 index 29d6383..0000000 --- a/RBBGCMuso/R/thread-2 +++ /dev/null @@ -1 +0,0 @@ -100 diff --git a/RBBGCMuso/R/thread-3 b/RBBGCMuso/R/thread-3 deleted file mode 100644 index 29d6383..0000000 --- a/RBBGCMuso/R/thread-3 +++ /dev/null @@ -1 +0,0 @@ -100 diff --git a/RBBGCMuso/R/thread-4 b/RBBGCMuso/R/thread-4 deleted file mode 100644 index 29d6383..0000000 --- a/RBBGCMuso/R/thread-4 +++ /dev/null @@ -1 +0,0 @@ -100 From 89e5f6933a63371479c90c8f833e16046b633a5a Mon Sep 17 00:00:00 2001 From: Hollos Roland Date: Thu, 24 Sep 2020 18:19:51 +0200 Subject: [PATCH 5/7] fixing meteorology issue --- RBBGCMuso/NAMESPACE | 11 ------- RBBGCMuso/man/calibrateMuso.Rd | 6 +++- RBBGCMuso/man/changemulline.Rd | 2 +- RBBGCMuso/man/getMeteoData1BGC.Rd | 52 ------------------------------- 4 files changed, 6 insertions(+), 65 deletions(-) delete mode 100644 RBBGCMuso/man/getMeteoData1BGC.Rd diff --git a/RBBGCMuso/NAMESPACE b/RBBGCMuso/NAMESPACE index 76126a4..9939ede 100644 --- a/RBBGCMuso/NAMESPACE +++ b/RBBGCMuso/NAMESPACE @@ -13,7 +13,6 @@ export(createSoilFile) export(getAnnualOutputList) export(getConstMatrix) export(getDailyOutputList) -export(getMeteoData1BGC) export(getyearlycum) export(getyearlymax) export(mtclim) @@ -54,8 +53,6 @@ importFrom(dplyr,mutate) importFrom(dplyr,select) importFrom(dplyr,summarize) importFrom(dplyr,tbl_df) -importFrom(ecmwfr,wf_request) -importFrom(ecmwfr,wf_set_key) importFrom(future,future) importFrom(ggplot2,aes) importFrom(ggplot2,aes_string) @@ -85,19 +82,11 @@ importFrom(limSolve,xsample) importFrom(lubridate,leap_year) importFrom(magrittr,'%<>%') importFrom(magrittr,'%>%') -importFrom(ncdf4,nc_close) -importFrom(ncdf4,nc_create) -importFrom(ncdf4,nc_open) -importFrom(ncdf4,ncdim_def) -importFrom(ncdf4,ncvar_def) -importFrom(ncdf4,ncvar_get) -importFrom(ncdf4,ncvar_put) importFrom(openxlsx,read.xlsx) importFrom(rmarkdown,pandoc_version) importFrom(rmarkdown,render) importFrom(scales,percent) importFrom(stats,approx) -importFrom(tcltk,tk_choose.dir) importFrom(tcltk,tk_choose.files) importFrom(tibble,rownames_to_column) importFrom(tidyr,gather) diff --git a/RBBGCMuso/man/calibrateMuso.Rd b/RBBGCMuso/man/calibrateMuso.Rd index 5710638..8530339 100644 --- a/RBBGCMuso/man/calibrateMuso.Rd +++ b/RBBGCMuso/man/calibrateMuso.Rd @@ -27,7 +27,11 @@ calibrateMuso( numCores = (parallel::detectCores() - 1), pb = txtProgressBar(min = 0, max = iterations, style = 3), maxLikelihoodEpc = TRUE, - pbUpdate = setTxtProgressBar + pbUpdate = setTxtProgressBar, + method = "GLUE", + lg = FALSE, + w = NULL, + ... ) } \description{ diff --git a/RBBGCMuso/man/changemulline.Rd b/RBBGCMuso/man/changemulline.Rd index 848e443..3676348 100644 --- a/RBBGCMuso/man/changemulline.Rd +++ b/RBBGCMuso/man/changemulline.Rd @@ -4,7 +4,7 @@ \alias{changemulline} \title{changemulline} \usage{ -changemulline(filePaths, calibrationPar, contents, src) +changemulline(filePaths, calibrationPar, contents, src, outFiles = filePaths) } \description{ The function uses the previous changspecline function to operate. diff --git a/RBBGCMuso/man/getMeteoData1BGC.Rd b/RBBGCMuso/man/getMeteoData1BGC.Rd deleted file mode 100644 index 6f4265e..0000000 --- a/RBBGCMuso/man/getMeteoData1BGC.Rd +++ /dev/null @@ -1,52 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/getMeteoData1.R -\name{getMeteoData1BGC} -\alias{getMeteoData1BGC} -\title{getMeteoData1BGC} -\usage{ -getMeteoData1BGC( - startYear = 2017, - endYear = 2017, - lon = 18.8, - lat = 47.3, - timeOut = 7200, - monthList = sprintf("\%02d", 1:12), - dayList = sprintf("\%02d", 1:31), - hourList = sprintf("\%02d:00", 0:23), - destDir = NULL, - apiFile = NULL, - fileDir = NULL -) -} -\arguments{ -\item{startYear}{Start year of the downloading data. It shall be greater than 1978.} - -\item{endYear}{End year of the downloading data. It shall be smaller than 2020.} - -\item{lon}{Geographical longitude of the selected grid point (negative values: W, positive values: E). It shall be between -180 and 180. The value is rounded to two decimal places.} - -\item{lat}{Geographical latitude of the selected grid point (negative values: S, positive values: N). It shall be between -90 and 90. The value is rounded to two decimal places.} - -\item{timeOut}{Time in seconds to wait to download from Copernicus.} - -\item{monthList}{Month selector (e.g. to download data for April then write "01", to download data for September and October then write c("04","05")} - -\item{dayList}{Day selector (e.g. to download data for the 1st day of the month(s) then write "01", to download data for the 2nd and 4th day of the month then write c("02","04")} - -\item{hourList}{Hour selector (e.g. to download data for 00 UTC then write "00:00", to download data for 01 UTC and 15 UTC then write c("01","15")} - -\item{destDir}{Directory where the wth file will be created. If it is NULL then destDir is the same as fileDir.} - -\item{apiFile}{Directory where the cdsapirc file is located} - -\item{fileDir}{Directory where the .ini and .mtcin files will be created} -} -\description{ -This function downloads hourly 2m air temperature and total precipitation values in NetCDF file format -at one grid point to create MTClim files which contain daily data. -Please note that, to download ERA5, you need Copernicus registration and API key. -According to default settings, data will be downloaded in the gridpoint of Martonvasar (Hungary), in 2017. -} -\author{ -Erzsebet Kristof, Roland Hollos -} From 33abd4dff3ffc214ab0fd420bcae45133eb848dd Mon Sep 17 00:00:00 2001 From: Hollos Roland Date: Sat, 26 Sep 2020 07:59:13 +0200 Subject: [PATCH 6/7] Removing C/C++ dependencies --- RBBGCMuso/DESCRIPTION | 2 - RBBGCMuso/NAMESPACE | 5 - RBBGCMuso/R/RcppExports.R | 52 - RBBGCMuso/R/calibrateMuso.R | 7 +- RBBGCMuso/man/changeMusoC.Rd | 16 - RBBGCMuso/man/mtclim.Rd | 15 - RBBGCMuso/man/musoRandomizer.Rd | 17 - RBBGCMuso/src/RcppExports.cpp | 89 -- RBBGCMuso/src/ini.h | 156 -- RBBGCMuso/src/mtc.cpp | 28 - RBBGCMuso/src/mtclim43.c | 2223 --------------------------- RBBGCMuso/src/mtclim43_constants.h | 29 - RBBGCMuso/src/mtclim43_parameters.h | 35 - RBBGCMuso/src/musoIO.cpp | 138 -- RBBGCMuso/src/musoRandomizer.cpp | 245 --- 15 files changed, 5 insertions(+), 3052 deletions(-) delete mode 100644 RBBGCMuso/R/RcppExports.R delete mode 100644 RBBGCMuso/man/changeMusoC.Rd delete mode 100644 RBBGCMuso/man/mtclim.Rd delete mode 100644 RBBGCMuso/man/musoRandomizer.Rd delete mode 100644 RBBGCMuso/src/RcppExports.cpp delete mode 100644 RBBGCMuso/src/ini.h delete mode 100644 RBBGCMuso/src/mtc.cpp delete mode 100644 RBBGCMuso/src/mtclim43.c delete mode 100644 RBBGCMuso/src/mtclim43_constants.h delete mode 100644 RBBGCMuso/src/mtclim43_parameters.h delete mode 100644 RBBGCMuso/src/musoIO.cpp delete mode 100644 RBBGCMuso/src/musoRandomizer.cpp diff --git a/RBBGCMuso/DESCRIPTION b/RBBGCMuso/DESCRIPTION index c267c32..0624fa2 100644 --- a/RBBGCMuso/DESCRIPTION +++ b/RBBGCMuso/DESCRIPTION @@ -31,8 +31,6 @@ Imports: openxlsx, ncdf4, tcltk -LinkingTo: Rcpp -SystemRequirements: C++11 Maintainer: Roland Hollo's RoxygenNote: 7.1.0 Suggests: knitr, diff --git a/RBBGCMuso/NAMESPACE b/RBBGCMuso/NAMESPACE index 9939ede..19adfe4 100644 --- a/RBBGCMuso/NAMESPACE +++ b/RBBGCMuso/NAMESPACE @@ -2,7 +2,6 @@ export(calibMuso) export(calibrateMuso) -export(changeMusoC) export(changemulline) export(checkMeteoBGC) export(cleanupMuso) @@ -15,7 +14,6 @@ export(getConstMatrix) export(getDailyOutputList) export(getyearlycum) export(getyearlymax) -export(mtclim) export(musoDate) export(musoGlue) export(musoMapping) @@ -23,7 +21,6 @@ export(musoMappingFind) export(musoMonte) export(musoQuickEffect) export(musoRand) -export(musoRandomizer) export(musoSensi) export(normalMuso) export(optiMuso) @@ -41,7 +38,6 @@ export(supportedMuso) export(updateMusoMapping) import(ggplot2) import(utils) -importFrom(Rcpp,evalCpp) importFrom(data.table,':=') importFrom(data.table,data.table) importFrom(data.table,fread) @@ -91,4 +87,3 @@ importFrom(tcltk,tk_choose.files) importFrom(tibble,rownames_to_column) importFrom(tidyr,gather) importFrom(tidyr,separate) -useDynLib(RBBGCMuso) diff --git a/RBBGCMuso/R/RcppExports.R b/RBBGCMuso/R/RcppExports.R deleted file mode 100644 index 3abb296..0000000 --- a/RBBGCMuso/R/RcppExports.R +++ /dev/null @@ -1,52 +0,0 @@ -# Generated by using Rcpp::compileAttributes() -> do not edit by hand -# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 - -#' mtclim -#' -#' This is the core mtclim function -#' @importFrom Rcpp evalCpp -#' @useDynLib RBBGCMuso -#' @param iniFile is the name of the inifile -#' @keywords internal -#' @export -mtclim <- function(iniFile) { - invisible(.Call('_RBBGCMuso_mtclim', PACKAGE = 'RBBGCMuso', iniFile)) -} - -getWritePositions <- function(a) { - .Call('_RBBGCMuso_getWritePositions', PACKAGE = 'RBBGCMuso', a) -} - -#' changeMusoC -#' -#' This function is fastly randomize values based on min and max values -#' @importFrom Rcpp evalCpp -#' @useDynLib RBBGCMuso -#' @param inFile is the big matrix -#' @param outFile is the small matrix -#' @export -changeMusoC <- function(inFile, outFile, inMat) { - invisible(.Call('_RBBGCMuso_changeMusoC', PACKAGE = 'RBBGCMuso', inFile, outFile, inMat)) -} - -randTypeOne <- function(m) { - .Call('_RBBGCMuso_randTypeOne', PACKAGE = 'RBBGCMuso', m) -} - -randTypeTwo <- function(m) { - .Call('_RBBGCMuso_randTypeTwo', PACKAGE = 'RBBGCMuso', m) -} - -#' musoRandomizer -#' -#' This function is fastly randomize values based on min and max values, -#' and row indexes. -#' @importFrom Rcpp evalCpp -#' @useDynLib RBBGCMuso -#' @param A is the big matrix -#' @param B is the small matrix -#' @export -musoRandomizer <- function(A, B) { - .Call('_RBBGCMuso_musoRandomizer', PACKAGE = 'RBBGCMuso', A, B) -} - diff --git a/RBBGCMuso/R/calibrateMuso.R b/RBBGCMuso/R/calibrateMuso.R index 2ac1992..89c643e 100644 --- a/RBBGCMuso/R/calibrateMuso.R +++ b/RBBGCMuso/R/calibrateMuso.R @@ -125,9 +125,12 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, }, "agromo"={ liks <- results[,sprintf("%s_likelihood",names(likelihood))] - epcIndexes <- value(fut[[1]], stdout = FALSE, signal=FALSE) + epcIndexes <- future::value(fut[[1]], stdout = FALSE, signal=FALSE) epcVals <- results[which.max(liks),1:length(epcIndexes)] - changemulline(filePaths= settings$epcInput[2], epcIndexes, epcVals, src = settings$epcInput[2], outFiles = "maxLikelihood_epc.epc") + epcPlace <- file.path(dirname(settings$inputFiles),settings$epc)[2] + changemulline(filePaths= epcPlace, epcIndexes, + epcVals, src = settings$epcInput[2], + outFiles = "maxLikelihood_epc.epc") names(epcVals) <- epcIndexes xdate <- as.Date(measuredData$date) meanM <- measuredData[,sprintf("mean.%s", names(likelihood))] diff --git a/RBBGCMuso/man/changeMusoC.Rd b/RBBGCMuso/man/changeMusoC.Rd deleted file mode 100644 index 270a7f1..0000000 --- a/RBBGCMuso/man/changeMusoC.Rd +++ /dev/null @@ -1,16 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{changeMusoC} -\alias{changeMusoC} -\title{changeMusoC} -\usage{ -changeMusoC(inFile, outFile, inMat) -} -\arguments{ -\item{inFile}{is the big matrix} - -\item{outFile}{is the small matrix} -} -\description{ -This function is fastly randomize values based on min and max values -} diff --git a/RBBGCMuso/man/mtclim.Rd b/RBBGCMuso/man/mtclim.Rd deleted file mode 100644 index 5adbaf2..0000000 --- a/RBBGCMuso/man/mtclim.Rd +++ /dev/null @@ -1,15 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{mtclim} -\alias{mtclim} -\title{mtclim} -\usage{ -mtclim(iniFile) -} -\arguments{ -\item{iniFile}{is the name of the inifile} -} -\description{ -This is the core mtclim function -} -\keyword{internal} diff --git a/RBBGCMuso/man/musoRandomizer.Rd b/RBBGCMuso/man/musoRandomizer.Rd deleted file mode 100644 index e731389..0000000 --- a/RBBGCMuso/man/musoRandomizer.Rd +++ /dev/null @@ -1,17 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{musoRandomizer} -\alias{musoRandomizer} -\title{musoRandomizer} -\usage{ -musoRandomizer(A, B) -} -\arguments{ -\item{A}{is the big matrix} - -\item{B}{is the small matrix} -} -\description{ -This function is fastly randomize values based on min and max values, -and row indexes. -} diff --git a/RBBGCMuso/src/RcppExports.cpp b/RBBGCMuso/src/RcppExports.cpp deleted file mode 100644 index 83ef4fc..0000000 --- a/RBBGCMuso/src/RcppExports.cpp +++ /dev/null @@ -1,89 +0,0 @@ -// Generated by using Rcpp::compileAttributes() -> do not edit by hand -// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 - -#include - -using namespace Rcpp; - -// mtclim -void mtclim(std::string iniFile); -RcppExport SEXP _RBBGCMuso_mtclim(SEXP iniFileSEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< std::string >::type iniFile(iniFileSEXP); - mtclim(iniFile); - return R_NilValue; -END_RCPP -} -// getWritePositions -IntegerVector getWritePositions(double a); -RcppExport SEXP _RBBGCMuso_getWritePositions(SEXP aSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type a(aSEXP); - rcpp_result_gen = Rcpp::wrap(getWritePositions(a)); - return rcpp_result_gen; -END_RCPP -} -// changeMusoC -void changeMusoC(std::string inFile, std::string outFile, NumericMatrix inMat); -RcppExport SEXP _RBBGCMuso_changeMusoC(SEXP inFileSEXP, SEXP outFileSEXP, SEXP inMatSEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< std::string >::type inFile(inFileSEXP); - Rcpp::traits::input_parameter< std::string >::type outFile(outFileSEXP); - Rcpp::traits::input_parameter< NumericMatrix >::type inMat(inMatSEXP); - changeMusoC(inFile, outFile, inMat); - return R_NilValue; -END_RCPP -} -// randTypeOne -NumericMatrix randTypeOne(NumericMatrix m); -RcppExport SEXP _RBBGCMuso_randTypeOne(SEXP mSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericMatrix >::type m(mSEXP); - rcpp_result_gen = Rcpp::wrap(randTypeOne(m)); - return rcpp_result_gen; -END_RCPP -} -// randTypeTwo -NumericMatrix randTypeTwo(NumericMatrix m); -RcppExport SEXP _RBBGCMuso_randTypeTwo(SEXP mSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericMatrix >::type m(mSEXP); - rcpp_result_gen = Rcpp::wrap(randTypeTwo(m)); - return rcpp_result_gen; -END_RCPP -} -// musoRandomizer -NumericMatrix musoRandomizer(NumericMatrix A, NumericMatrix B); -RcppExport SEXP _RBBGCMuso_musoRandomizer(SEXP ASEXP, SEXP BSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP); - Rcpp::traits::input_parameter< NumericMatrix >::type B(BSEXP); - rcpp_result_gen = Rcpp::wrap(musoRandomizer(A, B)); - return rcpp_result_gen; -END_RCPP -} - -static const R_CallMethodDef CallEntries[] = { - {"_RBBGCMuso_mtclim", (DL_FUNC) &_RBBGCMuso_mtclim, 1}, - {"_RBBGCMuso_getWritePositions", (DL_FUNC) &_RBBGCMuso_getWritePositions, 1}, - {"_RBBGCMuso_changeMusoC", (DL_FUNC) &_RBBGCMuso_changeMusoC, 3}, - {"_RBBGCMuso_randTypeOne", (DL_FUNC) &_RBBGCMuso_randTypeOne, 1}, - {"_RBBGCMuso_randTypeTwo", (DL_FUNC) &_RBBGCMuso_randTypeTwo, 1}, - {"_RBBGCMuso_musoRandomizer", (DL_FUNC) &_RBBGCMuso_musoRandomizer, 2}, - {NULL, NULL, 0} -}; - -RcppExport void R_init_RBBGCMuso(DllInfo *dll) { - R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); - R_useDynamicSymbols(dll, FALSE); -} diff --git a/RBBGCMuso/src/ini.h b/RBBGCMuso/src/ini.h deleted file mode 100644 index dc847f7..0000000 --- a/RBBGCMuso/src/ini.h +++ /dev/null @@ -1,156 +0,0 @@ -/* -ini.h -Peter Thornton, NTSG -8/6/93 - -Header file for functions that read *.ini files and either -store values in variables or open files for input or output -General *.ini file format: -value/string (whitespace) comment, not required (\n) -*/ - - - -/* structure definition for filename handling */ -typedef struct { - char name[128]; - FILE *ptr; -} file; - -/* function prototypes */ -int file_open (file *target, char mode); -int scan_value (file ini, void *var, char mode); -int scan_open (file ini,file *target,char mode); - - -/* file_open is the generic file opening routine using the file structure -defined above */ -int file_open (file *target, char mode) -/* Possible values for mode - 'r' for read binary - 'i' for read ascii - 'w' for write binary - 'o' for write ascii -*/ -{ - int ok=1; - - switch (mode) - { - case 'r': - if ((target->ptr = fopen(target->name,"rb")) == NULL) - { - printf("Can't open %s for binary read ... Exiting\n",target->name); - ok=0; - } - break; - - case 'i': - if ((target->ptr = fopen(target->name,"r")) == NULL) - { - printf("Can't open %s for ascii read ... Exiting\n",target->name); - ok=0; - } - break; - - case 'w': - if ((target->ptr = fopen(target->name,"wb")) == NULL) - { - printf("Can't open %s for binary write ... Exiting\n",target->name); - ok=0; - } - break; - - case 'o': - if ((target->ptr = fopen(target->name,"w")) == NULL) - { - printf("Can't open %s for ascii write ... Exiting\n",target->name); - ok=0; - } - break; - - default: - printf("Invalid mode specification for file_open ... Exiting\n"); - ok=0; - } - return(!ok); -} - - -/* scan_value is the generic ascii input function for use with text -initialization files. Reads the first whitespace delimited word on a line, -and discards the remainder of the line. Returns a value depending on the -specified scan type */ -int scan_value (file ini, void *var, char type) -/* Possible values for type - 'i' for integer - 'd' for double - 's' for string -*/ - -{ - int ok_scan; - int ok=1; - - switch (type) - { - case 'i': - ok_scan = fscanf(ini.ptr, "%d%*[^\n]",(int*)var); - if (ok_scan == 0 || ok_scan == EOF) - { - printf("Error reading int value from %s ... exiting\n",ini.name); - ok=0; - } - break; - - case 'd': - ok_scan = fscanf(ini.ptr, "%lf%*[^\n]",(double*)var); - if (ok_scan == 0 || ok_scan == EOF) - { - printf("Error reading double value from %s... exiting\n",ini.name); - ok=0; - } - break; - - case 's': - ok_scan = fscanf(ini.ptr, "%s%*[^\n]",(char*)var); - if (ok_scan == 0 || ok_scan == EOF) - { - printf("Error reading string value from %s... exiting\n",ini.name); - ok=0; - } - break; - - default: - printf("Invalid type specifier for scan_value ... Exiting\n"); - ok=0; - } - return(!ok); -} - -/* combines scan_value with file_open for reading a filename from an -initialization file and then opening it with a specified access mode */ -int scan_open (file ini,file *target,char mode) -/* Possible values for mode - 'r' for read binary - 'i' for read ascii - 'w' for write binary - 'o' for write ascii -*/ -{ - int ok=1; - - if (scan_value(ini,target->name,'s')) - { - printf("Error reading filename from %s... Exiting\n",ini.name); - ok=0; - } - if (ok) - { - if (file_open(target,mode)) - { - ok=0; - } - } - return(!ok); -} diff --git a/RBBGCMuso/src/mtc.cpp b/RBBGCMuso/src/mtc.cpp deleted file mode 100644 index 54e167e..0000000 --- a/RBBGCMuso/src/mtc.cpp +++ /dev/null @@ -1,28 +0,0 @@ -#include -#include -#include -#include -#include -#include - - -using namespace Rcpp; -using namespace std; -extern "C" { - void mtc(char*); - }; - -//' mtclim -//' -//' This is the core mtclim function -//' @importFrom Rcpp evalCpp -//' @useDynLib RBBGCMuso -//' @param iniFile is the name of the inifile -//' @keywords internal -//' @export -// [[Rcpp::export]] -void mtclim(std::string iniFile){ - char *y = new char[iniFile.length() + 1]; // Allocate memory for char array input - std::strcpy(y, iniFile.c_str()); // Copy c++ string to that input. - mtc(y); -} diff --git a/RBBGCMuso/src/mtclim43.c b/RBBGCMuso/src/mtclim43.c deleted file mode 100644 index a619cfe..0000000 --- a/RBBGCMuso/src/mtclim43.c +++ /dev/null @@ -1,2223 +0,0 @@ -/* - -This is a slitely modified version of mtclim prepared to be part of RBBGCMuso. -modification is done by: Roland Hollos (hollorol@gmail.com) - - -mtclim43.c - -MTCLIM -VERSION 4.3 - -Peter Thornton -NTSG, School of Forestry -University of Montana -1/20/00 - -*************************************** -** Questions or comments? Contact... ** -** Dr. Peter E. Thornton ** -** NTSG, School of Forestry ** -** University of Montana ** -** Missoula, MT 59812 ** -** email: peter@ntsg.umt.edu ** -** phone: 406-243-4326 ** -*************************************** - --------------------------------------------------------------------------- -Code History ------------- -Original code written by R.R. Nemani -Updated 4/1/1989 by J.C. Coughlan -Updated 12/23/1989 by Joe Glassy -Updated 1/4/1993 by Raymond Hunt (version 2.1) -Updated 3/26/1997 by Peter Thornton (version 3.0) -Updated 7/28/1997 by Peter Thornton (version 3.1) -Updated 5/7/1998 by Peter Thornton (version 4.0) -Updated 8/1/1998 by Peter Thornton (version 4.1) -Updated 4/20/1999 by Peter Thornton (version 4.2) -Updated 1/20/2000 by Peter Thornton (version 4.3) - -CHANGES FROM VERSION 4.2 TO VERSION 4.3 -These changes based on the results of: -Thornton, P.E., H. Hasenauer, and M.A. White, 2000. Simultaneous -estimation of daily solar radiation and humidity from observed temperature -and precipitation: an application over complex terrain in Austria. -For submission to Ag For Met. -1) ADDED THE SNOW ACCUMULATION/MELT MODEL FOR CORRECTING RADIATION -FOR INFLUENCE OF SNOWPACK. -2) ADDED ALBEDO CORRECTION TO HORIZON OBSTRUCTION CALCULATIONS. -3) HUMIDITY ALGORITHM NOW DEPENDS ON THE ANNUAL RATION PET/PRCP: -FOR VALUES < 2.5, USING THE TMIN=TDEW ALGORITHM, FOR VALUES > 2.5, -USING THE ARID CORRECTION FROM: -Kimball, J.S., S.W. Running, and R. Nemani, 1997. An improved method for -estimating surface humidity from daily minimum temperature. Ag For Met -85:87-98. -4) MOVED ALL MODEL PARAMETERS INTO THE PARAMETERS.H FILE - -CHANGES FROM VERSION 4.1 TO VERSION 4.2 -1) PUT THE SLOPE AND ASPECT CORRECTIONS TO RADIATION BACK IN. THEY -HAD BEEN REMOVED DURING THE DEVELOPMENT OF THE NEW RADIATION CODE. -This includes the estimation for diffuse radiation on sloping surfaces -during periods when the sun is above a flat horizon, but not providing -direct illumination to the slope. Site east and west horizon corrections -have also been reintroduced. - -CHANGES FROM VERSION 4.0 TO VERSION 4.1 -1) ADDITIONAL REVISION OF RADIATION CODE, FOLLOWING SUBMISSION OF AG FOR MET -MANUSCRIPT DESCRIBING ANALYSIS OF SAMSON DATA. The current code -follows exactly the algorithm detailed in -Thornton, P.E. and S.W. Running, 1999. An improved algorithm for estimating -incident daily solar radiation from measurements of temperature, humidity, -and precipitation. Ag For Met 93:211-228. - -changes from version 3.1 to version 4.0: -1) radiation code completely revamped, following analysis of samson -observations for the vemap2 project. -2) includes an iterative algorithm for estimating vapor pressure and -radiation -3) par output no longer an option, since the old par algorithm is suspect -4) 2-day tmin smoothing no longer an option -5) solar constant now calculated as an analytical function of yearday, -instead of the monthly array of values in earlier versions. removes the -discontinuities between months. -6) boxcar function now uses the previous n days of data, making it a -"pulled boxcar" instead of a "centered boxcar" - -Changes from version 2.0 to version 3.1 -Modified to include the following improvements to the original MTCLIM logic: -1) Improved vapor pressure calculation -2) Improved transmissivity calculation -3) Improved daylength calculation - -Some other differences between version 3.1 and previous versions of MTCLIM: -1) No english units option -2) No total or average radiation option -3) No threshold radiation option -5) No pre-rainy days correction for transmissivity -6) No radiation correction to air temperatures -7) No LAI corrections to air temperatures -8) Only one precipitation station allowed -9) Some parameters formerly in initialization file are now in mtclim_const.h - --------------------------------------------------------------------------- -UNITS ------ -Temperatures degrees C -Temp. lapse rates degrees C / 1000 m -Precipitation cm / day -Vapor pressure Pa -Vapor pressure deficit Pa -Radiation W/m2, average over daylight period -Daylength s (sunrise to sunset, flat horizons) -Elevation m -Latitude decimal degrees -Aspect decimal degrees -Slope decimal degrees -E/W horizons decimal degrees - --------------------------------------------------------------------------- -Files ------ - -Parameter initialization file -Input meteorological data file -Output meteorological data file (*.mtcout) - -Example initialization file: - ----top of init file------------------------------------------------------- -sample (this entire line written to outfiles) -other comments (this entire line discarded) - -IOFILES Keyword, don't edit this line -test.mtcin Name of input meteorological data file -test Prefix for output file - -CONTROL Keyword, don't edit this line -3 (int) Number of header lines in input file -365 (int) Number of days of data in input file -0 (int) Dewpoint temperature input? (0=NO, 1=YES) -1 (int) Humidity output flag (0=VPD, 1=VP) -1 (int) Year field in input file? (0=NO, 1=YES) - -PARAMETERS Keyword, don't edit this line -500.0 (double) Base station elevation, meters -50.0 (double) Base station annual precip isohyet, cm -48.0 (double) Site latitude, degrees (- for south) -1500.0 (double) Site elevation, meters -15.0 (double) Site slope, degrees -180.0 (double) Site aspect, degrees (0=N,90=E,180=S,270=W) -75.0 (double) Site annual precip isohyet, cm -2.0 (double) Site east horizon, degrees -5.0 (double) Site west horizon, degrees --6.0 (double) Lapse rate for max temperature, deg C/1000m --3.0 (double) Lapse rate for min temperature, deg C/1000m - -END Keyword, don't edit this line ----end of init file------------------------------------------------------- - -*.init FILE INFO -For all lines, except the header and comment lines, the parameter value on the -left can be followed by comments on the right, as long as there is whitespace -separating the value on the left from the following comment. Blank lines can be -inserted or deleted, but all keyword lines and parameter lines must be -included. The order and number of non-blank lines in this file is crucial. The -keywords are there as a rudimentary quality check on the format of the -initialization file, and they ensure that the proper number of lines are read. -They DON'T ensure that the parameters are in the proper order. - -NOTE: The dashed lines at the top and bottom of the example file shown above -are NOT part of the file. - - -INPUT FILE INFO -All input temperatures are in degrees Celcius, and precipitation is in -centimeters. -Input data file format: - - -. -. -. - -Example input data file... without year field, and without dewpoint temperature ----start of input data file -------------- -This is a header line, which is discarded -101 16.0 2.0 1.2 -102 17.0 3.0 0.1 -103 16.5 5.2 0.0 -104 20.1 6.4 0.0 ----end of input data file ---------------- - - - -OUTPUT FILE INFO -The output file is created by appending ".mtcout" to the output filename -prefix specified in the initialization file. Existing files are overwritten, -so be careful to rename files between runs if you want to save the results, -and for safety, don't use ".mtcout" as the extension for the input data file. - -****************************** -Input and Output CONTROL FLAGS -****************************** -There is a flag in the initialization file that controls the input of dewpoint -temperature information. If your input file does not contain dewpoint data, -set this flag to 0. Otherwise, if you have dewpoint temperature information -in your input file, set this flag to 1, and be sure that the dewpoint -temperature field is between the tmin and prcp fields, as specified under -the heading "INPUT FILE INFO", above. - -There is another flag in the initialization file that controls the output of -humidity information. When set to 0, humidity output is the vapor pressure -deficit from the saturated vapor pressure at the daylight average temperature -to the ambient vapor pressure, in Pa. When the flag is set to 1, the output is -simply the ambient vapor pressure in Pa. By using the ambient vapor pressure -(flag set to 1), you can avoid possible errors related to a difference between -the temperature chosen for saturation calculations (in this case tday) and the -actual temperature at any given time of day. - -Another flag controls the treatment of a year-field in the input and output -files. When the flag is set (1), it is assumed that the first field in the -input file contains the year, which can be any integer between -32000 and -32000 (approx). This field is simply copied line-for-line into the output -file, where it preceeds the standard yearday output field. - -Example output data file... (using VPD output, no PAR, no year field) ----start of output data file ----------------------------------- -MTCLIM v3.1 OUTPUT FILE : Mon Mar 24 14:50:51 1997 -MTCLIM version 3.1 testing - yday Tmax Tmin Tday prcp VPD srad daylen - (deg C) (deg C) (deg C) (cm) (Pa) (W m-2) (s) - 101 10.00 -0.50 7.39 1.80 439.51 379.95 47672 - 102 11.00 1.10 7.67 0.15 387.27 364.95 47880 - 103 10.50 2.80 6.84 0.00 243.78 302.15 48087 - 104 14.10 3.40 9.29 0.00 390.67 390.22 48293 ----end of output data file ------------------------------------- - - -*/ - -/********************* -** ** -** START OF CODE ** -** ** -*********************/ - -#include -#include -#include -#include -#include - -#include "ini.h" /* file structure and I/O routines */ -#include "mtclim43_constants.h" /* physical constants */ -#include "mtclim43_parameters.h" /* model parameters */ - -/**************************** -** ** -** STRUCTURE DEFINITIONS ** -** ** -****************************/ -typedef struct -{ - file init; /* initialization file */ - file inmet; /* input meteorological data file */ - file outmet; /* output meteorological data file */ - int nhead; /* number of header lines in input met file */ - int ndays; /* number of days of data in input file */ - int indewpt; /* input dewpoint temperature flag (0=NO, 1=YES) */ - int outhum; /* output humidity flag (0=VPD, 1=VP) */ - int inyear; /* input year flag (0=NO, 1=YES) */ - char systime[100]; /* system time at start of simulation */ - char header[200]; /* header string, written to all output files */ - char outprefix[200]; /* output filename prefix */ -} control_struct; - -typedef struct -{ - double base_elev; /* base elevation, meters */ - double base_isoh; /* base annual precip isohyet, cm */ - double site_lat; /* site latitude, dec. degrees (- for south) */ - double site_elev; /* site elevation, meters */ - double site_slp; /* site slope, degrees */ - double site_asp; /* site aspect, degrees */ - double site_isoh; /* site annual precip isohyet, cm */ - double site_ehoriz; /* site east horizon, degrees */ - double site_whoriz; /* site west horizon, degrees */ - double tmax_lr; /* maximum temperature lapse rate, deg C/1000m */ - double tmin_lr; /* minimum temperature lapse rate, deg C/1000m */ -} parameter_struct; - -typedef struct -{ - int *year; /* array of year values */ - int *yday; /* array of yearday values */ - double *tmax; /* array of base maximum temperature values */ - double *tmin; /* array of base minimum temperature values */ - double *prcp; /* array of base daily precipitation values */ - double *tdew; /* array of base dewpoint temperature values */ - double *s_tmax; /* array of site tmax values */ - double *s_tmin; /* array of site tmin values */ - double *s_tday; /* array of site daylight temperature values */ - double *s_prcp; /* array of site prcp values */ - double *s_hum; /* array of site humidity values (VPD or VP, Pa) */ - double *s_srad; /* array of site shortwave radiation values */ - double *s_dayl; /* array of site daylength values */ - double *s_swe; /* array of site snow water equivalent values (cm) */ -} data_struct; - -/******************************** -** ** -** FUNCTION PROTOTYPES ** -** ** -********************************/ -int read_init(control_struct *ctrl, parameter_struct *p); -int data_alloc(const control_struct *ctrl, data_struct *data); -int read_inmet(const control_struct *ctrl, data_struct *data); -int calc_tair(const control_struct *ctrl, const parameter_struct *p, - data_struct *data); -int calc_prcp(const control_struct *ctrl, const parameter_struct *p, - data_struct *data); -int snowpack(const control_struct *ctrl, const parameter_struct *p, - data_struct *data); -int calc_srad_humidity(const control_struct *ctrl, const parameter_struct *p, - data_struct *data); -int calc_srad_humidity_iterative(const control_struct *ctrl, - const parameter_struct *p, data_struct *data); -int write_out(const control_struct *ctrl, const data_struct *data); -int data_free(const control_struct *ctrl, data_struct *data); -double calc_pet(double rad, double ta, double pa, double dayl); -double atm_pres(double elev); -int pulled_boxcar(double *input,double *output,int n,int w,int w_flag); - -/******************************** -** ** -** START OF MAIN() ** -** ** -********************************/ -void mtc(char *argv) -{ - int argc=2; - /* variable declarations */ - control_struct ctrl; /* iofiles and program control variables */ - parameter_struct p; /* site and base parameters */ - data_struct data; /* site and base meteorological data arrays */ - - /* inital progress indicator */ - printf("Starting MTCLIM version 4.3\n"); - - /* read the name of the initialization file from the command line - and store as init.name */ - if (argc != 2) - { - printf("usage: \n"); - exit(1); - } - strcpy(ctrl.init.name, argv); - - /* read initialization file, open input files, do basic error checking - on input parameters, and open output file */ - if (read_init(&ctrl, &p)) - { - printf("Error in read_init()... exiting\n"); - exit(1); - } - printf("Completed read_init()\n"); - - /* allocate space in the data arrays for input and output data */ - if (data_alloc(&ctrl, &data)) - { - printf("Error in data_alloc()... exiting\n"); - exit(1); - } - printf("Completed data_alloc()\n"); - - /* read meteorological data from input file into data arrays */ - if (read_inmet(&ctrl, &data)) - { - printf("Error in read_inmet()... exiting\n"); - exit(1); - } - printf("Completed read_inmet()\n"); - - /* estimate daily air temperatures */ - if (calc_tair(&ctrl, &p, &data)) - { - printf("Error in calc_tair()... exiting\n"); - exit(1); - } - printf("Completed calc_tair()\n"); - - /* estimate daily precipitation */ - if (calc_prcp(&ctrl, &p, &data)) - { - printf("Error in calc_prcp()... exiting\n"); - exit(1); - } - printf("Completed calc_prcp()\n"); - - /* estimate daily snowpack */ - if (snowpack(&ctrl, &p, &data)) - { - printf("Error in snowpack()... exiting\n"); - exit(1); - } - printf("Completed snowpack()\n"); - - /* test for the presence of Tdew observations, and branch to the - appropriate srad and humidity algorithms */ - if (ctrl.indewpt) - { - /* estimate srad and humidity using real Tdew data */ - if (calc_srad_humidity(&ctrl, &p, &data)) - { - printf("Error in calc_srad_humidity()... exiting\n"); - exit(1); - } - printf("Completed calc_srad_humidity()\n"); - - } - else /* no dewpoint temperature data */ - { - /* estimate srad and humidity with iterative algorithm */ - if (calc_srad_humidity_iterative(&ctrl, &p, &data)) - { - printf("Error in calc_srad_humidity_iterative()... exiting\n"); - exit(1); - } - printf("Completed calc_srad_humidity_iterative()\n"); - } - - /* write output file */ - if (write_out(&ctrl, &data)) - { - printf("Error in write_out()... exiting\n"); - exit(1); - } - printf("Completed write_out()\n"); - - /* free previously allocated memory before returning */ - if (data_free(&ctrl, &data)) - { - printf("Error in data_free()... exiting\n"); - exit(1); - } - printf("Completed data_free()\n"); - printf("Output filename: %s\n",ctrl.outmet.name); - -} -/* end of main */ - -/**************************** -** ** -** START OF FUNCTION ** -** DEFINITIONS ** -** ** -****************************/ -/* read_init() reads data out of init file, opens input and output files */ -int read_init(control_struct *ctrl, parameter_struct *p) -{ - struct tm *tm_ptr; - time_t lt; - int ok = 1; - char temp[200]; - char key1[] = "IOFILES"; - char key2[] = "CONTROL"; - char key3[] = "PARAMETERS"; - char key4[] = "END"; - char keyword[80]; - char comment[200]; - - /* open the main init file for ascii read and check for errors */ - if (file_open(&(ctrl->init),'i')) - { - printf("Error opening init file: %s\n",ctrl->init.name); - ok = 0; - } - - /* get the system time */ - if (ok) - { - lt = time(NULL); - tm_ptr = localtime(<); - strcpy(ctrl->systime,asctime(tm_ptr)); - } - - /* read the header line and the comment line from init file */ - if (ok && fgets(ctrl->header, 200, ctrl->init.ptr)==NULL) - { - printf("Error reading header line\n"); - ok = 0; - } - if (ok && fgets(comment, 200, ctrl->init.ptr)==NULL) - { - printf("Error reading comment line\n"); - ok = 0; - } - - /* scan for the IOFILES keyword, exit if not next */ - if (ok && scan_value(ctrl->init, keyword, 's')) - { - printf("Error reading IOFILES keyword\n"); - ok=0; - } - if (ok && strcmp(keyword, key1)) - { - printf("Expecting keyword %s in file %s\n",key1,ctrl->init.name); - ok=0; - } - /* read input data filename and open the input data file for ascii read */ - if (ok && scan_open(ctrl->init, &ctrl->inmet, 'i')) - { - printf("Error opening input file \n"); - ok=0; - } - /* read the output filename prefix */ - if (ok && scan_value(ctrl->init, ctrl->outprefix, 's')) - { - printf("Error reading output file prefix\n"); - ok=0; - } - - /* scan for the CONTROL keyword, exit if not next */ - if (ok && scan_value(ctrl->init, keyword, 's')) - { - printf("Error reading CONTROL keyword\n"); - ok=0; - } - if (ok && strcmp(keyword, key2)) - { - printf("Expecting keyword %s in file %s\n", key2, ctrl->init.name); - ok=0; - } - /* number of header lines in input data file */ - if (ok && scan_value(ctrl->init, &ctrl->nhead, 'i')) - { - printf("Error reading number of header lines\n"); - ok=0; - } - /* total number of days in input data file */ - if (ok && scan_value(ctrl->init, &ctrl->ndays, 'i')) - { - printf("Error reading number of days\n"); - ok=0; - } - /* read the binary flag for dewpoint data in input file */ - if (ok && scan_value(ctrl->init, &ctrl->indewpt, 'i')) - { - printf("Error reading input dewpoint flag\n"); - ok=0; - } - /* read the binary flag for humidity output format - 0 = VPD, 1 = VP */ - if (ok && scan_value(ctrl->init, &ctrl->outhum, 'i')) - { - printf("Error reading output humidity flag\n"); - ok=0; - } - /* read the binary flag for year field in input file */ - if (ok && scan_value(ctrl->init, &ctrl->inyear, 'i')) - { - printf("Error reading year field flag\n"); - ok=0; - } - - /* scan for the PARAMETERS keyword, exit if not next */ - if (ok && scan_value(ctrl->init, keyword, 's')) - { - printf("Error reading PARAMETERS keyword\n"); - ok=0; - } - if (ok && strcmp(keyword, key3)) - { - printf("Expecting keyword %s in file %s\n", key3, ctrl->init.name); - ok=0; - } - /* base station elevation */ - if (ok && scan_value(ctrl->init, &p->base_elev, 'd')) - { - printf("Error reading base elevation\n"); - ok=0; - } - /* base station annual isohyet */ - if (ok && scan_value(ctrl->init, &p->base_isoh, 'd')) - { - printf("Error reading base station isohyet\n"); - ok=0; - } - /* prediction site latitude */ - if (ok && scan_value(ctrl->init, &p->site_lat, 'd')) - { - printf("Error reading prediction site latitude\n"); - ok=0; - } - /* prediction site elevation */ - if (ok && scan_value(ctrl->init, &p->site_elev, 'd')) - { - printf("Error reading prediction site elevation\n"); - ok=0; - } - /* prediction site slope */ - if (ok && scan_value(ctrl->init, &p->site_slp, 'd')) - { - printf("Error reading prediction site slope\n"); - ok=0; - } - /* prediction site aspect */ - if (ok && scan_value(ctrl->init, &p->site_asp, 'd')) - { - printf("Error reading prediction site aspect\n"); - ok=0; - } - /* predictions site annual isohyet */ - if (ok && scan_value(ctrl->init, &p->site_isoh, 'd')) - { - printf("Error reading prediction site isohyet\n"); - ok=0; - } - /* prediction site east horizon elevation angle */ - if (ok && scan_value(ctrl->init, &p->site_ehoriz, 'd')) - { - printf("Error reading prediction site east horizon\n"); - ok=0; - } - /* prediction site west horizon elevation angle */ - if (ok && scan_value(ctrl->init, &p->site_whoriz, 'd')) - { - printf("Error reading prediction site west horizon\n"); - ok=0; - } - /* maximum tempreature lapse rate */ - if (ok && scan_value(ctrl->init, &p->tmax_lr, 'd')) - { - printf("Error reading lapse rate for maximum temperature\n"); - ok=0; - } - /* minimum temperature lapse rate */ - if (ok && scan_value(ctrl->init, &p->tmin_lr, 'd')) - { - printf("Error reading lapse rate for minimum temperature\n"); - ok=0; - } - - /* scan for the END keyword, exit if not next */ - if (ok && scan_value(ctrl->init, keyword, 's')) - { - printf("Error reading END keyword\n"); - ok=0; - } - if (ok && strcmp(keyword, key4)) - { - printf("Expecting keyword %s in file %s\n", key4, ctrl->init.name); - ok=0; - } - - /* if all data read successfully, do basic error checking on parameters */ - if (ok) - { - /* error checking for control parameters */ - /* make sure total number of days is > 0 */ - if (ctrl->ndays <= 0) - { - printf("ERROR: number of days must be > 0 ...\n"); - ok=0; - } - /* make sure number of header lines is > 0 */ - if (ctrl->nhead < 0) - { - printf("ERROR: number of header lines must be >= 0 ...\n"); - ok=0; - } - /* check dewpoint input flag */ - if ((ctrl->indewpt < 0) || (ctrl->indewpt > 1)) - { - printf("WARNING: input dewpoint flag should be 0 or 1: assuming dewpoint input ...\n"); - ctrl->indewpt = 1; - } - /* check humidity output flag */ - if ((ctrl->outhum < 0) || (ctrl->outhum > 1)) - { - printf("WARNING: output humidity flag should be 0 or 1: assuming vapor pressure output ...\n"); - ctrl->outhum = 1; - } - if ((ctrl->inyear < 0) || (ctrl->inyear > 1)) - { - printf("WARNING: input year field flag should be 0 or 1 ...\n"); - ok=0; - } - - /* error checking for parameters */ - if (p->base_elev > 5000) - { - printf("WARNING: base station elev = %.1lf m: be sure to use meters, not feet\n",p->base_elev); - } - if (p->base_isoh <= 0.0) - { - printf("ERROR: base station isohyet must be > 0.0 ...\n"); - ok=0; - } - if ((p->site_lat < -90.0) || (p->site_lat > 90.0)) - { - printf("ERROR: site latitude must be in the range -90.0 - 90.0 degrees ...\n"); - ok=0; - } - if (p->site_elev > 5000) - { - printf("WARNING: site elev = %.1lf m: be sure to use meters, not feet\n",p->site_elev); - } - if (p->site_slp > 60.0) - { - printf("WARNING: site slope = %.1lf deg: be sure to use deg, not %%\n",p->site_slp); - } - if (p->site_slp < 0) - { - printf("ERROR: site slope must be >= 0.0\n"); - ok=0; - } - if ((p->site_asp < 0.0) || (p->site_asp > 360.0)) - { - printf("ERROR: site aspect must be in the range 0.0 - 360.0 degrees ...\n"); - ok=0; - } - if (p->site_isoh < 0.0) - { - printf("ERROR: site isohyet must be >= 0.0 ...\n"); - ok=0; - } - if ((p->site_ehoriz < 0.0) || (p->site_ehoriz > 180.0)) - { - printf("ERROR: site east horizon must be in the range 0.0 - 180.0 degrees ...\n"); - ok=0; - } - if ((p->site_whoriz < 0.0) || (p->site_whoriz > 180.0)) - { - printf("ERROR: site west horizon must be in the range 0.0 - 180.0 degrees ...\n"); - ok=0; - } - if (p->tmax_lr > 0.0) - { - printf("WARNING: Tmax lapse rate = %.2lf : This should typically be negative ...\n"); - } - if (p->tmin_lr > 0.0) - { - printf("WARNING: Tmin lapse rate = %.2lf : This should typically be negative ...\n"); - } - } /* end of basic error checking routine */ - - /* generate output filename based on outprefix and open file for ascii - write. Also write date, time, and header information into output file */ - if (ok) - { - strcpy(temp,ctrl->outprefix); - strcat(temp,POSTFIX); - strcpy(ctrl->outmet.name,temp); - if (file_open(&ctrl->outmet, 'o')) - { - printf("Error opening output file\n"); - ok=0; - } - } - if (ok) - { - fprintf(ctrl->outmet.ptr,"%s", ctrl->header); - fprintf(ctrl->outmet.ptr,"MTCLIM v4.3 OUTPUT FILE : %s", ctrl->systime); - } - - return (!ok); -} -/* end of read_init() */ - -/* data_alloc() allocates space for input and output data arrays */ -int data_alloc(const control_struct *ctrl, data_struct *data) -{ - int ok=1; - int ndays; - - ndays = ctrl->ndays; - - if (ok && ctrl->inyear && !(data->year = (int*) malloc(ndays * sizeof(int)))) - { - printf("Error allocating for year array\n"); - ok=0; - } - if (ok && !(data->yday = (int*) malloc(ndays * sizeof(int)))) - { - printf("Error allocating for yearday array\n"); - ok=0; - } - if (ok && !(data->tmax = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for tmax array\n"); - ok=0; - } - if (ok && !(data->tmin = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for tmin array\n"); - ok=0; - } - if (ok && !(data->prcp = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for prcp array\n"); - ok=0; - } - if (ok && ctrl->indewpt && !(data->tdew = (double*) malloc(ndays * - sizeof(double)))) - { - printf("Error allocating for input humidity array\n"); - ok=0; - } - if (ok && !(data->s_tmax = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for site Tmax array\n"); - ok=0; - } - if (ok && !(data->s_tmin = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for site Tmin array\n"); - ok=0; - } - if (ok && !(data->s_tday = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for site Tday array\n"); - ok=0; - } - if (ok && !(data->s_prcp = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for site prcp array\n"); - ok=0; - } - if (ok && !(data->s_hum = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for site VPD array\n"); - ok=0; - } - if (ok && !(data->s_srad = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for site radiation array\n"); - ok=0; - } - if (ok && !(data->s_dayl = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for site daylength array\n"); - ok=0; - } - if (ok && !(data->s_swe = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for site snowpack array\n"); - ok=0; - } - - return (!ok); -} -/* end of data_alloc() */ - -/* read_inmet() reads data from the meteorological input file into arrays */ -int read_inmet(const control_struct *ctrl, data_struct *data) -{ - int ok=1; - - int i,ndays; - int id,iy; - - char header[200]; - - id = ctrl->indewpt; - iy = ctrl->inyear; - - ndays = ctrl->ndays; - - /* read header lines and discard */ - for (i=0 ; inhead ; i++) - { - if (ok && fgets(header, 200, ctrl->inmet.ptr)==NULL) - { - printf("Error reading input header line\n"); - ok = 0; - } - } - - /* read lines of data, variable read format depending on flag */ - i = 0; - while (ok && iinmet.ptr, "%d",&(data->year[i])) != 1) - { - printf("Error reading year from input file line #%d \n", i+1); - ok=0; - } - } - if (id) /* dewpoint data in input file */ - { - if (fscanf(ctrl->inmet.ptr, "%d%lf%lf%lf%lf", - &(data->yday[i]), &(data->tmax[i]), &(data->tmin[i]), - &(data->tdew[i]), &(data->prcp[i])) != 5) - { - printf("Error reading line #%d from input file\n", i+1); - ok=0; - } - } - else /* no dewpoint data in input file */ - { - if (fscanf(ctrl->inmet.ptr, "%d%lf%lf%lf", - &(data->yday[i]), &(data->tmax[i]), &(data->tmin[i]), - &(data->prcp[i])) != 4) - { - printf("Error reading line #%d from input file\n", i+1); - ok=0; - } - } - - /* check for valid yearday */ - if (ok && ((data->yday[i] > 366) || data->yday[i] < 1)) - { - printf("Invalid yearday (%d) at line #%d\n",data->yday[i],i+1); - ok=0; - } - - i++; - - } /* end while */ - - /* check for gaps in yearday field */ - i = 1; - while (ok && iyday[i] != 1) /* not January 1 */ - { - if (data->yday[i] != (data->yday[i-1]+1)) - { - printf("Gap in yeardays at line #%d\n",i+1); - ok=0; - } - } - else /* January 1 */ - { - if ((data->yday[i-1] != 365) && (data->yday[i-1] != 366)) - { - printf("Gap in yeardays at line #%d\n",i+1); - ok=0; - } - } - - i++; - } - - /* check for tmax < tmin */ - i=0; - while (ok && itmax[i] < data->tmin[i]) - { - printf("ERROR: Tmax < Tmin on line #%d\n",i+1+ctrl->nhead); - ok=0; - } - i++; - } - - return (!ok); -} -/* end of read_inmet() */ - -/* calc_tair() calculates daily air temperatures */ -int calc_tair(const control_struct *ctrl, const parameter_struct *p, - data_struct *data) -{ - int ok=1; - int i,ndays; - double dz; - double tmean, tmax, tmin; - - ndays = ctrl->ndays; - /* calculate elevation difference in kilometers */ - dz = (p->site_elev - p->base_elev)/1000.0; - - /* apply lapse rate corrections to tmax and tmin */ - /* Since tmax lapse rate usually has a larger absolute value than tmin - lapse rate, it is possible at high elevation sites for these corrections - to result in tmin > tmax. Check for that occurrence and force - tmin = corrected tmax - 0.5 deg C. */ - for (i=0 ; is_tmax[i] = tmax = data->tmax[i] + (dz * p->tmax_lr); - data->s_tmin[i] = tmin = data->tmin[i] + (dz * p->tmin_lr); - - /* derived temperatures */ - tmean = (tmax + tmin)/2.0; - data->s_tday[i] = ((tmax - tmean)*TDAYCOEF) + tmean; - } - - return (!ok); -} -/* end of calc_tair() */ - -/* calc_prcp() calculates daily total precipitation */ -int calc_prcp(const control_struct *ctrl, const parameter_struct *p, - data_struct *data) -{ - int ok=1; - int i,ndays; - double ratio; - - ndays = ctrl->ndays; - - ratio = p->site_isoh / p->base_isoh; - - if (ok) - { - for (i=0 ; is_prcp[i] = data->prcp[i] * ratio; - } - } - - return (!ok); -} -/* end of calc_prcp() */ - -/* snowpack() estimates the accumulation and melt of snow for radiation -algorithm corrections */ -int snowpack(const control_struct *ctrl, const parameter_struct *p, - data_struct *data) -{ - int ok=1; - int i,ndays,count; - int start_yday,prev_yday; - double snowpack,newsnow,snowmelt,sum; - - ndays = ctrl->ndays; - - /* first pass to initialize SWE array */ - snowpack = 0.0; - for (i=0 ; is_tmin[i] <= SNOW_TCRIT) newsnow = data->s_prcp[i]; - else snowmelt = SNOW_TRATE * (data->s_tmin[i] - SNOW_TCRIT); - snowpack += newsnow - snowmelt; - if (snowpack < 0.0) snowpack = 0.0; - data->s_swe[i] = snowpack; - } - - /* use the first pass to set the initial snowpack conditions for the - first day of data */ - start_yday = data->yday[0]; - if (start_yday == 1) prev_yday = 365; - else prev_yday = start_yday-1; - count = 0; - sum = 0.0; - for (i=1 ; iyday[i] == start_yday || data->yday[i] == prev_yday) - { - count ++; - sum += data->s_swe[i]; - } - } - /* Proceed with correction if there are valid days to reinitialize - the snowpack estiamtes. Otherwise use the first-pass estimate. */ - if (count) - { - snowpack = sum/(double)count; - for (i=0 ; is_tmin[i] <= SNOW_TCRIT) newsnow = data->s_prcp[i]; - else snowmelt = SNOW_TRATE * (data->s_tmin[i] - SNOW_TCRIT); - snowpack += newsnow - snowmelt; - if (snowpack < 0.0) snowpack = 0.0; - data->s_swe[i] = snowpack; - } - } - - return (!ok); -} -/* end of snowpack() */ - -/* when dewpoint temperature observations are available, radiation and -humidity can be estimated directly */ -int calc_srad_humidity(const control_struct *ctrl, const parameter_struct *p, - data_struct *data) -{ - int ok=1; - int i,ndays; - double pva,pvs,vpd; - int ami,yday; - double ttmax0[366]; - double flat_potrad[366]; - double slope_potrad[366]; - double daylength[366]; - double *dtr, *sm_dtr; - double tmax,tmin; - double t1,t2; - double pratio; - double lat,coslat,sinlat,dt,dh,h; - double cosslp,sinslp,cosasp,sinasp; - double bsg1,bsg2,bsg3; - double decl,cosdecl,sindecl,cosegeom,sinegeom,coshss,hss; - double sc,dir_beam_topa; - double sum_flat_potrad, sum_slope_potrad, sum_trans; - double cosh,sinh; - double cza,cbsa,coszeh,coszwh; - double dir_flat_topa,am; - double trans1,trans2; - double t_tmax,b,t_fmax; - double t_final,pdif,pdir,srad1,srad2; - double sky_prop; - double avg_horizon, slope_excess; - double horizon_scalar, slope_scalar; - - /* optical airmass by degrees */ - double optam[21] = {2.90,3.05,3.21,3.39,3.69,3.82,4.07,4.37,4.72,5.12,5.60, - 6.18,6.88,7.77,8.90,10.39,12.44,15.36,19.79,26.96,30.00}; - - /* number of simulation days */ - ndays = ctrl->ndays; - - /* calculate humidity from Tdew observations */ - for (i=0 ; itdew[i] / (239.0 + data->tdew[i])); - if (ctrl->outhum) - { - /* output humidity as vapor pressure */ - data->s_hum[i] = pva; - } - else - { - /* output humidity as vapor pressure deficit */ - /* calculate saturation vapor pressure at tday */ - pvs = 610.7 * exp(17.38 * data->s_tday[i] / (239.0 + data->s_tday[i])); - /* calculate vpd */ - vpd = pvs-pva; - if (vpd < 0.0) vpd = 0.0; - data->s_hum[i] = vpd; - } - } - - /* estimate radiation using Tdew observations */ - /* allocate space for DTR and smoothed DTR arrays */ - if (!(dtr = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for DTR array\n"); - ok=0; - } - if (!(sm_dtr = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for smoothed DTR array\n"); - ok=0; - } - - /* calculate diurnal temperature range for transmittance calculations */ - for (i=0 ; itmax[i]; - tmin = data->tmin[i]; - if (tmax < tmin) tmax = tmin; - dtr[i] = tmax-tmin; - } - - /* smooth dtr array using a 30-day antecedent smoothing window */ - if (ndays >= 30) - { - if (pulled_boxcar(dtr, sm_dtr, ndays, 30, 0)) - { - printf("Error in boxcar smoothing, calc_srad_humidity()\n"); - ok=0; - } - } - else /* smoothing window width = ndays */ - { - if (pulled_boxcar(dtr, sm_dtr, ndays, ndays, 0)) - { - printf("Error in boxcar smoothing, calc_srad_humidity()\n"); - ok=0; - } - } - - /***************************************** - * * - * start of the main radiation algorithm * - * * - *****************************************/ - - /* STEP (1) calculate pressure ratio (site/reference) = f(elevation) */ - t1 = 1.0 - (LR_STD * p->site_elev)/T_STD; - t2 = G_STD / (LR_STD * (R/MA)); - pratio = pow(t1,t2); - - /* STEP (2) correct initial transmittance for elevation */ - trans1 = pow(TBASE,pratio); - - /* STEP (3) build 366-day array of ttmax0, potential rad, and daylength */ - - /* precalculate the transcendentals */ - lat = p->site_lat; - /* check for (+/-) 90 degrees latitude, throws off daylength calc */ - lat *= RADPERDEG; - if (lat > 1.5707) lat = 1.5707; - if (lat < -1.5707) lat = -1.5707; - coslat = cos(lat); - sinlat = sin(lat); - cosslp = cos(p->site_slp * RADPERDEG); - sinslp = sin(p->site_slp * RADPERDEG); - cosasp = cos(p->site_asp * RADPERDEG); - sinasp = sin(p->site_asp * RADPERDEG); - /* cosine of zenith angle for east and west horizons */ - coszeh = cos(1.570796 - (p->site_ehoriz * RADPERDEG)); - coszwh = cos(1.570796 - (p->site_whoriz * RADPERDEG)); - - /* sub-daily time and angular increment information */ - dt = SRADDT; /* set timestep */ - dh = dt / SECPERRAD; /* calculate hour-angle step */ - - /* begin loop through yeardays */ - for (i=0 ; i<365 ; i++) - { - /* calculate cos and sin of declination */ - decl = MINDECL * cos(((double)i + DAYSOFF) * RADPERDAY); - cosdecl = cos(decl); - sindecl = sin(decl); - - /* do some precalculations for beam-slope geometry (bsg) */ - bsg1 = -sinslp * sinasp * cosdecl; - bsg2 = (-cosasp * sinslp * sinlat + cosslp * coslat) * cosdecl; - bsg3 = (cosasp * sinslp * coslat + cosslp * sinlat) * sindecl; - - /* calculate daylength as a function of lat and decl */ - cosegeom = coslat * cosdecl; - sinegeom = sinlat * sindecl; - coshss = -(sinegeom) / cosegeom; - if (coshss < -1.0) coshss = -1.0; /* 24-hr daylight */ - if (coshss > 1.0) coshss = 1.0; /* 0-hr daylight */ - hss = acos(coshss); /* hour angle at sunset (radians) */ - /* daylength (seconds) */ - daylength[i] = 2.0 * hss * SECPERRAD; - - /* solar constant as a function of yearday (W/m^2) */ - sc = 1368.0 + 45.5*sin((2.0*PI*(double)i/365.25) + 1.7); - /* extraterrestrial radiation perpendicular to beam, total over - the timestep (J) */ - dir_beam_topa = sc * dt; - - sum_trans = 0.0; - sum_flat_potrad = 0.0; - sum_slope_potrad = 0.0; - - /* begin sub-daily hour-angle loop, from -hss to hss */ - for (h=-hss ; h 0.0) - { - /* when sun is above the ideal (flat) horizon, do all the - flat-surface calculations to determine daily total - transmittance, and save flat-surface potential radiation - for later calculations of diffuse radiation */ - - /* potential radiation for this time period, flat surface, - top of atmosphere */ - dir_flat_topa = dir_beam_topa * cza; - - /* determine optical air mass */ - am = 1.0/(cza + 0.0000001); - if (am > 2.9) - { - ami = (int)(acos(cza)/RADPERDEG) - 69; - if (ami < 0) ami = 0; - if (ami > 20) ami = 20; - am = optam[ami]; - } - - /* correct instantaneous transmittance for this optical - air mass */ - trans2 = pow(trans1,am); - - /* instantaneous transmittance is weighted by potential - radiation for flat surface at top of atmosphere to get - daily total transmittance */ - sum_trans += trans2 * dir_flat_topa; - - /* keep track of total potential radiation on a flat - surface for ideal horizons */ - sum_flat_potrad += dir_flat_topa; - - /* keep track of whether this time step contributes to - component 1 (direct on slope) */ - if ((h<0.0 && cza>coszeh && cbsa>0.0) || - (h>=0.0 && cza>coszwh && cbsa>0.0)) - { - /* sun between east and west horizons, and direct on - slope. this period contributes to component 1 */ - sum_slope_potrad += dir_beam_topa * cbsa; - } - - } /* end if sun above ideal horizon */ - - } /* end of sub-daily hour-angle loop */ - - /* calculate maximum daily total transmittance and daylight average - flux density for a flat surface and the slope */ - if (daylength[i]) - { - ttmax0[i] = sum_trans / sum_flat_potrad; - flat_potrad[i] = sum_flat_potrad / daylength[i]; - slope_potrad[i] = sum_slope_potrad / daylength[i]; - } - else - { - ttmax0[i] = 0.0; - flat_potrad[i] = 0.0; - slope_potrad[i] = 0.0; - } - } /* end of i=365 days loop */ - - /* force yearday 366 = yearday 365 */ - ttmax0[365] = ttmax0[364]; - flat_potrad[365] = flat_potrad[364]; - slope_potrad[365] = slope_potrad[364]; - daylength[365] = daylength[364]; - - /* STEP (4) calculate the sky proportion for diffuse radiation */ - /* uses the product of spherical cap defined by average horizon angle - and the great-circle truncation of a hemisphere. this factor does not - vary by yearday. */ - avg_horizon = (p->site_ehoriz + p->site_whoriz)/2.0; - horizon_scalar = 1.0 - sin(avg_horizon * RADPERDEG); - if (p->site_slp > avg_horizon) slope_excess = p->site_slp - avg_horizon; - else slope_excess = 0.0; - if (2.0*avg_horizon > 180.0) slope_scalar = 0.0; - else - { - slope_scalar = 1.0 - (slope_excess/(180.0 - 2.0*avg_horizon)); - if (slope_scalar < 0.0) slope_scalar = 0.0; - } - sky_prop = horizon_scalar * slope_scalar; - - /* STEP (5) final calculation of daily total radiation */ - for (i=0 ; iyday[i]-1; - pva = 610.7 * exp(17.38 * data->tdew[i] / (239.0 + data->tdew[i])); - t_tmax = ttmax0[yday] + ABASE * pva; - - /* b parameter from 30-day average of DTR */ - b = B0 + B1 * exp(-B2 * sm_dtr[i]); - - /* proportion of daily maximum transmittance */ - t_fmax = 1.0 - 0.9 * exp(-b * pow(dtr[i],C)); - - /* correct for precipitation if this is a rain day */ - if (data->prcp[i]) t_fmax *= RAIN_SCALAR; - - /* final daily total transmittance */ - t_final = t_tmax * t_fmax; - - /* estimate fraction of radiation that is diffuse, on an - instantaneous basis, from relationship with daily total - transmittance in Jones (Plants and Microclimate, 1992) - Fig 2.8, p. 25, and Gates (Biophysical Ecology, 1980) - Fig 6.14, p. 122. */ - pdif = -1.25*t_final + 1.25; - if (pdif > 1.0) pdif = 1.0; - if (pdif < 0.0) pdif = 0.0; - - /* estimate fraction of radiation that is direct, on an - instantaneous basis */ - pdir = 1.0 - pdif; - - /* the daily total radiation is estimated as the sum of the - following two components: - 1. The direct radiation arriving during the part of - the day when there is direct beam on the slope. - 2. The diffuse radiation arriving over the entire daylength - (when sun is above ideal horizon). - */ - - /* component 1 (direct) */ - srad1 = slope_potrad[yday] * t_final * pdir; - - /* component 2 (diffuse) */ - /* includes the effect of surface albedo in raising the diffuse - radiation for obstructed horizons */ - srad2 = flat_potrad[yday] * t_final * pdif * - (sky_prop + DIF_ALB*(1.0-sky_prop)); - - /* snow pack influence on radiation */ - if (data->s_swe[i] > 0.0) - { - /* snow correction in J/m2/day */ - sc = (1.32 + 0.096 * data->s_swe[i]) * 1e6; - /* convert to W/m2 and check for zero daylength */ - if (daylength[yday] > 0.0) sc /= daylength[yday]; - else sc = 0.0; - /* set a maximum correction of 100 W/m2 */ - if (sc > 100.0) sc = 100.0; - } - else sc = 0.0; - - /* save daily radiation and daylength */ - data->s_srad[i] = srad1 + srad2 + sc; - data->s_dayl[i] = daylength[yday]; - } - - /* free local array memory */ - free(dtr); - free(sm_dtr); - - return (!ok); -} /* end of calc_srad_humidity() */ - -/* without Tdew input data, an iterative estimation of shortwave radiation -and humidity is required */ -int calc_srad_humidity_iterative(const control_struct *ctrl, - const parameter_struct *p, data_struct *data) -{ - int ok=1; - int i,j,ndays; - int start_yday,end_yday,isloop; - int ami,yday; - double ttmax0[366]; - double flat_potrad[366]; - double slope_potrad[366]; - double daylength[366]; - double *dtr, *sm_dtr; - double *parray, *window, *t_fmax, *tdew; - double *save_pet; - double sum_prcp,ann_prcp,effann_prcp; - double sum_pet,ann_pet; - double tmax,tmin; - double t1,t2; - double pratio; - double lat,coslat,sinlat,dt,h,dh; - double cosslp,sinslp,cosasp,sinasp; - double bsg1,bsg2,bsg3; - double decl,cosdecl,sindecl,cosegeom,sinegeom,coshss,hss; - double sc,dir_beam_topa; - double sum_flat_potrad,sum_slope_potrad,sum_trans; - double cosh,sinh; - double cza,cbsa,coszeh,coszwh; - double dir_flat_topa,am; - double pva,t_tmax,b; - double tmink,pet,ratio,ratio2,ratio3,tdewk; - double pvs,vpd; - double trans1,trans2; - double t_final,pdif,pdir,srad1,srad2; - double pa; - double sky_prop; - double avg_horizon, slope_excess; - double horizon_scalar, slope_scalar; - - /* optical airmass by degrees */ - double optam[21] = {2.90,3.05,3.21,3.39,3.69,3.82,4.07,4.37,4.72,5.12,5.60, - 6.18,6.88,7.77,8.90,10.39,12.44,15.36,19.79,26.96,30.00}; - - /* number of simulation days */ - ndays = ctrl->ndays; - - /* local array memory allocation */ - /* allocate space for DTR and smoothed DTR arrays */ - if (!(dtr = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for DTR array\n"); - ok=0; - } - if (!(sm_dtr = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for smoothed DTR array\n"); - ok=0; - } - /* allocate space for effective annual precip array */ - if (!(parray = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for effective annual precip array\n"); - ok=0; - } - /* allocate space for the prcp totaling array */ - if (!(window = (double*) malloc((ndays+90)*sizeof(double)))) - { - printf("Error allocating for prcp totaling array\n"); - ok = 0; - } - /* allocate space for t_fmax */ - if (!(t_fmax = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for p_tt_max array\n"); - ok=0; - } - /* allocate space for Tdew array */ - if (!(tdew = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for Tdew array\n"); - ok=0; - } - /* allocate space for save_pet array */ - if (!(save_pet = (double*) malloc(ndays * sizeof(double)))) - { - printf("Error allocating for save_pet array\n"); - ok=0; - } - - /* calculate diurnal temperature range for transmittance calculations */ - for (i=0 ; itmax[i]; - tmin = data->tmin[i]; - if (tmax < tmin) tmax = tmin; - dtr[i] = tmax-tmin; - } - - /* smooth dtr array: After Bristow and Campbell, 1984 */ - if (ndays >= 30) /* use 30-day antecedent smoothing window */ - { - if (pulled_boxcar(dtr, sm_dtr, ndays, 30, 0)) - { - printf("Error in boxcar smoothing, calc_srad_humidity()\n"); - ok=0; - } - } - else /* smoothing window width = ndays */ - { - if (pulled_boxcar(dtr, sm_dtr, ndays, ndays, 0)) - { - printf("Error in boxcar smoothing, calc_srad_humidity()\n"); - ok=0; - } - } - - /* calculate the annual total precip for decision between - simple and arid-corrected humidity algorithm */ - sum_prcp = 0.0; - for (i=0 ; is_prcp[i]; - } - ann_prcp = (sum_prcp/(double)ndays) * 365.25; - if (ann_prcp == 0.0) ann_prcp = 1.0; - - /* Generate the effective annual precip, based on a 3-month - moving-window. Requires some special case handling for the - beginning of the record and for short records. */ - /* check if there are at least 90 days in this input file, if not, - use a simple total scaled to effective annual precip */ - if (ndays < 90) - { - sum_prcp = 0.0; - for (i=0 ; is_prcp[i]; - } - effann_prcp = (sum_prcp/(double)ndays) * 365.25; - /* if the effective annual precip for this period - is less than 8 cm, set the effective annual precip to 8 cm - to reflect an arid condition, while avoiding possible - division-by-zero errors and very large ratios (PET/Pann) */ - if (effann_prcp < 8.0) effann_prcp = 8.0; - for (i=0 ; iyday[0]; - end_yday = data->yday[ndays-1]; - if (start_yday != 1) - { - isloop = (end_yday == start_yday-1) ? 1 : 0; - } - else - { - isloop = (end_yday == 365 || end_yday == 366) ? 1 : 0; - } - - /* fill the first 90 days of window */ - for (i=0 ; i<90 ; i++) - { - if (isloop) window[i] = data->s_prcp[ndays-90+i]; - else window[i] = data->s_prcp[i]; - } - /* fill the rest of the window array */ - for (i=0 ; is_prcp[i]; - } - - /* for each day, calculate the effective annual precip from - scaled 90-day total */ - for (i=0 ; i= 90 */ - - /***************************************** - * * - * start of the main radiation algorithm * - * * - *****************************************/ - - /* before starting the iterative algorithm between humidity and - radiation, calculate all the variables that don't depend on - humidity so they only get done once. */ - - /* STEP (1) calculate pressure ratio (site/reference) = f(elevation) */ - t1 = 1.0 - (LR_STD * p->site_elev)/T_STD; - t2 = G_STD / (LR_STD * (R/MA)); - pratio = pow(t1,t2); - - /* STEP (2) correct initial transmittance for elevation */ - trans1 = pow(TBASE,pratio); - - /* STEP (3) build 366-day array of ttmax0, potential rad, and daylength */ - - /* precalculate the transcendentals */ - lat = p->site_lat; - /* check for (+/-) 90 degrees latitude, throws off daylength calc */ - lat *= RADPERDEG; - if (lat > 1.5707) lat = 1.5707; - if (lat < -1.5707) lat = -1.5707; - coslat = cos(lat); - sinlat = sin(lat); - cosslp = cos(p->site_slp * RADPERDEG); - sinslp = sin(p->site_slp * RADPERDEG); - cosasp = cos(p->site_asp * RADPERDEG); - sinasp = sin(p->site_asp * RADPERDEG); - /* cosine of zenith angle for east and west horizons */ - coszeh = cos(1.570796 - (p->site_ehoriz * RADPERDEG)); - coszwh = cos(1.570796 - (p->site_whoriz * RADPERDEG)); - - /* sub-daily time and angular increment information */ - dt = SRADDT; /* set timestep */ - dh = dt / SECPERRAD; /* calculate hour-angle step */ - - /* begin loop through yeardays */ - for (i=0 ; i<365 ; i++) - { - /* calculate cos and sin of declination */ - decl = MINDECL * cos(((double)i + DAYSOFF) * RADPERDAY); - cosdecl = cos(decl); - sindecl = sin(decl); - - /* do some precalculations for beam-slope geometry (bsg) */ - bsg1 = -sinslp * sinasp * cosdecl; - bsg2 = (-cosasp * sinslp * sinlat + cosslp * coslat) * cosdecl; - bsg3 = (cosasp * sinslp * coslat + cosslp * sinlat) * sindecl; - - /* calculate daylength as a function of lat and decl */ - cosegeom = coslat * cosdecl; - sinegeom = sinlat * sindecl; - coshss = -(sinegeom) / cosegeom; - if (coshss < -1.0) coshss = -1.0; /* 24-hr daylight */ - if (coshss > 1.0) coshss = 1.0; /* 0-hr daylight */ - hss = acos(coshss); /* hour angle at sunset (radians) */ - /* daylength (seconds) */ - daylength[i] = 2.0 * hss * SECPERRAD; - - /* solar constant as a function of yearday (W/m^2) */ - sc = 1368.0 + 45.5*sin((2.0*PI*(double)i/365.25) + 1.7); - /* extraterrestrial radiation perpendicular to beam, total over - the timestep (J) */ - dir_beam_topa = sc * dt; - - sum_trans = 0.0; - sum_flat_potrad = 0.0; - sum_slope_potrad = 0.0; - - /* begin sub-daily hour-angle loop, from -hss to hss */ - for (h=-hss ; h 0.0) - { - /* when sun is above the ideal (flat) horizon, do all the - flat-surface calculations to determine daily total - transmittance, and save flat-surface potential radiation - for later calculations of diffuse radiation */ - - /* potential radiation for this time period, flat surface, - top of atmosphere */ - dir_flat_topa = dir_beam_topa * cza; - - /* determine optical air mass */ - am = 1.0/(cza + 0.0000001); - if (am > 2.9) - { - ami = (int)(acos(cza)/RADPERDEG) - 69; - if (ami < 0) ami = 0; - if (ami > 20) ami = 20; - am = optam[ami]; - } - - /* correct instantaneous transmittance for this optical - air mass */ - trans2 = pow(trans1,am); - - /* instantaneous transmittance is weighted by potential - radiation for flat surface at top of atmosphere to get - daily total transmittance */ - sum_trans += trans2 * dir_flat_topa; - - /* keep track of total potential radiation on a flat - surface for ideal horizons */ - sum_flat_potrad += dir_flat_topa; - - /* keep track of whether this time step contributes to - component 1 (direct on slope) */ - if ((h<0.0 && cza>coszeh && cbsa>0.0) || - (h>=0.0 && cza>coszwh && cbsa>0.0)) - { - /* sun between east and west horizons, and direct on - slope. this period contributes to component 1 */ - sum_slope_potrad += dir_beam_topa * cbsa; - } - - } /* end if sun above ideal horizon */ - - } /* end of sub-daily hour-angle loop */ - - /* calculate maximum daily total transmittance and daylight average - flux density for a flat surface and the slope */ - if (daylength[i]) - { - ttmax0[i] = sum_trans / sum_flat_potrad; - flat_potrad[i] = sum_flat_potrad / daylength[i]; - slope_potrad[i] = sum_slope_potrad / daylength[i]; - } - else - { - ttmax0[i] = 0.0; - flat_potrad[i] = 0.0; - slope_potrad[i] = 0.0; - } - } /* end of i=365 days loop */ - - /* force yearday 366 = yearday 365 */ - ttmax0[365] = ttmax0[364]; - flat_potrad[365] = flat_potrad[364]; - slope_potrad[365] = slope_potrad[364]; - daylength[365] = daylength[364]; - - /* STEP (4) calculate the sky proportion for diffuse radiation */ - /* uses the product of spherical cap defined by average horizon angle - and the great-circle truncation of a hemisphere. this factor does not - vary by yearday. */ - avg_horizon = (p->site_ehoriz + p->site_whoriz)/2.0; - horizon_scalar = 1.0 - sin(avg_horizon * RADPERDEG); - if (p->site_slp > avg_horizon) slope_excess = p->site_slp - avg_horizon; - else slope_excess = 0.0; - if (2.0*avg_horizon > 180.0) slope_scalar = 0.0; - else - { - slope_scalar = 1.0 - (slope_excess/(180.0 - 2.0*avg_horizon)); - if (slope_scalar < 0.0) slope_scalar = 0.0; - } - sky_prop = horizon_scalar * slope_scalar; - - /* b parameter, and t_fmax not varying with Tdew, so these can be - calculated once, outside the iteration between radiation and humidity - estimates. Requires storing t_fmax in an array. */ - for (i=0 ; iprcp[i]) t_fmax[i] *= RAIN_SCALAR; - - } - - /* As a first approximation, calculate radiation assuming - that Tdew = Tmin */ - for (i=0 ; iyday[i]-1; - tdew[i] = data->s_tmin[i]; - pva = 610.7 * exp(17.38 * tdew[i] / (239.0 + tdew[i])); - t_tmax = ttmax0[yday] + ABASE * pva; - - /* final daily total transmittance */ - t_final = t_tmax * t_fmax[i]; - - /* estimate fraction of radiation that is diffuse, on an - instantaneous basis, from relationship with daily total - transmittance in Jones (Plants and Microclimate, 1992) - Fig 2.8, p. 25, and Gates (Biophysical Ecology, 1980) - Fig 6.14, p. 122. */ - pdif = -1.25*t_final + 1.25; - if (pdif > 1.0) pdif = 1.0; - if (pdif < 0.0) pdif = 0.0; - - /* estimate fraction of radiation that is direct, on an - instantaneous basis */ - pdir = 1.0 - pdif; - - /* the daily total radiation is estimated as the sum of the - following two components: - 1. The direct radiation arriving during the part of - the day when there is direct beam on the slope. - 2. The diffuse radiation arriving over the entire daylength - (when sun is above ideal horizon). - */ - - /* component 1 */ - srad1 = slope_potrad[yday] * t_final * pdir; - - /* component 2 (diffuse) */ - /* includes the effect of surface albedo in raising the diffuse - radiation for obstructed horizons */ - srad2 = flat_potrad[yday] * t_final * pdif * - (sky_prop + DIF_ALB*(1.0-sky_prop)); - - /* snow pack influence on radiation */ - if (data->s_swe[i] > 0.0) - { - /* snow correction in J/m2/day */ - sc = (1.32 + 0.096 * data->s_swe[i]) * 1e6; - /* convert to W/m2 and check for zero daylength */ - if (daylength[yday] > 0.0) sc /= daylength[yday]; - else sc = 0.0; - /* set a maximum correction of 100 W/m2 */ - if (sc > 100.0) sc = 100.0; - } - else sc = 0.0; - - /* save daily radiation and daylength */ - data->s_srad[i] = srad1 + srad2 + sc; - data->s_dayl[i] = daylength[yday]; - } - - /* estimate annual PET first, to decide which humidity algorithm - should be used */ - /* estimate air pressure at site */ - pa = atm_pres(p->site_elev); - sum_pet = 0.0; - for (i=0 ; is_srad[i],data->s_tday[i],pa,data->s_dayl[i]); - sum_pet += save_pet[i]; - } - ann_pet = (sum_pet/(double)ndays) * 365.25; - - /* humidity algorithm decision: - PET/prcp >= 2.5 -> arid correction - PET/prcp < 2.5 -> use tdew-tmin, which is already finished */ - printf("PET/PRCP = %.4lf\n",ann_pet/ann_prcp); - - if (ann_pet/ann_prcp >= 2.5) - { - printf("Using arid-climate humidity algorithm\n"); - /* Estimate Tdew using the initial estimate of radiation for PET */ - for (i=0 ; is_tmin[i] + 273.15; - pet = save_pet[i]; - - /* calculate ratio (PET/effann_prcp) and correct the dewpoint */ - ratio = pet/parray[i]; - ratio2 = ratio*ratio; - ratio3 = ratio2*ratio; - tdewk = tmink*(-0.127 + 1.121*(1.003 - 1.444*ratio + 12.312*ratio2 - - 32.766*ratio3) + 0.0006*(dtr[i])); - tdew[i] = tdewk - 273.15; - } - - /* Revise estimate of radiation using new Tdew */ - for (i=0 ; iyday[i]-1; - pva = 610.7 * exp(17.38 * tdew[i] / (239.0 + tdew[i])); - t_tmax = ttmax0[yday] + ABASE * pva; - - /* final daily total transmittance */ - t_final = t_tmax * t_fmax[i]; - - /* estimate fraction of radiation that is diffuse, on an - instantaneous basis, from relationship with daily total - transmittance in Jones (Plants and Microclimate, 1992) - Fig 2.8, p. 25, and Gates (Biophysical Ecology, 1980) - Fig 6.14, p. 122. */ - pdif = -1.25*t_final + 1.25; - if (pdif > 1.0) pdif = 1.0; - if (pdif < 0.0) pdif = 0.0; - - /* estimate fraction of radiation that is direct, on an - instantaneous basis */ - pdir = 1.0 - pdif; - - /* the daily total radiation is estimated as the sum of the - following two components: - 1. The direct radiation arriving during the part of - the day when there is direct beam on the slope. - 2. The diffuse radiation arriving over the entire daylength - (when sun is above ideal horizon). - */ - - /* component 1 */ - srad1 = slope_potrad[yday] * t_final * pdir; - - /* component 2 (diffuse) */ - /* includes the effect of surface albedo in raising the diffuse - radiation for obstructed horizons */ - srad2 = flat_potrad[yday] * t_final * pdif * - (sky_prop + DIF_ALB*(1.0-sky_prop)); - - /* snow pack influence on radiation */ - if (data->s_swe[i] > 0.0) - { - /* snow correction in J/m2/day */ - sc = (1.32 + 0.096 * data->s_swe[i]) * 1e6; - /* convert to W/m2 and check for zero daylength */ - if (daylength[yday] > 0.0) sc /= daylength[yday]; - else sc = 0.0; - /* set a maximum correction of 100 W/m2 */ - if (sc > 100.0) sc = 100.0; - } - else sc = 0.0; - - /* save daily radiation */ - data->s_srad[i] = srad1 + srad2 + sc; - } - - /* Revise estimate of Tdew using new radiation */ - for (i=0 ; is_tmin[i] + 273.15; - pet = calc_pet(data->s_srad[i],data->s_tday[i],pa,data->s_dayl[i]); - - /* calculate ratio (PET/effann_prcp) and correct the dewpoint */ - ratio = pet/parray[i]; - ratio2 = ratio*ratio; - ratio3 = ratio2*ratio; - tdewk = tmink*(-0.127 + 1.121*(1.003 - 1.444*ratio + 12.312*ratio2 - - 32.766*ratio3) + 0.0006*(dtr[i])); - tdew[i] = tdewk - 273.15; - } - } /* end of arid-correction humidity and radiation estimation */ - else - { - printf("Using Tdew=Tmin humidity algorithm\n"); - } - - /* now calculate vapor pressure from tdew */ - for (i=0 ; iouthum) - { - /* output humidity as vapor pressure (Pa) */ - data->s_hum[i] = pva; - } - else - { - /* output humidity as vapor pressure deficit (Pa) */ - /* calculate saturated VP at tday */ - pvs = 610.7 * exp(17.38 * data->s_tday[i]/(239.0+data->s_tday[i])); - vpd = pvs - pva; - if (vpd < 0.0) vpd = 0.0; - data->s_hum[i] = vpd; - } - } /* end for i = ndays loop */ - - /* free local array memory */ - free(dtr); - free(sm_dtr); - free(parray); - free(window); - free(t_fmax); - free(tdew); - free(save_pet); - - return (!ok); -} /* end of calc_srad_humidity_iterative() */ - -/* write_out writes the output file, governed by control flags */ -int write_out(const control_struct *ctrl, const data_struct *data) -{ - int ok=1; - int i,ndays; - char humstr[16]; - int iy; - - iy = ctrl->inyear; - if (ctrl->outhum) strcpy(humstr,"VP"); - else strcpy(humstr, "VPD"); - - ndays = ctrl->ndays; - - - /* write the column header lines */ - if (iy) fprintf(ctrl->outmet.ptr,"%6s","year"); - fprintf(ctrl->outmet.ptr,"%6s%8s%8s%8s%8s%9s%9s%8s\n", - "yday","Tmax","Tmin","Tday","prcp",humstr,"srad","daylen"); - if (iy) fprintf(ctrl->outmet.ptr,"%6s"," "); - fprintf(ctrl->outmet.ptr,"%6s%8s%8s%8s%8s%9s%9s%8s\n", - " ","(deg C)","(deg C)","(deg C)","(cm)","(Pa)","(W m-2)","(s)"); - - /* write each day's output */ - for (i=0 ; ioutmet.ptr,"%6d",data->year[i]); - fprintf(ctrl->outmet.ptr,"%6d%8.2lf%8.2lf%8.2lf%8.2lf%9.2lf%9.2lf%8.0lf\n", - data->yday[i],data->s_tmax[i],data->s_tmin[i],data->s_tday[i], - data->s_prcp[i],data->s_hum[i],data->s_srad[i],data->s_dayl[i]); - } - - fclose(ctrl->outmet.ptr); - return (!ok); -} - -/* data_free frees the memory previously allocated by data_alloc() */ -int data_free(const control_struct *ctrl, data_struct *data) -{ - int ok=1; - if (ctrl->inyear) free(data->year); - free(data->yday); - free(data->tmax); - free(data->tmin); - free(data->prcp); - if (ctrl->indewpt) free(data->tdew); - free(data->s_tmax); - free(data->s_tmin); - free(data->s_tday); - free(data->s_prcp); - free(data->s_hum); - free(data->s_srad); - free(data->s_dayl); - - return (!ok); -} - -/* calc_pet() calculates the potential evapotranspiration for aridity -corrections in calc_vpd(), according to Kimball et al., 1997 */ -double calc_pet(double rad, double ta, double pa, double dayl) -{ - /* input parameters and units : - double rad (W/m2) daylight average incident shortwave radiation - double ta (deg C) daylight average air temperature - double pa (Pa) air pressure - double dayl (s) daylength - */ - - double rnet; /* (W m-2) absorbed shortwave radiation avail. for ET */ - double lhvap; /* (J kg-1) latent heat of vaporization of water */ - double gamma; /* (Pa K-1) psychrometer parameter */ - double dt = 0.2; /* offset for saturation vapor pressure calculation */ - double t1, t2; /* (deg C) air temperatures */ - double pvs1, pvs2; /* (Pa) saturated vapor pressures */ - double pet; /* (kg m-2 day-1) potential evapotranspiration */ - double s; /* (Pa K-1) slope of saturated vapor pressure curve */ - - /* calculate absorbed radiation, assuming albedo = 0.2 and ground - heat flux = 10% of absorbed radiation during daylight */ - rnet = rad * 0.72; - - /* calculate latent heat of vaporization as a function of ta */ - lhvap = 2.5023e6 - 2430.54 * ta; - - /* calculate the psychrometer parameter: gamma = (cp pa)/(lhvap epsilon) - where: - cp (J/kg K) specific heat of air - epsilon (unitless) ratio of molecular weights of water and air - */ - gamma = CP * pa / (lhvap * EPS); - - /* estimate the slope of the saturation vapor pressure curve at ta */ - /* temperature offsets for slope estimate */ - t1 = ta+dt; - t2 = ta-dt; - - /* calculate saturation vapor pressures at t1 and t2, using formula from - Abbott, P.F., and R.C. Tabony, 1985. The estimation of humidity parameters. - Meteorol. Mag., 114:49-56. - */ - pvs1 = 610.7 * exp(17.38 * t1 / (239.0 + t1)); - pvs2 = 610.7 * exp(17.38 * t2 / (239.0 + t2)); - - /* calculate slope of pvs vs. T curve near ta */ - s = (pvs1-pvs2) / (t1-t2); - - /* calculate PET using Priestly-Taylor approximation, with coefficient - set at 1.26. Units of result are kg/m^2/day, equivalent to mm water/day */ - pet = (1.26 * (s/(s+gamma)) * rnet * dayl)/lhvap; - - /* return a value in centimeters/day, because this value is used in a ratio - to annual total precip, and precip units are centimeters */ - return (pet/10.0); -} - -/* atm_pres() calculates the atmospheric pressure as a function of elevation */ -double atm_pres(double elev) -{ - /* daily atmospheric pressure (Pa) as a function of elevation (m) */ - /* From the discussion on atmospheric statics in: - Iribane, J.V., and W.L. Godson, 1981. Atmospheric Thermodynamics, 2nd - Edition. D. Reidel Publishing Company, Dordrecht, The Netherlands. - (p. 168) - */ - - int ok=1; - double t1,t2; - double pa; - - t1 = 1.0 - (LR_STD * elev)/T_STD; - t2 = G_STD / (LR_STD * (R / MA)); - pa = P_STD * pow(t1,t2); - - return(pa); -} - -/* pulled_boxcar() calculates a moving average of antecedent values in an -array, using either a ramped (w_flag=1) or a flat (w_flag=0) weighting */ -int pulled_boxcar(double *input,double *output,int n,int w,int w_flag) -{ - int ok=1; - int i,j; - double *wt; - double total,sum_wt; - - if (w > n) { - printf("Boxcar longer than array...\n"); - printf("Resize boxcar and try again\n"); - ok=0; - } - - if (ok && !(wt = (double*) malloc(w * sizeof(double)))) - { - printf("Allocation error in boxcar()\n"); - ok=0; - } - - if (ok) - { - /* when w_flag != 0, use linear ramp to weight tails, - otherwise use constant weight */ - sum_wt = 0.0; - if (w_flag) - { - for (i=0 ; i -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -using namespace Rcpp; -using namespace std; -// [[Rcpp::plugins(cpp11)]] - -// [[Rcpp::export]] -IntegerVector getWritePositions(double a){ - //getWritePositions returns abstract rownumbers to rownumbers and other indexek - // getWritePositions(173.62) = c(173,7,3) // it supports up to 10 subvalues - IntegerVector outVec(3); // outVec is vector of rowIndex, colNumber, choosen Index - a = a * 100; - a = round(a); //without this line 155.92 ~= 155.9199999999, (int) 155.9199999999*100 = 15591 - outVec[0] = (int)a / 100; - outVec[1] = ((int)a /10) % 10 + 1; - outVec[2] = (int)a % 10; - - - return outVec; -} - -IntegerMatrix getPositions(NumericVector v){ - - int numVari = v.size(); - IntegerMatrix indexek(numVari,3); - IntegerVector positions(3); - for(int i = 0; i < numVari; ++i){ - positions = getWritePositions(v[i]); - indexek(i,_) = positions; - } - return indexek; -} - - -void goNextLine(std::ifstream& fin){ - char c='a'; - while((c!='\n')&&(fin.get(c))){} -} - -#define NEXT goNextLine(fin) -int fileChanger(std::string inFile, std::string outFile, IntegerVector linum, NumericVector num, IntegerVector colnum, IntegerVector colindex){ - std::ifstream fin(inFile); - if (!fin.is_open()) { - stop("Cannot open " + inFile + " for read"); - } - std::ofstream fot(outFile); - - if (!fot.is_open()) { - stop("Cannot open " + outFile + " for read"); - } - string tempString; - - int counter = 1; - int counterV = 0; - while (!fin.eof()) { - - - if(counter == linum[counterV]){ - if(colnum[counterV]==1){ - fot << num[counterV] << "\n"; - NEXT; - } else { - double * elements; - elements = new double [colnum[counterV]]; - if(linum[counterV]!=linum[counterV+1]){ - for(int i=0;i> elements[i]; - - if(i==colindex[counterV]){ - elements[i]=num[counterV]; - } - // std::cout << colnum[counterV] << " " << colindex[counterV] << " " << elements[i] << " "<< i <<"\n"; - // std::cout << colindex[counterV] << getWritePositions(155.92) <<"\n"; - // std::cout << "======================== \n"; - fot << elements[i] << '\t'; - } - } else { - int k=0; - for(int i=0;i> elements[i]; - if(i==colindex[counterV + k]){ - elements[i]=num[counterV + k]; - if(linum[counterV +k]==linum[counterV + k + 1]){ - ++k; - } - } - fot << elements[i] << '\t'; - - } - counterV = counterV + k; - } - - fot << "\n"; - delete [] elements; - NEXT; - } - ++counterV; - } else { - getline(fin,tempString); - fot << tempString << "\n"; - } - ++counter; - if(counter > 1000){ - stop("You modified a line which has not as many columns as you specified."); - } - } - fin.close(); - fot.close(); - return 0; -} - -//' changeMusoC -//' -//' This function is fastly randomize values based on min and max values -//' @importFrom Rcpp evalCpp -//' @useDynLib RBBGCMuso -//' @param inFile is the big matrix -//' @param outFile is the small matrix -//' @export -// [[Rcpp::export]] -void changeMusoC(std::string inFile, std::string outFile, NumericMatrix inMat){ - int numChanges = inMat.nrow(); - IntegerMatrix indexes(numChanges,3); - indexes = getPositions(inMat(_,0)); - fileChanger(inFile,outFile,indexes(_,0),inMat(_,1),indexes(_,1),indexes(_,2)); -} diff --git a/RBBGCMuso/src/musoRandomizer.cpp b/RBBGCMuso/src/musoRandomizer.cpp deleted file mode 100644 index f9e143a..0000000 --- a/RBBGCMuso/src/musoRandomizer.cpp +++ /dev/null @@ -1,245 +0,0 @@ -#include -#include -#include -#include -#include -#include - - -using namespace Rcpp; -using namespace std; -// [[Rcpp::plugins(cpp11)]] - - - -NumericMatrix randTypeZero(NumericMatrix m){ - /* - A typical matrix like m: - - | INDEX | DEPENDENCE | MIN | MAX | - |-------+------------+------+------| - | 21 | 0 | 0 | 364 | - | 57 | 0 | bkla | asdf | - - This randomization type is the easiest, - the function produces a matrix which first - column contains the indexes, and the second contains - random numbers, which was drawn from uniform distribution - with the corresponding min and max parameters, specified - by the matrix m. - - */ - int n=m.nrow()-1; - NumericMatrix M(n+1,2); - M(_,0)=m(_,0); - for(int i=0;i<=n;++i){ - double min=m(i,2); - double max=m(i,3); - M(i,1)=runif(1,min,max)[0]; - } - return M; -} - -// [[Rcpp::export]] -NumericMatrix randTypeOne(NumericMatrix m){ - /* - The stucture of matrix m here is the same - as in randTypeZero. - - - - */ - - - NumericVector dependence=m(_,2); - int n=m.nrow()-1; - NumericMatrix M(n+1,2); - M(_,0)=m(_,0); - M(0,1)=runif(1,m(0,2),m(0,3))[0]; - for(int i=1;i<=n;++i){ - int dep=m(i,1)-1; - double min=max(M(dep,1),m(i,2)); - double max=m(i,3); - M(i,1)=runif(1,min,max)[0]; - } - return M; -} - - - -IntegerVector orderDec(NumericVector v){ - //This function is order a vector decreasingly - Function f("order"); - return f(v,_["decreasing"]=1); -} - -// NumericMatrix randTypeTwo(NumericMatrix m){ -// int n=m.nrow()-1; -// int N=n-1; -// NumericMatrix mv=m(Range(0,(n-1)),_); -// NumericVector dependence=m(_,2); -// NumericMatrix M(n+1,2); -// M(_,0)=m(_,0); -// IntegerVector indexes=orderDec(mv(_,2)); -// NumericVector sorban=mv(_,2); -// sorban.sort(true); -// NumericVector sor=cumsum(sorban); -// sor.sort(true); -// for(int i=0;i<=N;++i){ -// if(i!=N){ -// mv((indexes[i]-1),3)-= sor[i+1]; -// } -// } - -// double rollingNumber=0; - -// for(int i=0;i<=N;++i){ -// double minimum=mv((indexes[i]-1),2); -// double maximum=mv((indexes[i]-1),3)-rollingNumber; -// M(i,1)=runif(1,minimum,maximum)[0]; -// rollingNumber+=M(i,1); -// // cout << "minimum:\t" << minimum << endl; -// // cout << "maximum:\t" << maximum << endl; -// // cout << "indexes:\t" << indexes[i] << endl; -// // cout << "rollingNumber:\t" << rollingNumber << endl; -// // cout << "choosen:\t" << M(i,1) < Date: Mon, 28 Sep 2020 07:02:39 +0200 Subject: [PATCH 7/7] filepath bug in calibration --- RBBGCMuso/R/calibrateMuso.R | 31 +++++++++---------------------- 1 file changed, 9 insertions(+), 22 deletions(-) diff --git a/RBBGCMuso/R/calibrateMuso.R b/RBBGCMuso/R/calibrateMuso.R index 89c643e..efea427 100644 --- a/RBBGCMuso/R/calibrateMuso.R +++ b/RBBGCMuso/R/calibrateMuso.R @@ -14,7 +14,7 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, naVal = NULL, postProcString = NULL, thread_prefix="thread", numCores = (parallel::detectCores()-1), pb = txtProgressBar(min=0, max=iterations, style=3), maxLikelihoodEpc=TRUE, - pbUpdate = setTxtProgressBar, method="GLUE",lg = FALSE, w=NULL, ...){ + pbUpdate = setTxtProgressBar, outputLoc="./", method="GLUE",lg = FALSE, w=NULL, ...){ future::plan(future::multisession) file.remove(list.files(path = settings$inputLoc, pattern="progress.txt", recursive = TRUE, full.names=TRUE)) @@ -42,14 +42,16 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, fut <- lapply(1:numCores, function(i) { # browser() future({ - tryCatch(musoSingleThread(measuredData, parameters, startDate, + tryCatch( + musoSingleThread(measuredData, parameters, startDate, endDate, formatString, dataVar, outLoc, preTag, settings, outVars, iterations = threadCount[i], skipSpinup, plotName, modifyOriginal, likelihood, uncertainity, - naVal, postProcString, i), error = function(e){ + naVal, postProcString, i) + , error = function(e){ writeLines(as.character(iterations),"progress.txt") }) @@ -109,7 +111,6 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, # | | / _ \| '_ ` _ \| '_ \| | '_ \ / _ \ # | |__| (_) | | | | | | |_) | | | | | __/ # \____\___/|_| |_| |_|_.__/|_|_| |_|\___| - resultFiles <- list.files(pattern="preservedCalib.*csv$",recursive=TRUE) res0 <- read.csv(grep("thread_1/",resultFiles, value=TRUE),stringsAsFactors=FALSE) resultFilesSans0 <- grep("thread_1/", resultFiles, value=TRUE, invert=TRUE) @@ -129,8 +130,8 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, epcVals <- results[which.max(liks),1:length(epcIndexes)] epcPlace <- file.path(dirname(settings$inputFiles),settings$epc)[2] changemulline(filePaths= epcPlace, epcIndexes, - epcVals, src = settings$epcInput[2], - outFiles = "maxLikelihood_epc.epc") + epcVals, src =epcPlace,# settings$epcInput[2], + outFiles = file.path(outputLoc, "maxLikelihood_epc.epc")) names(epcVals) <- epcIndexes xdate <- as.Date(measuredData$date) meanM <- measuredData[,sprintf("mean.%s", names(likelihood))] @@ -156,19 +157,6 @@ calibrateMuso <- function(measuredData, parameters = NULL, startDate = NULL, }, stop(sprintf("method: %s not found, please choose from {GLUE, agromo}. See more about this in the documentation of the function!", method)) ) - # Here starts maxLikelihoodAgroMo: parameters - - - # Here ends maxLikelihoodAgromo - - # return(epcVals) - # ____ _ _ _ _____ - # / ___| | | | | | ____| - # | | _| | | | | | _| - # | |_| | |__| |_| | |___ - # \____|_____\___/|_____| - - # musoGlue("preservedCalib.csv",w=w, lg = lg) } copyToThreadDirs <- function(prefix="thread", numcores=parallel::detectCores()-1, runDir="."){ @@ -195,6 +183,7 @@ musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, if(length(iniFiles)==1){ iniFiles <- rep(iniFiles, 2) } + iniFiles <- iniFiles[1:2] settings <- setupMuso(iniInput = iniFiles) # Exanding likelihood likelihoodFull <- as.list(rep(NA,length(dataVar))) @@ -305,7 +294,6 @@ musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, for(i in 2:(iterations+1)){ - tmp <- tryCatch(calibMuso(settings = settings, parameters = randValues[(i-1),], silent= TRUE, @@ -321,7 +309,6 @@ musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, musoCodeToIndex = musoCodeToIndex, uncert = uncert) } - partialResult[1:numParameters] <- randValues[(i-1),] write.table(x=partialResult, file="preservedCalib.csv", append=TRUE, row.names=FALSE, sep=",", col.names=FALSE) @@ -356,7 +343,7 @@ calcLikelihoodsAndRMSE <- function(dataVar, mod, mes, likelihoods, alignIndexes, # NOT COMPATIBLE WITH OLD MEASUREMENT DATA, mes have to be a matrix likelihoodRMSE <- sapply(names(dataVar),function(key){ - # browser() + # browser() modelled <- mod[alignIndexes$mod,musoCodeToIndex[key]] selected <- grep(sprintf("%s$", key), colnames(mes)) # browser()