hollorol/RBBGCMuso
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Modify many

Browse Source
This commit is contained in:
hollorol 2019-02-19 09:48:29 +01:00
parent 4ab2b60f49
commit 0d9d551fd9
10 changed files with 74 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
RBBGCMuso/DESCRIPTION
View File

@ -1,6 +1,6 @@
Package: RBBGCMuso Package: RBBGCMuso
Title: An R package for BiomeBGC-MuSo ecosystem modelling Title: An R package for BiomeBGC-MuSo ecosystem modelling
Version: 0.7.0.0 Version: 0.7.0.1
Authors@R: person("Roland", "Hollo's", , "hollorol@gmail.com", role = c("aut", "cre")) Authors@R: person("Roland", "Hollo's", , "hollorol@gmail.com", role = c("aut", "cre"))
Description: What the package does (one paragraph). Description: What the package does (one paragraph).
Depends: R (>= 3.3.2) Depends: R (>= 3.3.2)

2
RBBGCMuso/NAMESPACE
View File

@ -23,7 +23,7 @@ export(paramSweep)
export(plotMuso) export(plotMuso)
export(plotMusoWithData) export(plotMusoWithData)
export(randEpc) export(randEpc)
export(readMeasuredMuso) export(readObservedData)
export(runMuso) export(runMuso)
export(rungetMuso) export(rungetMuso)
export(saveAllMusoPlots) export(saveAllMusoPlots)

5
RBBGCMuso/R/assistantFunctions.R
View File

@ -146,7 +146,7 @@ readValuesFromFile <- function(epc, linums){
#' @importFrom data.table fread data.table #' @importFrom data.table fread data.table
#' @export #' @export
readMeasuredMuso <- function(inFile, readObservedData <- function(inFile,
naString = NULL, sep = ",", naString = NULL, sep = ",",
leapYearHandling = TRUE, leapYearHandling = TRUE,
convert.var = NULL, convert.var = NULL,
@ -181,8 +181,7 @@ readMeasuredMuso <- function(inFile,
} }
head(baseData) head(baseData)
if(!is.null(selVar)){ if(!is.null(selVar)){
baseData <- cbind.data.frame(baseData,convert.fun(baseData[,selVar])) baseData[,selVar] <-convert.fun(baseData[,selVar])
colnames(baseData)[ncol(baseData)]<- paste0("M",selVar)
} }
return(data.table(baseData)) return(data.table(baseData))

22
RBBGCMuso/R/calibration.R
View File

@ -28,7 +28,7 @@
#' @export #' @export
optiMuso <- function(measuredData, parameters = NULL, startDate, optiMuso <- function(measuredData, parameters = NULL, startDate,
endDate, formatString = "%Y-%m-%d", endDate, formatString = "%Y-%m-%d",
leapYear = TRUE, leapYearHandling = TRUE,
dataVar, outLoc = "./calib", dataVar, outLoc = "./calib",
preTag = "cal-", preTag = "cal-",
settings = NULL, settings = NULL,
@ -40,7 +40,8 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
likelihood = function(x, y){ likelihood = function(x, y){
exp(-sqrt(mean((x-y)^2))) exp(-sqrt(mean((x-y)^2)))
}, },
calPar = 3009) continious,
modelVar = 3009)
{ {
dataCol <- grep(dataVar, colnames(measuredData)) dataCol <- grep(dataVar, colnames(measuredData))
@ -76,6 +77,7 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
##reading the original epc file at the specified ##reading the original epc file at the specified
## row numbers ## row numbers
print("optiMuso is randomizing the epc parameters now...",quote = FALSE)
if(iterations < 3000){ if(iterations < 3000){
randVals <- musoRand(parameters = parameters,constrains = constrains, iterations = 3000) randVals <- musoRand(parameters = parameters,constrains = constrains, iterations = 3000)
randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),] randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),]
@ -94,18 +96,22 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
## csv files for each run ## csv files for each run
progBar <- txtProgressBar(1,iterations,style=3) progBar <- txtProgressBar(1,iterations,style=3)
colNumb <- which(settings$dailyVarCodes == calPar) colNumb <- which(settings$dailyVarCodes == modelVar)
settings$iniInput[2] %>% settings$iniInput[2] %>%
(function(x) paste0(dirname(x),"/",tools::file_path_sans_ext(basename(x)),"-tmp.",tools::file_ext(x))) %>% (function(x) paste0(dirname(x),"/",tools::file_path_sans_ext(basename(x)),"-tmp.",tools::file_ext(x))) %>%
unlink unlink
randValues <- randVals[[2]] randValues <- randVals[[2]]
settings$calibrationPar <- randVals[[1]] settings$calibrationPar <- randVals[[1]]
list2env(alignData(measuredData,dataCol = dataCol,modellSettings = settings,startDate = startDate,endDate = endDate,leapYear = FALSE),envir=environment()) list2env(alignData(measuredData,dataCol = dataCol,modellSettings = settings,startDate = startDate,endDate = endDate,leapYear = leapYearHandling, continious = continious),envir=environment())
modellOut <- numeric(iterations + 1) # single variable solution modellOut <- numeric(iterations + 1) # single variable solution
origModellOut <- calibMuso(settings=settings,silent=TRUE) rmse <- numeric(iterations + 1)
origModellOut <- calibMuso(settings=settings,silent=TRUE, skipSpinup = skipSpinup)
write.csv(x=origModellOut, file=paste0(pretag,1,".csv")) write.csv(x=origModellOut, file=paste0(pretag,1,".csv"))
modellOut[1] <- likelihood(measuredData,origModellOut[modIndex,colNumb]) modellOut[1] <- likelihood(measuredData,origModellOut[modIndex,colNumb])
print("Running the model with the random epc values...", quote = FALSE)
for(i in 2:(iterations+1)){ for(i in 2:(iterations+1)){
tmp <- tryCatch(calibMuso(settings = settings, tmp <- tryCatch(calibMuso(settings = settings,
parameters = randValues[(i-1),], parameters = randValues[(i-1),],
@ -113,10 +119,10 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
skipSpinup = skipSpinup)[modIndex,colNumb], error = function (e) NA) skipSpinup = skipSpinup)[modIndex,colNumb], error = function (e) NA)
modellOut[i]<- likelihood(measuredData,tmp) modellOut[i]<- likelihood(measuredData,tmp)
rmse[i] <- sqrt(mean((measuredData-tmp)^2))
write.csv(x=tmp, file=paste0(pretag,(i+1),".csv")) write.csv(x=tmp, file=paste0(pretag,(i+1),".csv"))
setTxtProgressBar(progBar,i) setTxtProgressBar(progBar,i)
} }
modellOut
paramLines <- parameters[,2] paramLines <- parameters[,2]
paramLines <- order(paramLines) paramLines <- order(paramLines)
randInd <- randVals[[1]][(randVals[[1]] %in% parameters[,2])] randInd <- randVals[[1]][(randVals[[1]] %in% parameters[,2])]
@ -128,7 +134,7 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
randValues[,randVals[[1]] %in% parameters[,2]][,randInd]) randValues[,randVals[[1]] %in% parameters[,2]][,randInd])
preservedCalib <- cbind(epcStrip, preservedCalib <- cbind(epcStrip,rmsr,
modellOut) modellOut)
colnames(preservedCalib) <- c(parameterNames[paramLines], "likelihood") colnames(preservedCalib) <- c(parameterNames[paramLines], "likelihood")
p<-list() p<-list()
@ -137,7 +143,7 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
p[[i]] <- ggplot(as.data.frame(preservedCalib),aes_string(colnames(preservedCalib)[i],"likelihood"))+geom_point(size=0.9) p[[i]] <- ggplot(as.data.frame(preservedCalib),aes_string(colnames(preservedCalib)[i],"likelihood"))+geom_point(size=0.9)
} }
ggsave(plotName,grid.arrange(grobs = p, ncol = floor(sqrt(ncol(preservedCalib)-1))),dpi = 600) ggsave(plotName,grid.arrange(grobs = p, ncol = floor(sqrt(ncol(preservedCalib)-1))),dpi = 3000)
write.csv(preservedCalib,"preservedCalib.csv") write.csv(preservedCalib,"preservedCalib.csv")
return(preservedCalib[preservedCalib[,"likelihood"]==max(preservedCalib[,"likelihood"]),]) return(preservedCalib[preservedCalib[,"likelihood"]==max(preservedCalib[,"likelihood"]),])
} }

23
RBBGCMuso/R/musoTime.R
View File

@ -11,7 +11,6 @@
#' @importFrom lubridate leap_year #' @importFrom lubridate leap_year
#' @export #' @export
musoDate <- function(startYear, endYears = NULL, numYears, combined = TRUE, leapYearHandling = FALSE, prettyOut = FALSE){ musoDate <- function(startYear, endYears = NULL, numYears, combined = TRUE, leapYearHandling = FALSE, prettyOut = FALSE){
if(is.null(endYears) & is.null(numYears)){ if(is.null(endYears) & is.null(numYears)){
@ -60,7 +59,8 @@ musoDate <- function(startYear, endYears = NULL, numYears, combined = TRUE, leap
#' This function align the data to the model and the model to the data #' This function align the data to the model and the model to the data
#' @importFrom lubridate leap_year #' @importFrom lubridate leap_year
#' @keywords internal #' @keywords internal
alignData <- function(mdata, dataCol, modellSettings = NULL, startDate, endDate, formatString = "%Y-%m-%d", leapYear = TRUE){ alignData <- function(mdata, dataCol, modellSettings = NULL, startDate=NULL, endDate=NULL, formatString = "%Y-%m-%d", leapYear = TRUE, continious = TRUE){
startDate <- as.Date(startDate, format = formatString) startDate <- as.Date(startDate, format = formatString)
endDate <- as.Date(endDate, format = formatString) endDate <- as.Date(endDate, format = formatString)
@ -70,15 +70,24 @@ alignData <- function(mdata, dataCol, modellSettings = NULL, startDate, endDate
modellSettings <- setupMuso() modellSettings <- setupMuso()
} }
if(continious){
dates <- seq(startDate, to = endDate, by= "day") dates <- seq(startDate, to = endDate, by= "day")
if(!leapYear){ } else{
dates <- dates[which(format(dates,"%m%d") != "0229")] dates <- do.call(c,lapply(seq(nrow(mdata)), function(i){ as.Date(paste0(mdata[i,1],sprintf("%02d",mdata[i,2]),mdata[i,3]),format = "%Y%m%d")}))
} }
if(!leapYear){
casualDays <- which(format(dates,"%m%d") != "0229")
#mdata <- mdata[casualDays,]
dates <- dates[casualDays]
}
mdata <- mdata[dates >= as.Date(paste0(modellSettings$startYear,"01","01"),format = "%Y%m%d"),] mdata <- mdata[dates >= as.Date(paste0(modellSettings$startYear,"01","01"),format = "%Y%m%d"),]
dates <- dates[dates >= as.Date(paste0(modellSettings$startYear,"01","01"),format = "%Y%m%d")] dates <- dates[dates >= as.Date(paste0(modellSettings$startYear,"01","01"),format = "%Y%m%d")]
goodInd <- which(!(leap_year(dates)& ## goodInd <- which(!(leap_year(dates)&
(format(dates,"%m") == "12")& ## (format(dates,"%m") == "12")&
(format(dates,"%d") == "31"))) ## (format(dates,"%d") == "31")))
if(leapYear){ if(leapYear){
goodInd <- which(!(leap_year(dates)& goodInd <- which(!(leap_year(dates)&
(format(dates,"%m") == "12")& (format(dates,"%m") == "12")&

6
RBBGCMuso/R/plotMuso.R
View File

@ -251,7 +251,7 @@ plotMuso <- function(settings = NULL, variable = 1,
#' @export #' @export
plotMusoWithData <- function(mdata, plotName=NULL, plotMusoWithData <- function(mdata, plotName=NULL,
startDate, endDate, startDate, endDate,
colour=c("black","blue"),dataVar, modelVar, settings = setupMuso(), silent = TRUE){ colour=c("black","blue"),dataVar, modelVar, settings = setupMuso(), silent = TRUE, continious = TRUE){
dataCol<- grep(paste0("^",dataVar,"$"), colnames(mdata)) dataCol<- grep(paste0("^",dataVar,"$"), colnames(mdata))
selVar <- grep(modelVar,(settings$dailyVarCodes))+4 selVar <- grep(modelVar,(settings$dailyVarCodes))+4
@ -259,14 +259,14 @@ plotMusoWithData <- function(mdata, plotName=NULL,
list2env(alignData(mdata, dataCol = dataCol, list2env(alignData(mdata, dataCol = dataCol,
modellSettings = settings, modellSettings = settings,
startDate = startDate, startDate = startDate,
endDate = endDate, leapYear = FALSE),envir=environment()) endDate = endDate, leapYear = FALSE, continious = continious),envir=environment())
## measuredData is created ## measuredData is created
baseData <- calibMuso(settings = settings, silent = silent, prettyOut = TRUE)[modIndex,] baseData <- calibMuso(settings = settings, silent = silent, prettyOut = TRUE)[modIndex,]
baseData[,1] <- as.Date(baseData[,1],format = "%d.%m.%Y") baseData[,1] <- as.Date(baseData[,1],format = "%d.%m.%Y")
selVarName <- colnames(baseData)[selVar] selVarName <- colnames(baseData)[selVar]
if(colnames(baseData) != unique(colnames(baseData))){ if(!all.equal(colnames(baseData),unique(colnames(baseData)))){
notUnique <- setdiff((unlist(settings$dailyVarCodes)),unique(unlist(settings$dailyVarCodes))) notUnique <- setdiff((unlist(settings$dailyVarCodes)),unique(unlist(settings$dailyVarCodes)))
stop(paste0("Error: daily output variable list in the ini file must contain unique numbers. Check your ini files! Not unique codes: ",notUnique)) stop(paste0("Error: daily output variable list in the ini file must contain unique numbers. Check your ini files! Not unique codes: ",notUnique))
} }

5
RBBGCMuso/man/alignData.Rd
View File

@ -4,8 +4,9 @@
\alias{alignData} \alias{alignData}
\title{alignData} \title{alignData}
\usage{ \usage{
alignData(mdata, dataCol, modellSettings = NULL, startDate, endDate, alignData(mdata, dataCol, modellSettings = NULL, startDate = NULL,
formatString = "\%Y-\%m-\%d", leapYear = TRUE) endDate = NULL, formatString = "\%Y-\%m-\%d", leapYear = TRUE,
continious = TRUE)
} }
\description{ \description{
This function align the data to the model and the model to the data This function align the data to the model and the model to the data

34
RBBGCMuso/man/optiMuso.Rd
View File

@ -4,22 +4,16 @@
\alias{optiMuso} \alias{optiMuso}
\title{optiMuso} \title{optiMuso}
\usage{ \usage{
optiMuso(measuredDataFile, parameters = NULL, sep = ",", startDate, optiMuso(measuredData, parameters = NULL, startDate, endDate,
endDate, formatString = "\%Y-\%m-\%d", formatString = "\%Y-\%m-\%d", leapYear = TRUE, dataVar,
naString = getOption("datatable.na.strings", "NA"), leapYear = TRUE, outLoc = "./calib", preTag = "cal-", settings = NULL,
filterCol = NULL, filterVal = 1, selVar, outLoc = "./calib", outVars = NULL, iterations = 30, skipSpinup = TRUE,
preTag = "cal-", settings = NULL, outVars = NULL, constrains = NULL, plotName = "calib.jpg", likelihood = function(x,
iterations = 30, skipSpinup = TRUE, constrains = NULL, y) { exp(-sqrt(mean((x - y)^2))) }, calPar = 3009)
plotName = "calib.jpg", likelihood = function(x, y) {
exp(-sqrt(mean((x - y)^2))) }, calPar = 3009)
} }
\arguments{ \arguments{
\item{measuredDataFile}{a}
\item{parameters}{b} \item{parameters}{b}
\item{sep}{c}
\item{startDate}{d} \item{startDate}{d}
\item{endDate}{e} \item{endDate}{e}
@ -28,12 +22,6 @@ optiMuso(measuredDataFile, parameters = NULL, sep = ",", startDate,
\item{leapYear}{b} \item{leapYear}{b}
\item{filterCol}{a}
\item{filterVal}{b}
\item{selVar}{c}
\item{outLoc}{c} \item{outLoc}{c}
\item{settings}{e} \item{settings}{e}
@ -50,6 +38,16 @@ optiMuso(measuredDataFile, parameters = NULL, sep = ",", startDate,
\item{calPar}{a} \item{calPar}{a}
\item{measuredDataFile}{a}
\item{sep}{c}
\item{filterCol}{a}
\item{filterVal}{b}
\item{selVar}{c}
\item{pretag}{a} \item{pretag}{a}
} }
\description{ \description{

15
RBBGCMuso/man/readMeasuredMuso.Rd
View File

@ -1,15 +0,0 @@
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/assistantFunctions.R
\name{readMeasuredMuso}
\alias{readMeasuredMuso}
\title{readMeasuredMuso}
\usage{
readMeasuredMuso(inFile, naString = getOption("datatable.na.strings",
"NA"), sep = ",", leapYearHandling = TRUE, convert.var = NULL,
convert.scalar = 1, convert.fun = (function(x) { x *
convert.scalar }), convert.file = NULL, filterCol = NULL,
filterVal = 1, selVar = NULL)
}
\description{
MuSo data reader
}

15
RBBGCMuso/man/readObservedData.Rd Normal file
View File

@ -0,0 +1,15 @@
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/assistantFunctions.R
\name{readObservedData}
\alias{readObservedData}
\title{readMeasuredMuso}
\usage{
readObservedData(inFile, naString = NULL, sep = ",",
leapYearHandling = TRUE, convert.var = NULL, convert.scalar = 1,
convert.fun = (function(x) { x * convert.scalar }),
convert.file = NULL, filterCol = NULL, filterVal = 1,
selVar = NULL)
}
\description{
MuSo data reader
}