shypEstFun: Wrapper function for the Estimation of Soil Hydrologic...
In spsh: Estimation and Prediction of Parameters of Various Soil Hydraulic Property Models
Description
Usage
Arguments
Details
Value
Examples
Description
Estimates model parameters of implemented soil hydraulic property functions.
This function sets up the parameter estimation, given a set of arguments, and enables minimisation of (weighted)
sum of squared residuals, assuming independent and identically distributed model residuals.
More information on the options is given in the Details
Usage
1
2
3
4
5
6
7
8
9
10
11
12
13
14
shypEstFun(
shpmodel = "01110",
parL,
retdata,
condata,
ivap = NULL,
hclip = FALSE,
weightmethod = "none",
LikModel = "rss",
ALG = "DE",
set.itermax = 200,
ALGoptions = NULL,
lhs.query = FALSE
)
Arguments
shpmodel
Character specifying the soil hydraulic property model. Currently valid models as documented in shypFun and are:
01110 constrained unimodal van Genuchten-Mualem.
01210 constrained bimodal van Genuchten-Mualem.
01310 constrained trimodal van Genuchten-Mualem.
02110 unimodel Kosugi 2 parametric-Mualem model (Kosugi, 1996)
03110 unimodel Fredlund-Xing-Mualem model, with the contraint of m = 1-1/n (Fredlund D.G., and A. Xing, 1994)
parL
list of 4 vectors with named vectors, the order in the list is sensitive.
p vector with length l of model specific initial parameters, has to coincide with the chosen soil hydraulic property model.
psel vector with length l identifying which parameters are to be estimated
plo vector of lower bounds (non-transformed parameter boundaries)
pup vector of upper bounds (non-transformed parameters boundaries)
retdata
A dataframe or matrix with 2 columns. The first with log10 values of pressure head values in [cm] and the second with volumetric water contents in [cm cm-3].
condata
A dataframe or matrix with 2 columns. The first with log10 values of pressure head values in [cm] and the second with hydraulic conductivity values log10[cm d-1].
ivap
Specification if isothermal vapour conductivity after Saito et al. (2006) is accounted, defaults to NULL and no isothermal vapour conducitvity is considered. See Details.
hclip
Implemented for future development reasons and is not yet functional. Specification if the hydraulic conductivity model should be 'clipped', i.e. constrained to a maxium pore diamater as introduced by Iden et al. (2015), defaults to FALSE.
weightmethod
Specification of weight method. The implemented methods are
none no weights are considered, i.e. no measurement error assumed
range normalization of observations to the intervall [0,1]
fix1 fixed scalar weight for THETA is 0.05^2 and weight for log10K is 1
est1 Two scalar model weights (\$1/sigma^2) are treated as free parameters to be estimated by inversion, one for THETA and one for log10K
Alternatively, a list of vectors can be provided specifying the user given model weights (\$1/sigma^2). Either as skalar for each data class, or a vector with the same length as the number of data points given for each of the measurements in the respective data class.
The length of the list has to coincide with the data groups.
LikModel
Specification of inverse modelling type. Has to be specified but implemented for future compatability)
rss Default for the optimisation algorithm DEoptim. resFun returns skalar sum of squared (weighted) residuals
-2loglik Specified if ALG == -2*log-likelihood value, which is minimised assuming Gaussian and i.i.d (weighted) residuals
ALG
Select global optimisation algorithm or a Markov chain Monte Carlos (MCMC) sampler.
DE Default for the optimisation algortihm DEoptim. resFun returns a skalar sum of squared (weighted) residuals if LikModel == "rss".
modMCMC Default for the DRAM (Delayed Rejection Adaption Metropolis) algrothim implemented in modMCMC of the FME package.
resFun returns a skalar -2loglik and LikModel = "-2logLik" has to be specified.
set.itermax
Integer specifying the maximum number of iterations default = 200.
ALGoptions
A list with named entries setting the algorithm options. Each list element name is required to be identical with the names
as documented in the respective algortihm help DEoptim.control and modMCMC.
set.itermax overrides the maximum iterations argument.
lhs.query
default FALSE, TRUE will produce a Latin Hypercube Sample for the initial population when using DEoptim.
Details
Several in-built methods for weighting the (multi-) objective function residuals are available, they may be specified, or estimated as nuisance parameters for the two data groups. More details see weightFun.
Weights are the inverse of the squared standard deviation of the residuals (variance).
Generally, soil hydraulic property model parameters are estimated as transformed parameters: log10 for alpha_i, Ks, and log10 for n_i-1, Kc, Knc
For model codes in ivap please refer to KvapFun.
Parallel computing for package DEoptim is not supported. And the optional arguments in modMCMC are not supported.
Value
list returns the result of the optimisation algrorithm or MCMC sampler and all settings.
settings
a list with output of the optimisation and summary of settings:
weigth the list with weights for the retention and conductivity data.
parL the list of initial and selected model parameters, and upper and lower bounds.
transL list of parameter transformation rules used
shpmodel the used soil hydraulic property model
ivap isothermal vapour conductivity model
hclip for future compatability
LikModel the adopted method to calculate the objective function value
data a list of 2 objects with a) retention data and b) conductivity data used for the parameter estimation.
out
result of algorithm function DEoptim or modMCMC
Examples
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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
## Not run:
data("shpdata1")
retdata <- shpdata1$TS1$wrc
condata <- shpdata1$TS1$hcc
condata <- condata[!is.na(condata[,1]),]
weightmethod <- "range"
ivap <- NULL
set.itermax <- 1
LikModel <- "rss" # ALTERNATIVE OPTION: LikModel = "-2logLik"
ALG <- "DE" # ALTERNATIVE OPTION: ALG = "modMCMC"
parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ks"=100,"tau"=.5),
"psel" = c(1, 1, 1, 1, 1, 1),
"plo"= c(0.001 , 0.2 , 0.001 , 1.1, 1, -2),
"pup"= c(0.3 , 0.8 , .1, 11 , 1e4, 10))
out <- shypEstFun(shpmodel = "01110",
parL = parL,
retdata = retdata, condata = condata,
ivap = ivap,
hclip = FALSE,
weightmethod = weightmethod,
LikModel = LikModel,
ALG = ALG,
set.itermax = set.itermax,
lhs.query = FALSE)
\dontshow{
}
\donttest{
data("shpdata1")
retdata <- ret <- shpdata1$TS1$wrc
condata <- con <- shpdata1$TS1$hcc
condata <- condata[!is.na(condata[,1]),]
---
# 1 SET VARIABLES --------------------
# VARIABLES FOR PLOTTING
{pF <- seq(-3, 6.8, length = 201)
h <- 10^pF
ticksatmin <- -2
tcllen <- 0.4
ticksat <- seq(ticksatmin,5,1)
pow <- ticksatmin:6
# VARIABLES FOR THE FITTING ALGORITHM
weightmethod = "range"
ivap = NULL
set.itermax = 3e1
LikModel = "rss" # ALTERNATIVE OPTION: LikModel = "-2logLik"
ALG = "DE" # ALTERNATIVE OPTION: ALG = "modMCMC"
shpmodel.v <- c("01110", "01110FM")
plot.query = FALSE
no.shps <- length(shpmodel.v)
# initialising lists
out.L <- vector("list", no.shps)
gof.L <- vector("list", no.shps)
}
# Run comparison
for (i in 1:2) {
shpmodel = shpmodel.v[i]
# INITIAL PARAMETERS, BOUNDS, and SELECTED PARAMETERS FOR FITTING
switch(shpmodel,
"01110" = {
# van Genuchten-Mualem Model parameters
parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ks"=100,"tau"=.5),
"psel" = c(1, 1, 1, 1, 1, 1),
"plo"= c(0.001 , 0.2 , 0.001 , 1.1, 1, -2),
"pup"= c(0.3 , 0.8 , .1, 11 , 1e4, 10)
)
},
"01110FM" = {
# van Genuchten-Mualem Model parameters + BRUNSWICK MODEL
parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ksc"=100,
"tau"=.5,"Ksnc"=1e-4,"a"=1.5,"h0"=6.8),
"psel" = c( 1,1, 1 ,1 , 1,1,1, 0, 0),
"plo"= c(0.001 , 0.1 , 0.001 , 1.1, 1,0,1e-6 , 1, 6.5),
"pup"= c(0.35, 0.7 , .1, 11 , 1e4,10 ,1e0, 2, 6.9)
)
},
stop("Enter a meaningful shpmodel")
)
out <- shypEstFun(shpmodel = shpmodel,
parL = parL,
retdata = retdata, condata = condata,
ivap = ivap,
hclip = FALSE,
weightmethod = weightmethod,
LikModel = LikModel,
ALG = ALG,
set.itermax = set.itermax
,lhs.query = FALSE)
out$model <- shpmodel.v[i]
out.L[[i]] <- out
# Calculate the soil hydraulic properties for the plot
if(ALG == "DE"){
p <- out$out$optim$phattrans
}
if(ALG == "modMCMC"){
p <- out$out$bestpartrans
}
if(weightmethod =="est1"){
np <- length(p)
p <- p[-c(np-1, np)]
if(ALG =="modMCMC"){
parL$p[which(parL$psel==1)] <- p
p <- parL$p
}
}
if(plot.query==TRUE){
shyp.L<-shypFun(p,h,shpmodel=shpmodel.v[i],ivap.query=ivap)
if(shpmodel == c("01110")){
wrc<-shyp.L$theta
hcc<-log10(shyp.L$Kh)
# PLOT THE WATER RETENTION CURVE
par(mfrow=c(1,2),tcl=tcllen)
plot(retdata,ylim=c(.0,.50),xlim=c(0,6.8),ylab="",xlab="",
col="darkgrey",axes=FALSE,main="WaterRetentionCurve",cex=2)
lines(log10(abs(h)),wrc,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2,tcl=tcllen)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
mtext("vol.watercontent[-]",2,cex=1.2,line=2.8)
box()
# PLOT THE MEASURED HYDRAULIC CONDUCTIVITY CURVE
plot(condata,ylim=c(-8,2),xlim=c(0,6.8),ylab="",xlab="",col="darkgrey",
axes=FALSE,main="HydraulicConductivityCurve",cex=2)
lines(log10(abs(h)),hcc,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("log10K[cm/d]",2,cex=1.2,line=2.8)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
box()
par(mfrow=c(1,1))
}
if(shpmodel == "01110FM"){
wrc<-shyp.L$theta
wrccap<-shyp.L$thetacap
wrcnc<-shyp.L$thetanc
hcc<-log10(shyp.L$Kh)
hcccap<-log10(shyp.L$Kcap)
hccnc<-log10(shyp.L$Knc)
hcvap<-log10(shyp.L$Kvap)
par(mfrow=c(1,2),tcl=tcllen)
plot(retdata,ylim=c(.0,.50),xlim=c(0,6.8),ylab="",xlab="",
col="darkgrey",axes=FALSE,main="WaterRetentionCurve",cex=2)
lines(log10(h),wrccap,col="brown",lwd=2)
lines(log10(h),wrcnc,col="brown",lwd=2,lty=2)
lines(log10(h),wrc,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2,tcl=tcllen)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
mtext("vol.watercontent[-]",2,cex=1.2,line=2.8)
box()
# PLOT THE HYDRAULIC CONDUCTIVITY CURVE
plot(condata,ylim=c(-8,max(max(condata[,1]),max(hcc)))
,xlim=c(0,6.8),ylab="",xlab="",col="darkgrey",
axes=FALSE,main="HydraulicConductivityCurve",cex=2)
lines(log10(h),hcccap,col="brown",lwd=2)
lines(log10(h),hccnc,col="brown",lwd=2,lty=2)
lines(log10(h),hcc,col="darkblue",lwd=2)
lines(log10(h),hcvap,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("log10K[cm/d]",2,cex=1.2,line=2.8)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
box()
par(mfrow=c(1,1))
}
}
phattrans.m <- out$out$optim$phattrans
gof.L[[i]]<-gofFun(phattrans.m,shpmodel=shpmodel.v[i],retdata=retdata,condata=condata,
out.L[[i]]$settings$weight,parL$psel,ivap.query=NULL,hclip.query=FALSE)
}
statstab3 <- cbind("th_rmse" = unlist(lapply(lapply(gof.L, `[[`, "th"), '[[', "rmse")),
"log10Kh_rmse" = unlist(lapply(lapply(gof.L, `[[`, "log10Kh"), '[[', "rmse"))
)
}
## End(Not run)
spsh documentation built on April 14, 2020, 6:37 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Examples
Description
Estimates model parameters of implemented soil hydraulic property functions. This function sets up the parameter estimation, given a set of arguments, and enables minimisation of (weighted) sum of squared residuals, assuming independent and identically distributed model residuals. More information on the options is given in the Details
Usage
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | shypEstFun(
shpmodel = "01110",
parL,
retdata,
condata,
ivap = NULL,
hclip = FALSE,
weightmethod = "none",
LikModel = "rss",
ALG = "DE",
set.itermax = 200,
ALGoptions = NULL,
lhs.query = FALSE
)
|
Arguments
shpmodel |
| |||||||||
parL |
| |||||||||
retdata |
A dataframe or matrix with 2 columns. The first with log10 values of pressure head values in [cm] and the second with volumetric water contents in [cm cm-3]. | |||||||||
condata |
A dataframe or matrix with 2 columns. The first with log10 values of pressure head values in [cm] and the second with hydraulic conductivity values log10[cm d-1]. | |||||||||
ivap |
Specification if isothermal vapour conductivity after Saito et al. (2006) is accounted, defaults to | |||||||||
hclip |
Implemented for future development reasons and is not yet functional. Specification if the hydraulic conductivity model should be 'clipped', i.e. constrained to a maxium pore diamater as introduced by Iden et al. (2015), defaults to | |||||||||
weightmethod |
Specification of weight method. The implemented methods are
Alternatively, a list of vectors can be provided specifying the user given model weights (\$1/sigma^2). Either as skalar for each data class, or a vector with the same length as the number of data points given for each of the measurements in the respective data class. The length of the list has to coincide with the data groups. | |||||||||
LikModel |
Specification of inverse modelling type. Has to be specified but implemented for future compatability)
| |||||||||
ALG |
Select global optimisation algorithm or a Markov chain Monte Carlos (MCMC) sampler.
| |||||||||
set.itermax |
Integer specifying the maximum number of iterations | |||||||||
ALGoptions |
A list with named entries setting the algorithm options. Each list element name is required to be identical with the names
as documented in the respective algortihm help DEoptim.control and modMCMC. | |||||||||
lhs.query |
default |
Details
Several in-built methods for weighting the (multi-) objective function residuals are available, they may be specified, or estimated as nuisance parameters for the two data groups. More details see weightFun.
Weights are the inverse of the squared standard deviation of the residuals (variance).
Generally, soil hydraulic property model parameters are estimated as transformed parameters: log10 for alpha_i, Ks, and log10 for n_i-1, Kc, Knc
For model codes in ivap please refer to KvapFun.
Parallel computing for package DEoptim is not supported. And the optional arguments in modMCMC are not supported.
Value
list returns the result of the optimisation algrorithm or MCMC sampler and all settings.
settings |
a
| ||||||||||||||||
out |
result of algorithm function |
Examples
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 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 | ## Not run:
data("shpdata1")
retdata <- shpdata1$TS1$wrc
condata <- shpdata1$TS1$hcc
condata <- condata[!is.na(condata[,1]),]
weightmethod <- "range"
ivap <- NULL
set.itermax <- 1
LikModel <- "rss" # ALTERNATIVE OPTION: LikModel = "-2logLik"
ALG <- "DE" # ALTERNATIVE OPTION: ALG = "modMCMC"
parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ks"=100,"tau"=.5),
"psel" = c(1, 1, 1, 1, 1, 1),
"plo"= c(0.001 , 0.2 , 0.001 , 1.1, 1, -2),
"pup"= c(0.3 , 0.8 , .1, 11 , 1e4, 10))
out <- shypEstFun(shpmodel = "01110",
parL = parL,
retdata = retdata, condata = condata,
ivap = ivap,
hclip = FALSE,
weightmethod = weightmethod,
LikModel = LikModel,
ALG = ALG,
set.itermax = set.itermax,
lhs.query = FALSE)
\dontshow{
}
\donttest{
data("shpdata1")
retdata <- ret <- shpdata1$TS1$wrc
condata <- con <- shpdata1$TS1$hcc
condata <- condata[!is.na(condata[,1]),]
---
# 1 SET VARIABLES --------------------
# VARIABLES FOR PLOTTING
{pF <- seq(-3, 6.8, length = 201)
h <- 10^pF
ticksatmin <- -2
tcllen <- 0.4
ticksat <- seq(ticksatmin,5,1)
pow <- ticksatmin:6
# VARIABLES FOR THE FITTING ALGORITHM
weightmethod = "range"
ivap = NULL
set.itermax = 3e1
LikModel = "rss" # ALTERNATIVE OPTION: LikModel = "-2logLik"
ALG = "DE" # ALTERNATIVE OPTION: ALG = "modMCMC"
shpmodel.v <- c("01110", "01110FM")
plot.query = FALSE
no.shps <- length(shpmodel.v)
# initialising lists
out.L <- vector("list", no.shps)
gof.L <- vector("list", no.shps)
}
# Run comparison
for (i in 1:2) {
shpmodel = shpmodel.v[i]
# INITIAL PARAMETERS, BOUNDS, and SELECTED PARAMETERS FOR FITTING
switch(shpmodel,
"01110" = {
# van Genuchten-Mualem Model parameters
parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ks"=100,"tau"=.5),
"psel" = c(1, 1, 1, 1, 1, 1),
"plo"= c(0.001 , 0.2 , 0.001 , 1.1, 1, -2),
"pup"= c(0.3 , 0.8 , .1, 11 , 1e4, 10)
)
},
"01110FM" = {
# van Genuchten-Mualem Model parameters + BRUNSWICK MODEL
parL<-list("p"=c("thr"=0.05,"ths"=0.45,"alf1"=0.01,"n"=2,"Ksc"=100,
"tau"=.5,"Ksnc"=1e-4,"a"=1.5,"h0"=6.8),
"psel" = c( 1,1, 1 ,1 , 1,1,1, 0, 0),
"plo"= c(0.001 , 0.1 , 0.001 , 1.1, 1,0,1e-6 , 1, 6.5),
"pup"= c(0.35, 0.7 , .1, 11 , 1e4,10 ,1e0, 2, 6.9)
)
},
stop("Enter a meaningful shpmodel")
)
out <- shypEstFun(shpmodel = shpmodel,
parL = parL,
retdata = retdata, condata = condata,
ivap = ivap,
hclip = FALSE,
weightmethod = weightmethod,
LikModel = LikModel,
ALG = ALG,
set.itermax = set.itermax
,lhs.query = FALSE)
out$model <- shpmodel.v[i]
out.L[[i]] <- out
# Calculate the soil hydraulic properties for the plot
if(ALG == "DE"){
p <- out$out$optim$phattrans
}
if(ALG == "modMCMC"){
p <- out$out$bestpartrans
}
if(weightmethod =="est1"){
np <- length(p)
p <- p[-c(np-1, np)]
if(ALG =="modMCMC"){
parL$p[which(parL$psel==1)] <- p
p <- parL$p
}
}
if(plot.query==TRUE){
shyp.L<-shypFun(p,h,shpmodel=shpmodel.v[i],ivap.query=ivap)
if(shpmodel == c("01110")){
wrc<-shyp.L$theta
hcc<-log10(shyp.L$Kh)
# PLOT THE WATER RETENTION CURVE
par(mfrow=c(1,2),tcl=tcllen)
plot(retdata,ylim=c(.0,.50),xlim=c(0,6.8),ylab="",xlab="",
col="darkgrey",axes=FALSE,main="WaterRetentionCurve",cex=2)
lines(log10(abs(h)),wrc,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2,tcl=tcllen)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
mtext("vol.watercontent[-]",2,cex=1.2,line=2.8)
box()
# PLOT THE MEASURED HYDRAULIC CONDUCTIVITY CURVE
plot(condata,ylim=c(-8,2),xlim=c(0,6.8),ylab="",xlab="",col="darkgrey",
axes=FALSE,main="HydraulicConductivityCurve",cex=2)
lines(log10(abs(h)),hcc,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("log10K[cm/d]",2,cex=1.2,line=2.8)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
box()
par(mfrow=c(1,1))
}
if(shpmodel == "01110FM"){
wrc<-shyp.L$theta
wrccap<-shyp.L$thetacap
wrcnc<-shyp.L$thetanc
hcc<-log10(shyp.L$Kh)
hcccap<-log10(shyp.L$Kcap)
hccnc<-log10(shyp.L$Knc)
hcvap<-log10(shyp.L$Kvap)
par(mfrow=c(1,2),tcl=tcllen)
plot(retdata,ylim=c(.0,.50),xlim=c(0,6.8),ylab="",xlab="",
col="darkgrey",axes=FALSE,main="WaterRetentionCurve",cex=2)
lines(log10(h),wrccap,col="brown",lwd=2)
lines(log10(h),wrcnc,col="brown",lwd=2,lty=2)
lines(log10(h),wrc,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2,tcl=tcllen)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
mtext("vol.watercontent[-]",2,cex=1.2,line=2.8)
box()
# PLOT THE HYDRAULIC CONDUCTIVITY CURVE
plot(condata,ylim=c(-8,max(max(condata[,1]),max(hcc)))
,xlim=c(0,6.8),ylab="",xlab="",col="darkgrey",
axes=FALSE,main="HydraulicConductivityCurve",cex=2)
lines(log10(h),hcccap,col="brown",lwd=2)
lines(log10(h),hccnc,col="brown",lwd=2,lty=2)
lines(log10(h),hcc,col="darkblue",lwd=2)
lines(log10(h),hcvap,col="darkblue",lwd=2)
legend("topright",c("observed","simulated"),pch=c(1,NA),
lty=c(NA,1),lwd=2,bty="n",cex=1.3,col=c("darkgrey","darkblue"))
axis(1,at=pow,labels=parse(text=paste('10^',(pow),sep="")),tcl=tcllen)
axis(2)
axis(4,labels=NA)
axis(3,labels=NA)
mtext("log10K[cm/d]",2,cex=1.2,line=2.8)
mtext("pressurehead|h|[cm]",1,cex=1.2,line=2.8)
box()
par(mfrow=c(1,1))
}
}
phattrans.m <- out$out$optim$phattrans
gof.L[[i]]<-gofFun(phattrans.m,shpmodel=shpmodel.v[i],retdata=retdata,condata=condata,
out.L[[i]]$settings$weight,parL$psel,ivap.query=NULL,hclip.query=FALSE)
}
statstab3 <- cbind("th_rmse" = unlist(lapply(lapply(gof.L, `[[`, "th"), '[[', "rmse")),
"log10Kh_rmse" = unlist(lapply(lapply(gof.L, `[[`, "log10Kh"), '[[', "rmse"))
)
}
## End(Not run)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.