| shypEstFun {spsh} | R Documentation |
Wrapper function for the Estimation of Soil Hydrologic Property Model Parameters
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
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
## 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)