-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy path.Rhistory
176 lines (176 loc) · 8.19 KB
/
.Rhistory
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
getwd()
setwd("~/Rpackages/SOS/R/Model_March2017")
setwd('../..')
setwd("~/Documents/SOS")
LakeName = 'Toolik'
##### LOAD PACKAGES ########################
library(lubridate)
library(LakeMetabolizer)
library(dplyr)
library(FME)
library(zoo)
library(plyr)
##### LOAD FUNCTIONS #######################
source("./R/Model_March2017/SOS_Sedimentation.R")
source("./R/Model_March2017/SOS_SWGW.R")
source("./R/Model_March2017/SOS_GPP.R")
source("./R/Model_March2017/SOS_Resp.R")
source("./R/Model_March2017/modelDOC_7.R")
source("./R/Model_March2017/SOS_fixToolik.R")
##### INPUT FILE NAMES ################
TimeSeriesFile <- paste('./',LakeName,'Lake/',LakeName,'TS.csv',sep='')
RainFile <- paste('./',LakeName,'Lake/',LakeName,'Rain.csv',sep='')
ParameterFile <- paste('./',LakeName,'Lake/','ConfigurationInputs',LakeName,'.txt',sep='')
FreeParFile <- paste('./R/FMEresults/',LakeName,'_fitpars.csv',sep='')
ValidationFileDOC <- paste('./',LakeName,'Lake/',LakeName,'ValidationDOC.csv',sep='')
ValidationFileDO <- paste('./',LakeName,'Lake/',LakeName,'ValidationDO.csv',sep='')
timestampFormat = '%Y-%m-%d'
InflowQFile = paste0('./R/Loadest/',LakeName,'/','observed.csv')
InflowDOCFile = paste0('./R/Loadest/',LakeName,'/loadestComp.csv')
##### READ MAIN INPUT FILE #################
RawData <- read.csv(TimeSeriesFile,header=T,stringsAsFactors = F) #Read main data file with GLM outputs (physical input) and NPP input
RawData$datetime <- as.POSIXct(strptime(RawData$datetime,timestampFormat),tz="GMT") #Convert time to POSIX
cc = which(complete.cases(RawData))
RawData = RawData[cc[1]:tail(cc,1),]
# Read LOADEST files
inflowQ = read.csv(InflowQFile,stringsAsFactors = F)
inflowQ$Dates = as.POSIXct(strptime(inflowQ$Dates,timestampFormat),tz="GMT") #Convert time to POSIX
inflowDOC = read.csv(InflowDOCFile,stringsAsFactors = F)
inflowDOC$date = as.POSIXct(strptime(inflowDOC$date,timestampFormat),tz="GMT") #Convert time to POSIX
RawData = RawData %>% left_join(inflowQ,by=c('datetime'='Dates')) %>%
left_join(inflowDOC,by=c('datetime'='date')) %>%
dplyr::select(datetime,Volume,FlowIn = 'Flow',HypoTemp:Chla,SW_DOC='fit',Secchi) %>%
mutate(FlowOut = FlowIn)
# Fill time-series gaps (linear interpolation)
ts_new <- data.frame(datetime = seq(RawData$datetime[1],RawData$datetime[nrow(RawData)],by="day")) #Interpolate gapless time-series
InputData <- merge(RawData,ts_new,all=T)
for (col in 2:ncol(InputData)){
InputData[,col] <- na.approx(InputData[,col],na.rm = T)}
InputData$Chla[InputData$Chla == 0] = 0.0001
##### READ RAIN FILE #######################
RainData <- read.csv(RainFile,header=T,stringsAsFactors = F) #Read daily rain file (units=mm) Read separately and added back to main file to avoid issues of linear interpolation with rain data in length units
RainData$datetime <- as.POSIXct(strptime(RainData$datetime,timestampFormat,tz='GMT'))
InputData$Rain <- RainData$Rain[RainData$datetime %in% InputData$datetime] #Plug daily rain data into InputData file to integrate with original code.
#### For TOOLIK ONLY #### (dealing with ice season)
if (LakeName=='Toolik') {
InputData = fixToolik(InputData,LakeName)
}
#DOC Validation Output Setup
ValidationDataDOC <- read.csv(ValidationFileDOC,header=T,stringsAsFactors = F)
outlier.limit = (mean(ValidationDataDOC$DOC,na.rm=T) + 3*(sd(ValidationDataDOC$DOC,na.rm=T))) # Calculate mean + 3 SD of DOC column
ValidationDataDOC = ValidationDataDOC %>% mutate(datetime = (as.POSIXct(strptime(datetime,timestampFormat),tz="GMT"))) %>%
dplyr::filter(complete.cases(.)) %>%
dplyr::filter(DOC <= outlier.limit) %>%
group_by(datetime) %>%
summarise_all(mean) %>%
mutate(datetime = as.Date(datetime))
#DO Validation Output Setup
ValidationDataDO <- read.csv(ValidationFileDO,header=T)
k <- 0.7 #m/d
PhoticDepth <- data.frame(datetime = InputData$datetime,PhoticDepth = log(100)/(1.7/InputData$Secchi))
ValidationDataDO = ValidationDataDO %>% mutate(datetime = as.POSIXct(strptime(datetime,timestampFormat),tz="GMT")) %>%
dplyr::filter(complete.cases(.)) %>%
inner_join(PhoticDepth,by='datetime') %>%
mutate(DO_sat = o2.at.sat.base(temp = wtr)) %>%
mutate(Flux = k*(DO_con - DO_sat)/(0.5*PhoticDepth)) %>% #g/m3/d
mutate(datetime = as.Date(datetime))
##### READ PARAMETER FILE ##################
parameters <- read.table(file = ParameterFile,header=TRUE,comment.char="#",stringsAsFactors = F)
for (i in 1:nrow(parameters)){ # assign parameters
assign(parameters[i,1],parameters[i,2])
}
freeParams <- read.table(file = FreeParFile,header=TRUE,comment.char="#",stringsAsFactors = F)
freeParamsNames <- c('DOCR_RespParam','DOCL_RespParam','BurialFactor_R','BurialFactor_L')
for (i in 1:4){ # assign parameters
assign(freeParamsNames[i],freeParams[i,1])
}
###### Run Period and Time Step Setup #####
TimeStep <- as.numeric(InputData$datetime[2]-InputData$datetime[1]) #days
steps <- nrow(InputData)
#################### OPTIMIZATION ROUTINE ############################################
pars = c(DOCR_RespParam,DOCL_RespParam,BurialFactor_R,BurialFactor_L)
DOC_DO_diff <- function(pars){
# DOC model
modeled = modelDOC(pars[1],pars[2],pars[3],pars[4])
joinDOC = inner_join(ValidationDataDOC,modeled,by='datetime')
resDOC = joinDOC$DOC - joinDOC$DOC_conc
joinDO = inner_join(ValidationDataDO,modeled,by='datetime')
resDO = joinDO$DO_con - joinDO$MetabOxygen
lengthScale = length(resDO)/length(resDOC)
return(c(resDOC,resDO/lengthScale))
}
testACF <- function(pars){
# DOC model
modeled = modelDOC(pars[1],pars[2],pars[3],pars[4])
joinDOC = inner_join(ValidationDataDOC,modeled,by='datetime')
resDOC = joinDOC$DOC - joinDOC$DOC_conc
acf(resDOC)
joinDO = inner_join(ValidationDataDO,modeled,by='datetime')
resDO = joinDO$DO_con - joinDO$MetabOxygen
acf(resDO)
}
par(mfrow=c(2,1),mgp=c(1.5,0.5,0))
testACF(pars)
## PLOTTING and GOF
# #Goodness of fit
# library(hydroGOF)
# rmse(c(joinDOC$DOC,joinDO$DO_con), c(joinDOC$DOC_conc,joinDO$MetabOxygen.oxy_conc)) #Harp 0.43 Trout 0.427 Monona 0.62 Vanern 0.32
# NSE(c(joinDOC$DOC,joinDO$DO_con), c(joinDOC$DOC_conc,joinDO$MetabOxygen.oxy_conc)) #Harp 0.09 Trtou -0.015 Monona 0.279 Vanern -0.03
# Starting parameters cannot be negative, because of bounds we set
parStart = pars
lowerBound = c(0.00003,0.0003,0,0)
upperBound = c(0.03,3,1,1)
parStart[(parStart - lowerBound) < 0] = lowerBound[(parStart - lowerBound) < 0]
parStart[(upperBound - parStart) < 0] = upperBound[(upperBound - parStart) < 0]
names(parStart) = c('DOCR_RespParam','DOCL_RespParam','BurialFactor_R','BurialFactor_L')
# For difficult problems it may be efficient to perform some iterations with Pseudo, which will bring the algorithm
# near the vicinity of a (the) minimum, after which the default algorithm (Marq) is used to locate the minimum more precisely.
Fit2 <- modFit(f = DOC_DO_diff, p=parStart,method = 'Pseudo',
lower= lowerBound,
upper= upperBound)
Fit2par = Fit2$par # 0.001612197 0.018551566 0.999999949 0.000000000 (Trout)
Fit2par
pars
fitTest(pars,plot=F)
fitTest(Fit2$par,plot=F)
#Test Fits
fitTest <- function(pars,plot=F){
# DOC model
modeled = modelDOC(pars[1],pars[2],pars[3],pars[4])
joinDOC = inner_join(ValidationDataDOC,modeled,by='datetime')
joinDO = inner_join(ValidationDataDO,modeled,by='datetime')
# PLOTTING and GOF
if (plot == T){
png(paste0('R/FMEresults/',LakeName,'FMEfit.png'),width = 6,height = 8,units = 'in',res = 300)
}
par(mar=c(3,3,3,1),mgp=c(1.5,0.5,0),mfrow=c(2,1),cex=0.8)
plot(joinDOC$datetime,joinDOC$DOC,type='o',xlab='Date',ylab = 'DOC (mg/L)',pch=16,main=LakeName)
lines(joinDOC$datetime,joinDOC$DOCwc,type='o',col='grey50',pch=16)
lines(joinDOC$datetime,joinDOC$DOC_conc,type='o',col='red3',pch=16)
legend('bottomleft',legend = c('ObsSurf','ObsWC','Mod'),col=c('black','grey50','red3'),pch=16)
plot(joinDO$datetime,joinDO$DO_con,xlab='Date',type='o',ylab = 'DO (mg/L)',pch=16,main=LakeName)
lines(joinDO$datetime,joinDO$MetabOxygen,type='o',col='red3',pch=16)
if (plot == T){
dev.off()
}
#Goodness of fit
library(hydroGOF)
print(paste('RMSE = ',rmse(c(joinDOC$DOC,joinDO$DO_con), c(joinDOC$DOC_conc,joinDO$MetabOxygen))))
print(paste('NSE = ',NSE(c(joinDOC$DOC,joinDO$DO_con), c(joinDOC$DOC_conc,joinDO$MetabOxygen))))
}
fitTest(pars,plot=F)
fitTest(Fit2$par,plot=F)
Fit2$par
pars
pars[2]
pars[3] = 0
fitTest(pars,plot=F)
pars[2] = 0.4
fitTest(pars,plot=F)
pars[4] = 1
fitTest(pars,plot=F)
pars[4] = 0
pars[1] = 0.001
fitTest(pars,plot=F)
pars[1] = 0.003
fitTest(pars,plot=F)