Last active
August 29, 2015 14:16
-
-
Save ougx/45bbfd8fce7b6a90ac19 to your computer and use it in GitHub Desktop.
Create the SWAT files for PEST
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # Produced by Michael Ou, 2/11/2015 | |
| # using PEST for SWAT calibration | |
| # grouping the same parameter for all HRUs | |
| # change of parameter is accommodated by multiplying a relative rate which is operated by PEST | |
| # for uniform parameters (same value for all HRUs), we replace it directly. | |
| # observations are donwloaded automatically from USGS | |
| # r__CN2.mgt -0.2 0.2 | |
| # v__ALPHA_BF.gw 0.0 1.0 | |
| # v__GW_DELAY.gw 30.0 450.0 | |
| # v__GWQMN.gw 0.0 2.0 | |
| # | |
| # | |
| # v__GW_REVAP.gw 0.0 0.2 | |
| # v__ESCO.hru 0.8 1.0 | |
| # v__CH_N2.rte 0.0 0.3 | |
| # v__CH_K2.rte 5.0 130.0 | |
| # v__ALPHA_BNK.rte 0.0 1.0 | |
| # r__SOL_AWC(1).sol -0.2 0.4 | |
| # r__SOL_K(1).sol -0.8 0.8 | |
| # r__SOL_BD(1).sol -0.5 0.6 | |
| # v__SFTMP.bsn -5.0 5.0 | |
| # | |
| # | |
| # | |
| # -------------------------------------------------- | |
| # example of parameterization: | |
| # | |
| # r__Precipitation(){1990001-2000265}.pcp 0 0.4 | |
| # r__AUTO_EFF{[],11,AUTO_NSTRS=0.6}.mgt 0.2 0.6 | |
| # r__BLAI{22}.crop.dat 2 8 | |
| # | |
| # | |
| # | |
| # ----------------------------------------- | |
| # | |
| # r__CN2.mgt________1 -0.2 0.2 | |
| # r__SOL_AWC(1).sol________1 -0.2 0.1 | |
| # r__SOL_K(1).sol________1 -0.8 0.8 | |
| # r__SOL_BD(1).sol________1 -0.5 0.6 | |
| # a__GWQMN.gw________1 0.0 25.0 | |
| # a__GW_REVAP.gw________1 -0.1 0.0 | |
| # v__REVAPMN.gw________1 0.0 10.0 | |
| # a__ESCO.hru________1 0.0 0.2 | |
| # r__HRU_SLP.hru________1 0.0 0.2 | |
| # r__OV_N.hru________1 -0.2 0.0 | |
| # r__SLSUBBSN.hru________1 0.0 0.2 | |
| # | |
| # | |
| # | |
| # r__CN2.mgt________3-6 -0.2 0.2 | |
| # r__SOL_AWC(1).sol________3-6 -0.2 0.1 | |
| # r__SOL_K(1).sol________3-6 -0.8 0.8 | |
| # r__SOL_BD(1).sol________3-6 -0.5 0.6 | |
| # a__GWQMN.gw________3-6 0.0 25.0 | |
| # a__GW_REVAP.gw________3-6 -0.1 0.0 | |
| # v__REVAPMN.gw________3-6 0.0 10.0 | |
| # a__ESCO.hru________3-6 0.0 0.2 | |
| # r__HRU_SLP.hru________3-6 0.0 0.2 | |
| # r__OV_N.hru________3-6 -0.2 0.0 | |
| # r__SLSUBBSN.hru________3-6 0.0 0.2 | |
| # | |
| # | |
| # r__CN2.mgt________7-20 -0.2 0.2 | |
| # r__SOL_AWC(1).sol________7-20 -0.2 0.1 | |
| # r__SOL_K(1).sol________7-20 -0.8 0.8 | |
| # r__SOL_BD(1).sol________7-20 -0.5 0.6 | |
| # a__GWQMN.gw________7-20 0.0 25.0 | |
| # a__GW_REVAP.gw________7-20 -0.1 0.0 | |
| # v__REVAPMN.gw________7-20 0.0 10.0 | |
| # a__ESCO.hru________7-20 0.0 0.2 | |
| # r__HRU_SLP.hru________7-20 0.0 0.2 | |
| # r__OV_N.hru________7-20 -0.2 0.0 | |
| # r__SLSUBBSN.hru________7-20 0.0 0.2 | |
| # | |
| # | |
| # v__ALPHA_BF.gw 0.0 1.0 | |
| # v__GW_DELAY.gw 30.0 450.0 | |
| # v__CH_N2.rte 0.0 0.3 | |
| # v__CH_K2.rte 5.0 130.0 | |
| # v__ALPHA_BNK.rte 0.0 1.0 | |
| # v__SFTMP.bsn -5.0 5.0 | |
| library("stringr") | |
| library("gdata") | |
| setwd("D:\\@Workspace\\Okalahoma_Lake_Cr\\CN_46_HRU\\PESTCN") | |
| NPAR = 0 | |
| NOBS = 0 | |
| NPARGP = 0 | |
| NPRIOR = 0 | |
| NOBSGP = 0 | |
| # PARAMETER GROUP | |
| # PARGPNME INCTYP DERINC DERINCLB FORCEN DERINCMUL DERMTHD | |
| #123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc | |
| # 1 2 3 4 5 6 7 8 9 10 | |
| PARGPNME <- c( "SLSUBBSN", "HRU_SLP", "OV_N", "ESCO", "CN2", "BD", "AWC", "KSAT", "GW_DELAY", "ALPHA_BF", | |
| "GWQMN", "RCHRG_DP", "GW_SPYLD", "KS", "QR", "QS", "LAMBDA", "HBUB", "WCET") | |
| NPARGP = length(PARGPNME) | |
| #INCTYP is a character variable which can assume the values "relative", "absolute" or "rel_to_max". | |
| INCTYP<-rep("relative", NPARGP) | |
| #DERINC is the increment used for forward-difference calculation of derivatives | |
| DERINC <- c( 0.1, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.5, 0.1, | |
| 0.1, 0.1, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05) | |
| #DERINCLB is an absolute lower bound can be placed on parameter increments | |
| DERINCLB <-rep(0.01, NPARGP) | |
| #The character variable FORCEN (an abbreviation of "FORward/CENtral") determines whether derivatives for group members are calculated using | |
| # forward differences, one of the variants of the central difference method, of whether both alternatives are used in the course of an optimisation run. | |
| # It must assume one of the values "always_2", "always_3" or "switch". | |
| # Experience has shown that in most instances the most appropriate value for FORCEN is "switch". | |
| FORCEN = rep("switch", NPARGP) | |
| #DERINCMUL Whenever the central method is employed for derivatives calculation, DERINC is multiplied by DERINCMUL | |
| # Experience shows that a value between 1.0 and 2.0 is usually satisfactory. | |
| DERINCMUL = rep(2.0, NPARGP) | |
| # | |
| DERMTHD = rep("parabolic", NPARGP) | |
| ParaGroup <- data.frame(PARGPNME,INCTYP,DERINC,DERINCLB,FORCEN,DERINCMUL,DERMTHD) | |
| row.names(ParaGroup) <- PARGPNME | |
| # Parameter Data | |
| # PARNME PARTRANS PARCHGLIM PARVAL1 PARLBND PARUBND PARGP SCALE OFFSET DERCOM | |
| #123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc | |
| # 1 2 3 4 5 6 7 8 9 10 | |
| PARNME <- PARGPNME | |
| NPAR = length(PARNME) | |
| #PARTRANS is a character variable which must assume one of four values, viz. "none", "log", "fixed" or "tied". | |
| PARTRANS<-rep("none", NPAR) | |
| #This character variable is used to designate whether an adjustable parameter is relative-limited or factor-limited; | |
| # PARCHGLIM must be provided with one of two possible values, viz. "relative" or "factor". | |
| PARCHGLIM<-rep("factor", NPAR) | |
| #PARVAL1, a real variable, is a parameter's initial value. | |
| PARVAL1 <- rep(1, NPAR) | |
| #PARLBND and PARUBND: These two real variables represent a parameter's lower and upper bounds respectively. | |
| PARLBND <- c(0.01, 0.05, 0.1, 0.05, 0.8, 0.1, 0.3, 0.001, 0.1, 0.1, | |
| 0.01, 0.1, 0.05, 0.0001, 0.1, 0.1, 0.1, 0.01, 0.1) | |
| PARUBND <-c( 10.0, 10.0, 10.0, 10.0, 1.3, 3.0, 3.0, 1.0e3, 20, 10, | |
| 1e5, 10.0, 20.0, 1e4, 10, 10, 10, 100, 3) | |
| #PARGP is the name of the group to which a parameter belongs. | |
| PARGP = PARGPNME | |
| #SCALE and OFFSET | |
| #Just before a parameter value is written to a model input file (be it for initial determination of the objective function, | |
| # derivatives calculation or parameter upgrade), it is multiplied by the real variable SCALE, after which the real variable OFFSET is added. | |
| SCALE = rep(1.0, NPAR) | |
| OFFSET = rep(0.0, NPAR) | |
| #DERCOM: Unless using PEST's external derivatives functionality (see Chapter 9), this variable should be set to 1. | |
| DERCOM = rep(1, NPAR) | |
| # PARNME PARTRANS PARCHGLIM PARVAL1 PARLBND PARUBND PARGP SCALE OFFSET DERCOM | |
| ParaData <- data.frame(PARNME, PARTRANS, PARCHGLIM, PARVAL1, PARLBND, PARUBND, PARGP, SCALE, OFFSET, DERCOM) | |
| row.names(ParaData) <- PARNME | |
| # real values | |
| #123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc | |
| # 1 2 3 4 5 6 7 8 9 10 | |
| # "SLSUBBSN", "HRU_SLP", "OV_N", "ESCO", "CN2", "BD", "AWC", "KSAT", "GW_DELAY", "ALPHA_BF", | |
| # "GWQMN", "RCHRG_DP", "GW_SPYLD", "KS", "QR", "QS", "LAMBDA", "HBUB", "WCET" | |
| REALLBND <-c(10.0, 0.001, 0.008, 0.8, 35, 0.9, 0.05, 1e-4, 1, 0.01, | |
| 1, 0.01, 0.01, 1e-4, 1e-4, 0.1, 0.04, 0.1, 0.01) | |
| REALUBND <-c(1000, 0.3, 0.6, 1.0, 98, 2.0, 0.5, 1e4, 1000, 0.99, | |
| 1e5, 0.99, 0.5, 5e3, 0.2, 0.6, 1.2, 5000, 0.3) | |
| RealPara <- data.frame(PARNME, REALLBND, REALUBND) | |
| row.names(RealPara) <- PARNME | |
| #read number of sub | |
| nsub = 0 | |
| pname <- "" | |
| pgroup <- "" | |
| pvalue <- 0.0 | |
| pline <- "" | |
| tplname <- "" #this is the template file for par2par but not the pest | |
| ntpl <- 0 | |
| np <- 0 | |
| filename<-paste0("fig.fig") | |
| fig <- paste(readLines(filename), sep = "\n") | |
| fig <- chartr("\t"," ", fig) #delete tab | |
| ZBOT = 0 | |
| BD = 0 | |
| AWC = 0 | |
| KSAT = 0 | |
| isub = 0 | |
| ihru = 0 | |
| i = 0 | |
| hrutot = 0 | |
| for (isub in 1:length(fig)){ | |
| sline <- unlist(strsplit(fig[isub]," ")) | |
| sline <- sline[sline != ""] | |
| if (sline[1] == "subbasin") { | |
| nsub = nsub + 1 | |
| sline <- unlist(strsplit(fig[isub+1]," ")) | |
| sline <- sline[sline != ""] | |
| #### sub file ### | |
| filename<-sline[1] | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| sub <- paste(readLines(filename), sep = "\n") | |
| #CH_N1 | |
| #np=np+1 | |
| #pvalue[np]<-as.numeric(substr(sub[29], 1, 16)) | |
| #pname[np]<- paste0("CH_N1",isub) | |
| #pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.008 ","\t0.6","\tCH_N1 "," \t1.0 \t0.0 \t1") | |
| #substr(sub[29], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", sub), tplname[ntpl], sep = "\n") | |
| nhru = as.numeric(substr(sub[53], 1, 16)) | |
| hrutot = hrutot + nhru | |
| iisub = nsub | |
| for (ihru in 1:nhru){ | |
| ### hru file ### | |
| filename<-substr(sub[61+ihru], 1, 13) | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| hru <- paste(readLines(filename), sep = "\n") | |
| #SLSUBBSN | |
| np=np+1 | |
| pgroup[np] <- "SLSUBBSN" | |
| pvalue[np]<-as.numeric(substr(hru[3], 1, 16)) | |
| pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t10.0 ","\t1000.0","\t",pgroup[np]," \t1.0 \t0.0 \t1") | |
| substr(hru[3], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #HRU_SLP | |
| np=np+1 | |
| pgroup[np] <- "HRU_SLP" | |
| pvalue[np]<-as.numeric(substr(hru[4], 1, 16)) | |
| pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.005 ","\t0.2","\t",pgroup[np]," \t1.0 \t0.0 \t1") | |
| substr(hru[4], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #OV_N | |
| np=np+1 | |
| pgroup[np] <- "OV_N" | |
| pvalue[np]<-as.numeric(substr(hru[5], 1, 16)) | |
| pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.008 ","\t0.6","\t",pgroup[np]," \t1.0 \t0.0 \t1") | |
| substr(hru[5], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #ESCO | |
| np=np+1 | |
| pgroup[np] <- "ESCO" | |
| pvalue[np]<-as.numeric(substr(hru[11], 1, 16)) | |
| pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.8 ","\t1.0 ","\tESCO "," \t1.0 \t0.0 \t1") | |
| substr(hru[11], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| # #SURLAG | |
| # np=np+1 | |
| # pvalue[np]<-as.numeric(substr(hru[11], 1, 16)) | |
| # pname[np]<- paste0("SURLAG",iisub,ihru) | |
| # pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0 ","\t0.6","\tSURLAG "," \t1.0 \t0.0 \t1") | |
| # substr(hru[5], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", hru), tplname[ntpl], sep = "\n") | |
| ### mgt ### | |
| filename<-substr(sub[61+ihru], 14, 26) | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| mgt <- paste(readLines(filename), sep = "\n") | |
| #CN2 | |
| np=np+1 | |
| pgroup[np] <- "CN2" | |
| pvalue[np]<-as.numeric(substr(mgt[11], 1, 16)) | |
| pname[np]<- paste0("CN2",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t35 ","\t98","\tCN2 "," \t1.0 \t0.0 \t1") | |
| substr(mgt[11], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", mgt), tplname[ntpl], sep = "\n") | |
| ### sol ### | |
| filename<-substr(sub[61+ihru], 27, 39) | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| sol <- paste(readLines(filename), sep = "\n") | |
| NLAY <- (str_length(sol[8]) - 27) / 12 | |
| INDEX <- 28 + ((1:(NLAY+1))-1) * 12 | |
| ii <- 1:NLAY | |
| for (i in ii){ | |
| #BD | |
| np=np+1 | |
| pgroup[np] <- "BD" | |
| pvalue[np]<-as.numeric(substr(sol[9], INDEX[i], INDEX[i]+11)) | |
| pname[np]<- paste0("BD",iisub,ihru,i) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.9 ","\t2.0","\tBD "," \t1.0 \t0.0 \t1") | |
| substr(sol[9], INDEX[i], INDEX[i]+11) <- paste0("$",substr(paste0(pname[np]," "),1,10),"$") | |
| #AWC | |
| np=np+1 | |
| pgroup[np] <- "AWC" | |
| pvalue[np]<-as.numeric(substr(sol[10], INDEX[i], INDEX[i]+11)) | |
| pname[np]<- paste0("AWC",iisub,ihru,i) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.05 ","\t0.8","\tAWC "," \t1.0 \t0.0 \t1") | |
| substr(sol[10], INDEX[i], INDEX[i]+11) <- paste0("$",substr(paste0(pname[np]," "),1,10),"$") | |
| #KSAT | |
| np=np+1 | |
| pgroup[np] <- "KSAT" | |
| pvalue[np]<- as.numeric(substr(sol[11], INDEX[i], INDEX[i]+11)) | |
| pname[np]<- paste0("KSAT",iisub,ihru,i) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.0001 ","\t2000","\tKSAT "," \t1.0 \t0.0 \t1") | |
| substr(sol[11], INDEX[i], INDEX[i]+11) <- paste0("$",substr(paste0(pname[np]," "),1,10),"$") | |
| } | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", sol), tplname[ntpl], sep = "\n") | |
| # ### usf ### | |
| # filename<-paste0(substr(sub[61+ihru], 27, 35),".usf") | |
| # filename <-gsub("^\\s+|\\s+$", "", filename) | |
| # usf <- paste(readLines(filename), sep = "\n") | |
| # #skip the comment lines | |
| # | |
| # for (ncmt in 1:length(usf)){ | |
| # if (substr(usf[ncmt],1,1) != "#") break() | |
| # } | |
| # | |
| # #WCET | |
| # np=np+1 | |
| # pgroup[np] <- "WCET" | |
| # pvalue[np]<-as.numeric(substr(usf[5+ncmt], 1, 16)) | |
| # pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| # substr(usf[5+ncmt], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| # | |
| # | |
| # ntpl = ntpl + 1 | |
| # tplname[ntpl] <- paste0(filename,".tpl") | |
| # writeLines(c("ptf $", usf), tplname[ntpl], sep = "\n") | |
| ### chm ### | |
| ### gw ### | |
| filename<-substr(sub[61+ihru], 53, 65) | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| gw <- paste(readLines(filename), sep = "\n") | |
| #GW_DELAY | |
| np=np+1 | |
| pgroup[np] <- "GW_DELAY" | |
| pvalue[np]<-as.numeric(substr(gw[4], 1, 16)) | |
| pname[np]<- paste0("GW_DELAY",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t20 ","\t450","\tGW_DELAY"," \t1.0 \t0.0 \t1") | |
| substr(gw[4], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #ALPHA_BF | |
| np=np+1 | |
| pgroup[np] <- "ALPHA_BF" | |
| pvalue[np]<-as.numeric(substr(gw[5], 1, 16)) | |
| pname[np]<- paste0("ALPHA_BF",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t0.001 ","\t0.999","\tALPHA_BF"," \t1.0 \t0.0 \t1") | |
| substr(gw[5], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #GWQMN | |
| np=np+1 | |
| pgroup[np] <- "GWQMN" | |
| pvalue[np]<-as.numeric(substr(gw[6], 1, 16)) | |
| if (pvalue[np]<0.1) pvalue[np] <- 1.0 | |
| pname[np]<- paste0("GWQMN",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t0.01 ","\t10000.0","\tGWQMN"," \t1.0 \t0.0 \t1") | |
| substr(gw[6], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #RCHRG_DP | |
| np=np+1 | |
| pgroup[np] <- "RCHRG_DP" | |
| pvalue[np]<-as.numeric(substr(gw[9], 1, 16)) | |
| pname[np]<- paste0("RCHRG_DP",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t0.001 ","\t0.9","\tRCHRG_DP","\t1.0 \t0.0 \t1") | |
| substr(gw[9], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #GW_SPYLD | |
| np=np+1 | |
| pgroup[np] <- "GW_SPYLD" | |
| pvalue[np]<-as.numeric(substr(gw[11], 1, 16)) | |
| pname[np]<- paste0("GW_SPYLD",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t0.01 ","\t0.6","\tGW_SPYLD","\t1.0 \t0.0 \t1") | |
| substr(gw[11], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", gw), tplname[ntpl], sep = "\n") | |
| } | |
| } | |
| } | |
| #translate soil mat file | |
| # filename <- "soil.mat" | |
| # mat <- paste(readLines(filename), sep = "\n") | |
| # mat <- chartr("\t"," ", mat) | |
| # for (ncmt in 1:length(mat)){ | |
| # if (substr(mat[ncmt],1,1) != "#") break() | |
| # } | |
| # | |
| # sline <- unlist(strsplit(mat[ncmt]," ")) | |
| # sline <- sline[sline != ""] | |
| # | |
| # nmat <- as.numeric(sline[1]) | |
| # for (i in 1:(nmat)){ | |
| # sline <- unlist(strsplit(mat[i+ncmt]," ")) | |
| # sline <- sline[sline != ""] | |
| # | |
| # #KS | |
| # np=np+1 | |
| # pgroup[np] <- "KS" | |
| # pvalue[np]<-as.numeric(sline[1]) | |
| # pname[np]<- paste0(pgroup[np],i) | |
| # | |
| # | |
| # #QR | |
| # np=np+1 | |
| # pgroup[np] <- "QR" | |
| # pvalue[np]<-as.numeric(sline[2]) | |
| # pname[np]<- paste0(pgroup[np],i) | |
| # | |
| # | |
| # #QS | |
| # np=np+1 | |
| # pgroup[np] <- "QS" | |
| # pvalue[np]<-as.numeric(sline[3]) | |
| # pname[np]<- paste0(pgroup[np],i) | |
| # | |
| # | |
| # #LAMBDA | |
| # np=np+1 | |
| # pgroup[np] <- "LAMBDA" | |
| # pvalue[np]<-as.numeric(sline[4]) | |
| # pname[np]<- paste0(pgroup[np],i) | |
| # | |
| # | |
| # #ETA calculated by par2par | |
| # #np=np+1 | |
| # #pvalue[np]<-as.numeric(sline[5]) | |
| # #pname[np]<- paste0("LAMBDA",i) | |
| # #pline[np]<- paste0(pname[np],"\ttied "," \tfactor ","\t",pvalue[np]," \t10.0 ","\t1000.0","\tETA","\t1.0 \t0.0 \t1") | |
| # | |
| # #HBUB | |
| # np=np+1 | |
| # pgroup[np] <- "HBUB" | |
| # pvalue[np]<- -as.numeric(sline[6]) | |
| # pname[np]<- paste0(pgroup[np],i) | |
| # | |
| # mat[i+ncmt] <- paste0(paste0("$",substr(paste0(pname[np-4]," "),1,14),"$\t "), | |
| # paste0("$",substr(paste0(pname[np-3]," "),1,14),"$\t "), | |
| # paste0("$",substr(paste0(pname[np-2]," "),1,14),"$\t "), | |
| # paste0("$",substr(paste0(pname[np-1]," "),1,14),"$\t "), | |
| # paste0("$",substr(paste0("ETA",i," "),1,14),"$\t "), | |
| # paste0("$",substr(paste0(pname[np]," "),1,14),"$")) | |
| # } | |
| # | |
| # ntpl = ntpl + 1 | |
| # tplname[ntpl] <- paste0(filename,".tpl") | |
| # writeLines(c("ptf $", mat), tplname[ntpl], sep = "\n") | |
| # paragroup | |
| groups <- unique(pgroup) | |
| ngroup = length(groups) | |
| NPAR = ngroup | |
| NPARGP = ngroup | |
| ParaGroup <- subset(ParaGroup, PARGPNME %in% groups) | |
| #write.table(ParaGroup, "!!paragroup.dat", quote = FALSE, sep = " \t", row.names = FALSE, col.names = FALSE) | |
| write.fwf(ParaGroup, "!!paragroup.dat", quote = FALSE, rownames = FALSE, colnames = FALSE) | |
| # * parameter data for PEST (parameter groups) | |
| ParaData <- subset(ParaData, PARNME %in% groups) | |
| write.fwf(ParaData, "!!parameter.dat", quote = FALSE, rownames = FALSE, colnames = FALSE) | |
| #write.table(ParaData, "!!parameter.dat", quote = FALSE, row.names = FALSE, col.names = FALSE) | |
| ppestline = "" | |
| #for par2par tpl | |
| p2pline <- "" | |
| l=0 | |
| p2pline[(l=l+1)] <- "ptf $" | |
| p2pline[(l=l+1)] <- "* parameter data" | |
| for (gg in groups){ | |
| p2pline[(l=l+1)] <- paste(substr(paste0(gg," "),1,20), "=", paste0("$",substr(paste0(gg," "),1,14),"$")) | |
| } | |
| for (i in 1:np) { | |
| p2pline[(l=l+1)] <- paste(substr(paste0(pname[i]," "),1,20), "=", "max ( min (", | |
| substr(paste0(pgroup[i]," "),1,10), | |
| "*", substr(paste0(pvalue[i]," "),1,10), ",", | |
| substr(paste0(RealPara[pgroup[i],"REALUBND"]," "),1,10), "),", RealPara[pgroup[i],"REALLBND"], ")") | |
| if (pgroup[i] == "QS"){ | |
| p2pline[l] <- paste(substr(paste0(pname[i]," "),1,20), "=", "max ( min (", | |
| substr(paste0(pgroup[i]," "),1,10), | |
| "*", substr(paste0(pvalue[i]," "),1,10), ",", | |
| substr(paste0(RealPara[pgroup[i],"REALUBND"]," "),1,10), "),", pname[i-1], " + 0.1 )") | |
| } | |
| if (pgroup[i] == "LAMBDA"){ | |
| p2pline[(l=l+1)] <- paste(sprintf("%-20s", gsub("LAMBDA", "ETA", pname[i])), "= 2 /", pname[i],"+ 3") | |
| } | |
| if (pgroup[i] == "HBUB"){ | |
| p2pline[l] <- paste(substr(paste0(pname[i]," "),1,20), "=", "- max ( min (", | |
| substr(paste0(pgroup[i]," "),1,10), | |
| "*", substr(paste0(pvalue[i]," "),1,10), ",", | |
| substr(paste0(RealPara[pgroup[i],"REALUBND"]," "),1,10), "),", RealPara[pgroup[i],"REALLBND"], ")") | |
| } | |
| # if (pgroup[i] == "WCET"){ | |
| # p2pline[l] <- paste(substr(paste0(pname[i]," "),1,20), "=", "max ( min (", | |
| # substr(paste0(pgroup[i]," "),1,10), | |
| # "*", substr(paste0(pvalue[i]," "),1,10), ",", | |
| # substr(paste0(gsub("WCET", "QR", pname[i])," "),1,10), "),", gsub("WCET", "QS", pname[i]), ")") | |
| # } | |
| } | |
| p2pline[(l=l+1)] <- "* template and model input files" | |
| for (i in 1:ntpl){ | |
| p2pline[(l=l+1)] <- paste(tplname[i], substr(tplname[i], 1, str_length(tplname[i])-4)) | |
| } | |
| p2pline[(l=l+1)] <- "* control data" | |
| p2pline[(l=l+1)] <- "single point" | |
| writeLines(p2pline, "!!par2par.tpl", sep = "\n") | |
| #create observation data and instruction file | |
| nhru = hrutot | |
| hrus <- c(17,hrutot) | |
| NOBSGP = length(hrus) | |
| OBGNME <- paste0("flow", hrus) | |
| lad <- c(hrus[1],diff(hrus)) | |
| cio <- paste(readLines("file.cio"), sep = "\n") | |
| NBYR <- as.numeric(substr(cio[8], 1, 16)) | |
| IYR <- as.numeric(substr(cio[9], 1, 16)) | |
| IDAF <- as.numeric(substr(cio[10], 1, 16)) | |
| IDAL <- as.numeric(substr(cio[11], 1, 16)) | |
| NYSKIP <- as.numeric(substr(cio[60], 1, 16)) | |
| beginDate <- as.Date(IDAF-1, as.Date(paste0(IYR+NYSKIP,"-1-1"))) | |
| endDate <- as.Date(IDAL-1, as.Date(paste0(IYR+NBYR-1,"-1-1"))) | |
| ndays <- endDate - beginDate + 1 | |
| #processing streamflow data | |
| library("dataRetrieval") | |
| obsites <- c("07325840","07325850") | |
| runoff.obs<-readNWISdv(obsites, "00060",beginDate, endDate) | |
| runoff.obs$X_00060_00003 <- runoff.obs$X_00060_00003 * 0.3048^3 | |
| runoff.site<-readNWISsite(obsites) | |
| NOBS = nrow(runoff.obs) | |
| head(runoff.obs) | |
| tail(runoff.obs) | |
| nsites<-c(table(runoff.obs$site_no)) | |
| runoff.obs$hru <- rep(hrus, times = nsites) | |
| runoff.obs$weight <- 1.0 | |
| runoff.obs$weight[runoff.obs$hrus < max(hrus)] <- 1/length(hrus) | |
| runoff.obs$weight[runoff.obs$Date < as.Date("2006-10-1")] <- 0.0 | |
| runoffob <- data.frame(OBSNME = paste0("Q", sprintf("%02i", runoff.obs$hru),"D",format(runoff.obs$Date, "%y%m%d")), | |
| OBSVAL = runoff.obs$X_00060_00003, | |
| WEIGHT = runoff.obs$weight, | |
| OBGNME = paste0("flow", runoff.obs$hru)) | |
| head(runoffob) | |
| #convert to cms | |
| write.fwf(runoffob, "!!OBData.txt", quote = FALSE, rownames = FALSE, colnames = FALSE) | |
| obins <- "" | |
| l = 0 | |
| obins[(l=l+1)] <- "pif $" | |
| obins[(l=l+1)] <- "l9" | |
| minDate <- tapply(runoff.obs$Date, runoff.obs$hru, min) | |
| maxDate <- tapply(runoff.obs$Date, runoff.obs$hru, max) | |
| for (dd in 0:(ndays-1)){ | |
| for (i in 1:length(hrus)){ | |
| currentD <- as.Date(dd, beginDate) | |
| if (currentD >= minDate[i] & currentD <= maxDate[i]) | |
| obins[(l=l+1)] <- paste0("l",lad[i]," [Q", sprintf("%02i", hrus[i]),"D",format(currentD, "%y%m%d"), "]50:61") | |
| else | |
| obins[(l=l+1)] <- paste0("l",lad[i]) | |
| } | |
| } | |
| writeLines(obins, "!!OB.ins", sep = "\n") | |
| pst <- "" | |
| l = 0 | |
| # pcf | |
| pst[(l=l+1)] <- paste0("pcf") | |
| # * control data | |
| pst[(l=l+1)] <- paste0("* control data") | |
| # RSTFLE PESTMODE | |
| pst[(l=l+1)] <- paste0(" restart estimation") | |
| # NPAR NOBS NPARGP NPRIOR NOBSGP | |
| pst[(l=l+1)] <- paste(" ", NPAR, NOBS, NPARGP, NPRIOR, NOBSGP) | |
| # NTPLFLE NINSFLE PRECIS DPOINT NUMCOM JACFILE MESSFILE | |
| pst[(l=l+1)] <- paste0(" 1 1 single point 1 0 0") | |
| # RLAMBDA1 RLAMFAC PHIRATSUF PHIREDLAM NUMLAM | |
| pst[(l=l+1)] <- paste0(" 0.1 2.0 0.3 0.03 10") | |
| # RELPARMAX FACPARMAX FACORIG | |
| pst[(l=l+1)] <- paste0(" 5.0 5.0 0.001") | |
| # PHIREDSWH | |
| pst[(l=l+1)] <- paste0(" 0.1") | |
| # NOPTMAX PHIREDSTP NPHISTP NPHINORED RELPARSTP NRELPAR | |
| pst[(l=l+1)] <- paste0(" 30 .005 4 4 .005 4") | |
| # ICOV ICOR IEIG | |
| pst[(l=l+1)] <- paste0(" 1 1 1") | |
| # * singular value decomposition | |
| pst[(l=l+1)] <- paste0("* singular value decomposition") | |
| pst[(l=l+1)] <- paste0(" 1") | |
| pst[(l=l+1)] <- paste0(" 10 5.0000000E-07") | |
| pst[(l=l+1)] <- paste0(" 1") | |
| # * parameter groups | |
| pst[(l=l+1)] <- paste0("* parameter groups") | |
| # PARGPNME INCTYP DERINC DERINCLB FORCEN DERINCMUL DERMTHD | |
| # (one such line for each of the NPARGP parameter groups) | |
| pst[(l=l+1)] <- paste0(readLines("!!paragroup.dat"), collapse = "\n") | |
| # * parameter data | |
| pst[(l=l+1)] <- paste0("* parameter data") | |
| # PARNME PARTRANS PARCHGLIM PARVAL1 PARLBND PARUBND PARGP SCALE OFFSET DERCOM | |
| # (one such line for each of the NPAR parameters) | |
| pst[(l=l+1)] <- paste0(readLines("!!parameter.dat"), collapse = "\n") | |
| # PARNME PARTIED | |
| # (one such line for each tied parameter) | |
| # * observation groups | |
| pst[(l=l+1)] <- paste0("* observation groups") | |
| # OBGNME | |
| # (one such line for each observation group) | |
| pst[(l=l+1)] <- paste0(OBGNME, collapse = "\n") | |
| # * observation data | |
| pst[(l=l+1)] <- paste0("* observation data") | |
| # OBSNME OBSVAL WEIGHT OBGNME | |
| # (one such line for each of the NOBS observations) | |
| pst[(l=l+1)] <- paste0(readLines("!!OBData.txt"), collapse = "\n") | |
| # * model command line | |
| pst[(l=l+1)] <- paste0("* model command line") | |
| # write the command which PEST must use to run the model | |
| pst[(l=l+1)] <- paste0(" !!model_predict.bat") | |
| # * model input/output | |
| pst[(l=l+1)] <- paste0("* model input/output") | |
| # TEMPFLE INFLE | |
| # (one such line for each model input file containing parameters) | |
| pst[(l=l+1)] <- paste0(" !!par2par.tpl !!par2par.dat") | |
| # INSFLE OUTFLE | |
| # (one such line for each model output file containing observations) | |
| pst[(l=l+1)] <- paste0(" !!OB.ins output.rch") | |
| # * prior information | |
| pst[(l=l+1)] <- paste0("* prior information") | |
| # PILBL PIFAC * PARNME + PIFAC * log(PARNME) ... = PIVAL WEIGHT OBGNME | |
| # (one such line for each of the NPRIOR articles of prior information) | |
| writeLines(pst, "!!calib1.pst", sep = "\n") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # Produced by Michael Ou, 2/11/2015 | |
| # using PEST for SWAT calibration | |
| # grouping the same parameter for all HRUs | |
| # change of parameter is accommodated by multiplying a relative rate which is operated by PEST | |
| # for uniform parameters (same value for all HRUs), we replace it directly. | |
| # r__CN2.mgt -0.2 0.2 | |
| # v__ALPHA_BF.gw 0.0 1.0 | |
| # v__GW_DELAY.gw 30.0 450.0 | |
| # v__GWQMN.gw 0.0 2.0 | |
| # | |
| # | |
| # v__GW_REVAP.gw 0.0 0.2 | |
| # v__ESCO.hru 0.8 1.0 | |
| # v__CH_N2.rte 0.0 0.3 | |
| # v__CH_K2.rte 5.0 130.0 | |
| # v__ALPHA_BNK.rte 0.0 1.0 | |
| # r__SOL_AWC(1).sol -0.2 0.4 | |
| # r__SOL_K(1).sol -0.8 0.8 | |
| # r__SOL_BD(1).sol -0.5 0.6 | |
| # v__SFTMP.bsn -5.0 5.0 | |
| # | |
| # | |
| # | |
| # -------------------------------------------------- | |
| # example of parameterization: | |
| # | |
| # r__Precipitation(){1990001-2000265}.pcp 0 0.4 | |
| # r__AUTO_EFF{[],11,AUTO_NSTRS=0.6}.mgt 0.2 0.6 | |
| # r__BLAI{22}.crop.dat 2 8 | |
| # | |
| # | |
| # | |
| # ----------------------------------------- | |
| # | |
| # r__CN2.mgt________1 -0.2 0.2 | |
| # r__SOL_AWC(1).sol________1 -0.2 0.1 | |
| # r__SOL_K(1).sol________1 -0.8 0.8 | |
| # r__SOL_BD(1).sol________1 -0.5 0.6 | |
| # a__GWQMN.gw________1 0.0 25.0 | |
| # a__GW_REVAP.gw________1 -0.1 0.0 | |
| # v__REVAPMN.gw________1 0.0 10.0 | |
| # a__ESCO.hru________1 0.0 0.2 | |
| # r__HRU_SLP.hru________1 0.0 0.2 | |
| # r__OV_N.hru________1 -0.2 0.0 | |
| # r__SLSUBBSN.hru________1 0.0 0.2 | |
| # | |
| # | |
| # | |
| # r__CN2.mgt________3-6 -0.2 0.2 | |
| # r__SOL_AWC(1).sol________3-6 -0.2 0.1 | |
| # r__SOL_K(1).sol________3-6 -0.8 0.8 | |
| # r__SOL_BD(1).sol________3-6 -0.5 0.6 | |
| # a__GWQMN.gw________3-6 0.0 25.0 | |
| # a__GW_REVAP.gw________3-6 -0.1 0.0 | |
| # v__REVAPMN.gw________3-6 0.0 10.0 | |
| # a__ESCO.hru________3-6 0.0 0.2 | |
| # r__HRU_SLP.hru________3-6 0.0 0.2 | |
| # r__OV_N.hru________3-6 -0.2 0.0 | |
| # r__SLSUBBSN.hru________3-6 0.0 0.2 | |
| # | |
| # | |
| # r__CN2.mgt________7-20 -0.2 0.2 | |
| # r__SOL_AWC(1).sol________7-20 -0.2 0.1 | |
| # r__SOL_K(1).sol________7-20 -0.8 0.8 | |
| # r__SOL_BD(1).sol________7-20 -0.5 0.6 | |
| # a__GWQMN.gw________7-20 0.0 25.0 | |
| # a__GW_REVAP.gw________7-20 -0.1 0.0 | |
| # v__REVAPMN.gw________7-20 0.0 10.0 | |
| # a__ESCO.hru________7-20 0.0 0.2 | |
| # r__HRU_SLP.hru________7-20 0.0 0.2 | |
| # r__OV_N.hru________7-20 -0.2 0.0 | |
| # r__SLSUBBSN.hru________7-20 0.0 0.2 | |
| # | |
| # | |
| # v__ALPHA_BF.gw 0.0 1.0 | |
| # v__GW_DELAY.gw 30.0 450.0 | |
| # v__CH_N2.rte 0.0 0.3 | |
| # v__CH_K2.rte 5.0 130.0 | |
| # v__ALPHA_BNK.rte 0.0 1.0 | |
| # v__SFTMP.bsn -5.0 5.0 | |
| library("stringr") | |
| library("gdata") | |
| setwd("D:\\@Workspace\\Okalahoma_Lake_Cr\\CN_46_HRU\\PESTUN") | |
| NPAR = 0 | |
| NOBS = 0 | |
| NPARGP = 0 | |
| NPRIOR = 0 | |
| NOBSGP = 0 | |
| # PARAMETER GROUP | |
| # PARGPNME INCTYP DERINC DERINCLB FORCEN DERINCMUL DERMTHD | |
| #123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc | |
| # 1 2 3 4 5 6 7 8 9 10 | |
| PARGPNME <- c( "SLSUBBSN", "HRU_SLP", "OV_N", "ESCO", "CN2", "BD", "AWC", "KSAT", "GW_DELAY", "ALPHA_BF", | |
| "GWQMN", "RCHRG_DP", "GW_SPYLD", "KS", "QR", "QS", "LAMBDA", "HBUB", "WCET") | |
| NPARGP = length(PARGPNME) | |
| #INCTYP is a character variable which can assume the values "relative", "absolute" or "rel_to_max". | |
| INCTYP<-rep("relative", NPARGP) | |
| #DERINC is the increment used for forward-difference calculation of derivatives | |
| DERINC <- c( 0.1, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.5, 0.1, | |
| 0.1, 0.1, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05) | |
| #DERINCLB is an absolute lower bound can be placed on parameter increments | |
| DERINCLB <-rep(0.01, NPARGP) | |
| #The character variable FORCEN (an abbreviation of "FORward/CENtral") determines whether derivatives for group members are calculated using | |
| # forward differences, one of the variants of the central difference method, of whether both alternatives are used in the course of an optimisation run. | |
| # It must assume one of the values "always_2", "always_3" or "switch". | |
| # Experience has shown that in most instances the most appropriate value for FORCEN is "switch". | |
| FORCEN = rep("switch", NPARGP) | |
| #DERINCMUL Whenever the central method is employed for derivatives calculation, DERINC is multiplied by DERINCMUL | |
| # Experience shows that a value between 1.0 and 2.0 is usually satisfactory. | |
| DERINCMUL = rep(2.0, NPARGP) | |
| # | |
| DERMTHD = rep("parabolic", NPARGP) | |
| ParaGroup <- data.frame(PARGPNME,INCTYP,DERINC,DERINCLB,FORCEN,DERINCMUL,DERMTHD) | |
| row.names(ParaGroup) <- PARGPNME | |
| # Parameter Data | |
| # PARNME PARTRANS PARCHGLIM PARVAL1 PARLBND PARUBND PARGP SCALE OFFSET DERCOM | |
| #123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc | |
| # 1 2 3 4 5 6 7 8 9 10 | |
| PARNME <- PARGPNME | |
| NPAR = length(PARNME) | |
| #PARTRANS is a character variable which must assume one of four values, viz. "none", "log", "fixed" or "tied". | |
| PARTRANS<-rep("none", NPAR) | |
| #This character variable is used to designate whether an adjustable parameter is relative-limited or factor-limited; | |
| # PARCHGLIM must be provided with one of two possible values, viz. "relative" or "factor". | |
| PARCHGLIM<-rep("factor", NPAR) | |
| #PARVAL1, a real variable, is a parameter's initial value. | |
| PARVAL1 <- rep(1, NPAR) | |
| #PARLBND and PARUBND: These two real variables represent a parameter's lower and upper bounds respectively. | |
| PARLBND <- c(0.01, 0.05, 0.1, 0.05, 0.8, 0.1, 0.3, 0.001, 0.1, 0.1, | |
| 0.01, 0.1, 0.05, 0.0001, 0.1, 0.1, 0.1, 0.01, 0.1) | |
| PARUBND <-c( 10.0, 10.0, 10.0, 10.0, 1.3, 3.0, 3.0, 1.0e3, 20, 10, | |
| 1e5, 10.0, 20.0, 1e4, 10, 10, 10, 100, 3) | |
| #PARGP is the name of the group to which a parameter belongs. | |
| PARGP = PARGPNME | |
| #SCALE and OFFSET | |
| #Just before a parameter value is written to a model input file (be it for initial determination of the objective function, | |
| # derivatives calculation or parameter upgrade), it is multiplied by the real variable SCALE, after which the real variable OFFSET is added. | |
| SCALE = rep(1.0, NPAR) | |
| OFFSET = rep(0.0, NPAR) | |
| #DERCOM: Unless using PEST's external derivatives functionality (see Chapter 9), this variable should be set to 1. | |
| DERCOM = rep(1, NPAR) | |
| # PARNME PARTRANS PARCHGLIM PARVAL1 PARLBND PARUBND PARGP SCALE OFFSET DERCOM | |
| ParaData <- data.frame(PARNME, PARTRANS, PARCHGLIM, PARVAL1, PARLBND, PARUBND, PARGP, SCALE, OFFSET, DERCOM) | |
| row.names(ParaData) <- PARNME | |
| # real values | |
| #123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc123456789abc | |
| # 1 2 3 4 5 6 7 8 9 10 | |
| # "SLSUBBSN", "HRU_SLP", "OV_N", "ESCO", "CN2", "BD", "AWC", "KSAT", "GW_DELAY", "ALPHA_BF", | |
| # "GWQMN", "RCHRG_DP", "GW_SPYLD", "KS", "QR", "QS", "LAMBDA", "HBUB", "WCET" | |
| REALLBND <-c(10.0, 0.001, 0.008, 0.8, 35, 0.9, 0.05, 1e-4, 1, 0.01, | |
| 1, 0.01, 0.01, 1e-4, 1e-4, 0.1, 0.04, 0.1, 0.01) | |
| REALUBND <-c(1000, 0.3, 0.6, 1.0, 98, 2.0, 0.5, 1e4, 1000, 0.99, | |
| 1e5, 0.99, 0.5, 5e3, 0.2, 0.6, 1.2, 5000, 0.3) | |
| RealPara <- data.frame(PARNME, REALLBND, REALUBND) | |
| row.names(RealPara) <- PARNME | |
| #read number of sub | |
| nsub = 0 | |
| pname <- "" | |
| pgroup <- "" | |
| pvalue <- 0.0 | |
| pline <- "" | |
| tplname <- "" #this is the template file for par2par but not the pest | |
| ntpl <- 0 | |
| np <- 0 | |
| filename<-paste0("fig.fig") | |
| fig <- paste(readLines(filename), sep = "\n") | |
| fig <- chartr("\t"," ", fig) #delete tab | |
| ZBOT = 0 | |
| BD = 0 | |
| AWC = 0 | |
| KSAT = 0 | |
| isub = 0 | |
| ihru = 0 | |
| i = 0 | |
| hrutot = 0 | |
| for (isub in 1:length(fig)){ | |
| sline <- unlist(strsplit(fig[isub]," ")) | |
| sline <- sline[sline != ""] | |
| if (sline[1] == "subbasin") { | |
| nsub = nsub + 1 | |
| sline <- unlist(strsplit(fig[isub+1]," ")) | |
| sline <- sline[sline != ""] | |
| #### sub file ### | |
| filename<-sline[1] | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| sub <- paste(readLines(filename), sep = "\n") | |
| #CH_N1 | |
| #np=np+1 | |
| #pvalue[np]<-as.numeric(substr(sub[29], 1, 16)) | |
| #pname[np]<- paste0("CH_N1",isub) | |
| #pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.008 ","\t0.6","\tCH_N1 "," \t1.0 \t0.0 \t1") | |
| #substr(sub[29], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", sub), tplname[ntpl], sep = "\n") | |
| nhru = as.numeric(substr(sub[53], 1, 16)) | |
| hrutot = hrutot + nhru | |
| iisub = nsub | |
| for (ihru in 1:nhru){ | |
| ### hru file ### | |
| filename<-substr(sub[61+ihru], 1, 13) | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| hru <- paste(readLines(filename), sep = "\n") | |
| #SLSUBBSN | |
| np=np+1 | |
| pgroup[np] <- "SLSUBBSN" | |
| pvalue[np]<-as.numeric(substr(hru[3], 1, 16)) | |
| pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t10.0 ","\t1000.0","\t",pgroup[np]," \t1.0 \t0.0 \t1") | |
| substr(hru[3], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #HRU_SLP | |
| np=np+1 | |
| pgroup[np] <- "HRU_SLP" | |
| pvalue[np]<-as.numeric(substr(hru[4], 1, 16)) | |
| pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.005 ","\t0.2","\t",pgroup[np]," \t1.0 \t0.0 \t1") | |
| substr(hru[4], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #OV_N | |
| np=np+1 | |
| pgroup[np] <- "OV_N" | |
| pvalue[np]<-as.numeric(substr(hru[5], 1, 16)) | |
| pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.008 ","\t0.6","\t",pgroup[np]," \t1.0 \t0.0 \t1") | |
| substr(hru[5], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #ESCO | |
| np=np+1 | |
| pgroup[np] <- "ESCO" | |
| pvalue[np]<-as.numeric(substr(hru[11], 1, 16)) | |
| pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.8 ","\t1.0 ","\tESCO "," \t1.0 \t0.0 \t1") | |
| substr(hru[11], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| # #SURLAG | |
| # np=np+1 | |
| # pvalue[np]<-as.numeric(substr(hru[11], 1, 16)) | |
| # pname[np]<- paste0("SURLAG",iisub,ihru) | |
| # pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0 ","\t0.6","\tSURLAG "," \t1.0 \t0.0 \t1") | |
| # substr(hru[5], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", hru), tplname[ntpl], sep = "\n") | |
| ### mgt ### | |
| filename<-substr(sub[61+ihru], 14, 26) | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| mgt <- paste(readLines(filename), sep = "\n") | |
| #CN2 | |
| np=np+1 | |
| pgroup[np] <- "CN2" | |
| pvalue[np]<-as.numeric(substr(mgt[11], 1, 16)) | |
| pname[np]<- paste0("CN2",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t35 ","\t98","\tCN2 "," \t1.0 \t0.0 \t1") | |
| substr(mgt[11], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", mgt), tplname[ntpl], sep = "\n") | |
| ### sol ### | |
| # filename<-substr(sub[61+ihru], 27, 39) | |
| # filename <-gsub("^\\s+|\\s+$", "", filename) | |
| # sol <- paste(readLines(filename), sep = "\n") | |
| # | |
| # NLAY <- (str_length(sol[8]) - 27) / 12 | |
| # INDEX <- 28 + ((1:(NLAY+1))-1) * 12 | |
| # ii <- 1:NLAY | |
| # for (i in ii){ | |
| # #BD | |
| # np=np+1 | |
| # pgroup[np] <- "BD" | |
| # pvalue[np]<-as.numeric(substr(sol[9], INDEX[i], INDEX[i]+11)) | |
| # pname[np]<- paste0("BD",iisub,ihru,i) | |
| # pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.9 ","\t2.0","\tBD "," \t1.0 \t0.0 \t1") | |
| # substr(sol[9], INDEX[i], INDEX[i]+11) <- paste0("$",substr(paste0(pname[np]," "),1,10),"$") | |
| # | |
| # | |
| # #AWC | |
| # np=np+1 | |
| # pgroup[np] <- "AWC" | |
| # pvalue[np]<-as.numeric(substr(sol[10], INDEX[i], INDEX[i]+11)) | |
| # pname[np]<- paste0("AWC",iisub,ihru,i) | |
| # pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.05 ","\t0.8","\tAWC "," \t1.0 \t0.0 \t1") | |
| # substr(sol[10], INDEX[i], INDEX[i]+11) <- paste0("$",substr(paste0(pname[np]," "),1,10),"$") | |
| # | |
| # | |
| # #KSAT | |
| # np=np+1 | |
| # pgroup[np] <- "KSAT" | |
| # pvalue[np]<- as.numeric(substr(sol[11], INDEX[i], INDEX[i]+11)) | |
| # pname[np]<- paste0("KSAT",iisub,ihru,i) | |
| # pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t",pvalue[np]," \t0.0001 ","\t2000","\tKSAT "," \t1.0 \t0.0 \t1") | |
| # substr(sol[11], INDEX[i], INDEX[i]+11) <- paste0("$",substr(paste0(pname[np]," "),1,10),"$") | |
| # } | |
| # | |
| # | |
| # ntpl = ntpl + 1 | |
| # tplname[ntpl] <- paste0(filename,".tpl") | |
| # writeLines(c("ptf $", sol), tplname[ntpl], sep = "\n") | |
| ### usf ### | |
| filename<-paste0(substr(sub[61+ihru], 27, 35),".usf") | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| usf <- paste(readLines(filename), sep = "\n") | |
| #skip the comment lines | |
| for (ncmt in 1:length(usf)){ | |
| if (substr(usf[ncmt],1,1) != "#") break() | |
| } | |
| #WCET | |
| np=np+1 | |
| pgroup[np] <- "WCET" | |
| pvalue[np]<-as.numeric(substr(usf[5+ncmt], 1, 16)) | |
| pname[np]<- paste0(pgroup[np],iisub,ihru) | |
| substr(usf[5+ncmt], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", usf), tplname[ntpl], sep = "\n") | |
| ### chm ### | |
| ### gw ### | |
| filename<-substr(sub[61+ihru], 53, 65) | |
| filename <-gsub("^\\s+|\\s+$", "", filename) | |
| gw <- paste(readLines(filename), sep = "\n") | |
| #GW_DELAY | |
| np=np+1 | |
| pgroup[np] <- "GW_DELAY" | |
| pvalue[np]<-as.numeric(substr(gw[4], 1, 16)) | |
| pname[np]<- paste0("GW_DELAY",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t20 ","\t450","\tGW_DELAY"," \t1.0 \t0.0 \t1") | |
| substr(gw[4], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #ALPHA_BF | |
| np=np+1 | |
| pgroup[np] <- "ALPHA_BF" | |
| pvalue[np]<-as.numeric(substr(gw[5], 1, 16)) | |
| pname[np]<- paste0("ALPHA_BF",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t0.001 ","\t0.999","\tALPHA_BF"," \t1.0 \t0.0 \t1") | |
| substr(gw[5], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #GWQMN | |
| np=np+1 | |
| pgroup[np] <- "GWQMN" | |
| pvalue[np]<-as.numeric(substr(gw[6], 1, 16)) | |
| if (pvalue[np]<0.1) pvalue[np] <- 1.0 | |
| pname[np]<- paste0("GWQMN",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t0.01 ","\t10000.0","\tGWQMN"," \t1.0 \t0.0 \t1") | |
| substr(gw[6], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #RCHRG_DP | |
| np=np+1 | |
| pgroup[np] <- "RCHRG_DP" | |
| pvalue[np]<-as.numeric(substr(gw[9], 1, 16)) | |
| pname[np]<- paste0("RCHRG_DP",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t0.001 ","\t0.9","\tRCHRG_DP","\t1.0 \t0.0 \t1") | |
| substr(gw[9], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| #GW_SPYLD | |
| np=np+1 | |
| pgroup[np] <- "GW_SPYLD" | |
| pvalue[np]<-as.numeric(substr(gw[11], 1, 16)) | |
| pname[np]<- paste0("GW_SPYLD",iisub,ihru) | |
| pline[np]<- paste0(pname[np],"\tlog ","\tfactor ","\t", pvalue[np]," \t0.01 ","\t0.6","\tGW_SPYLD","\t1.0 \t0.0 \t1") | |
| substr(gw[11], 1, 16) <- paste0("$",substr(paste0(pname[np]," "),1,14),"$") | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", gw), tplname[ntpl], sep = "\n") | |
| } | |
| } | |
| } | |
| #translate soil mat file | |
| filename <- "soil.mat" | |
| mat <- paste(readLines(filename), sep = "\n") | |
| mat <- chartr("\t"," ", mat) | |
| for (ncmt in 1:length(mat)){ | |
| if (substr(mat[ncmt],1,1) != "#") break() | |
| } | |
| sline <- unlist(strsplit(mat[ncmt]," ")) | |
| sline <- sline[sline != ""] | |
| nmat <- as.numeric(sline[1]) | |
| for (i in 1:(nmat)){ | |
| sline <- unlist(strsplit(mat[i+ncmt]," ")) | |
| sline <- sline[sline != ""] | |
| #KS | |
| np=np+1 | |
| pgroup[np] <- "KS" | |
| pvalue[np]<-as.numeric(sline[1]) | |
| pname[np]<- paste0(pgroup[np],i) | |
| #QR | |
| np=np+1 | |
| pgroup[np] <- "QR" | |
| pvalue[np]<-as.numeric(sline[2]) | |
| pname[np]<- paste0(pgroup[np],i) | |
| #QS | |
| np=np+1 | |
| pgroup[np] <- "QS" | |
| pvalue[np]<-as.numeric(sline[3]) | |
| pname[np]<- paste0(pgroup[np],i) | |
| #LAMBDA | |
| np=np+1 | |
| pgroup[np] <- "LAMBDA" | |
| pvalue[np]<-as.numeric(sline[4]) | |
| pname[np]<- paste0(pgroup[np],i) | |
| #ETA calculated by par2par | |
| #np=np+1 | |
| #pvalue[np]<-as.numeric(sline[5]) | |
| #pname[np]<- paste0("LAMBDA",i) | |
| #pline[np]<- paste0(pname[np],"\ttied "," \tfactor ","\t",pvalue[np]," \t10.0 ","\t1000.0","\tETA","\t1.0 \t0.0 \t1") | |
| #HBUB | |
| np=np+1 | |
| pgroup[np] <- "HBUB" | |
| pvalue[np]<- -as.numeric(sline[6]) | |
| pname[np]<- paste0(pgroup[np],i) | |
| mat[i+1] <- paste0(paste0("$",substr(paste0(pname[np-4]," "),1,14),"$\t "), | |
| paste0("$",substr(paste0(pname[np-3]," "),1,14),"$\t "), | |
| paste0("$",substr(paste0(pname[np-2]," "),1,14),"$\t "), | |
| paste0("$",substr(paste0(pname[np-1]," "),1,14),"$\t "), | |
| paste0("$",substr(paste0("ETA",i," "),1,14),"$\t "), | |
| paste0("$",substr(paste0(pname[np]," "),1,14),"$")) | |
| } | |
| ntpl = ntpl + 1 | |
| tplname[ntpl] <- paste0(filename,".tpl") | |
| writeLines(c("ptf $", mat), tplname[ntpl], sep = "\n") | |
| # paragroup | |
| groups <- unique(pgroup) | |
| ngroup = length(groups) | |
| NPAR = ngroup | |
| NPARGP = ngroup | |
| ParaGroup <- subset(ParaGroup, PARGPNME %in% groups) | |
| #write.table(ParaGroup, "!!paragroup.dat", quote = FALSE, sep = " \t", row.names = FALSE, col.names = FALSE) | |
| write.fwf(ParaGroup, "!!paragroup.dat", quote = FALSE, rownames = FALSE, colnames = FALSE) | |
| # * parameter data for PEST (parameter groups) | |
| ParaData <- subset(ParaData, PARNME %in% groups) | |
| write.fwf(ParaData, "!!parameter.dat", quote = FALSE, rownames = FALSE, colnames = FALSE) | |
| #write.table(ParaData, "!!parameter.dat", quote = FALSE, row.names = FALSE, col.names = FALSE) | |
| ppestline = "" | |
| #for par2par tpl | |
| p2pline <- "" | |
| l=0 | |
| p2pline[(l=l+1)] <- "ptf $" | |
| p2pline[(l=l+1)] <- "* parameter data" | |
| for (gg in groups){ | |
| p2pline[(l=l+1)] <- paste(substr(paste0(gg," "),1,20), "=", paste0("$",substr(paste0(gg," "),1,14),"$")) | |
| } | |
| for (i in 1:np) { | |
| p2pline[(l=l+1)] <- paste(substr(paste0(pname[i]," "),1,20), "=", "max ( min (", | |
| substr(paste0(pgroup[i]," "),1,10), | |
| "*", substr(paste0(pvalue[i]," "),1,10), ",", | |
| substr(paste0(RealPara[pgroup[i],"REALUBND"]," "),1,10), "),", RealPara[pgroup[i],"REALLBND"], ")") | |
| if (pgroup[i] == "QS"){ | |
| p2pline[l] <- paste(substr(paste0(pname[i]," "),1,20), "=", "max ( min (", | |
| substr(paste0(pgroup[i]," "),1,10), | |
| "*", substr(paste0(pvalue[i]," "),1,10), ",", | |
| substr(paste0(RealPara[pgroup[i],"REALUBND"]," "),1,10), "),", pname[i-1], " + 0.1 )") | |
| } | |
| if (pgroup[i] == "LAMBDA"){ | |
| p2pline[(l=l+1)] <- paste(sprintf("%-20s", gsub("LAMBDA", "ETA", pname[i])), "= 2 /", pname[i],"+ 3") | |
| } | |
| if (pgroup[i] == "HBUB"){ | |
| p2pline[l] <- paste(substr(paste0(pname[i]," "),1,20), "=", "- max ( min (", | |
| substr(paste0(pgroup[i]," "),1,10), | |
| "*", substr(paste0(pvalue[i]," "),1,10), ",", | |
| substr(paste0(RealPara[pgroup[i],"REALUBND"]," "),1,10), "),", RealPara[pgroup[i],"REALLBND"], ")") | |
| } | |
| } | |
| p2pline[(l=l+1)] <- "* template and model input files" | |
| for (i in 1:ntpl){ | |
| p2pline[(l=l+1)] <- paste(tplname[i], substr(tplname[i], 1, str_length(tplname[i])-4)) | |
| } | |
| p2pline[(l=l+1)] <- "* control data" | |
| p2pline[(l=l+1)] <- "single point" | |
| writeLines(p2pline, "!!par2par.tpl", sep = "\n") | |
| #create observation data and instruction file | |
| nhru = hrutot | |
| hrus <- c(17,hrutot) | |
| NOBSGP = length(hrus) | |
| OBGNME <- paste0("flow", hrus) | |
| lad <- c(hrus[1],diff(hrus)) | |
| cio <- paste(readLines("file.cio"), sep = "\n") | |
| NBYR <- as.numeric(substr(cio[8], 1, 16)) | |
| IYR <- as.numeric(substr(cio[9], 1, 16)) | |
| IDAF <- as.numeric(substr(cio[10], 1, 16)) | |
| IDAL <- as.numeric(substr(cio[11], 1, 16)) | |
| NYSKIP <- as.numeric(substr(cio[60], 1, 16)) | |
| beginDate <- as.Date(IDAF-1, as.Date(paste0(IYR+NYSKIP,"-1-1"))) | |
| endDate <- as.Date(IDAL-1, as.Date(paste0(IYR+NBYR-1,"-1-1"))) | |
| ndays <- endDate - beginDate + 1 | |
| #processing streamflow data | |
| library("dataRetrieval") | |
| obsites <- c("07325840","07325850") | |
| runoff.obs<-readNWISdv(obsites, "00060",beginDate, endDate) | |
| runoff.obs$X_00060_00003 <- runoff.obs$X_00060_00003 * 0.3048^3 | |
| runoff.site<-readNWISsite(obsites) | |
| NOBS = nrow(runoff.obs) | |
| head(runoff.obs) | |
| tail(runoff.obs) | |
| nsites<-c(table(runoff.obs$site_no)) | |
| runoff.obs$hru <- rep(hrus, times = nsites) | |
| runoff.obs$weight <- 1.0 | |
| runoff.obs$weight[runoff.obs$Date < as.Date("2006-1-1")] <- 0.0 | |
| runoffob <- data.frame(OBSNME = paste0("Q", sprintf("%02i", runoff.obs$hru),"D",format(runoff.obs$Date, "%y%m%d")), | |
| OBSVAL = runoff.obs$X_00060_00003, | |
| WEIGHT = runoff.obs$weight, | |
| OBGNME = paste0("flow", runoff.obs$hru)) | |
| head(runoffob) | |
| #convert to cms | |
| write.fwf(runoffob, "!!OBData.txt", quote = FALSE, rownames = FALSE, colnames = FALSE) | |
| obins <- "" | |
| l = 0 | |
| obins[(l=l+1)] <- "pif $" | |
| obins[(l=l+1)] <- "l9" | |
| minDate <- tapply(runoff.obs$Date, runoff.obs$hru, min) | |
| maxDate <- tapply(runoff.obs$Date, runoff.obs$hru, max) | |
| for (dd in 0:(ndays-1)){ | |
| for (i in 1:length(hrus)){ | |
| currentD <- as.Date(dd, beginDate) | |
| if (currentD >= minDate[i] & currentD <= maxDate[i]) | |
| obins[(l=l+1)] <- paste0("l",lad[i]," [Q", sprintf("%02i", hrus[i]),"D",format(currentD, "%y%m%d"), "]50:61") | |
| else | |
| obins[(l=l+1)] <- paste0("l",lad[i]) | |
| } | |
| } | |
| writeLines(obins, "!!OB.ins", sep = "\n") | |
| pst <- "" | |
| l = 0 | |
| # pcf | |
| pst[(l=l+1)] <- paste0("pcf") | |
| # * control data | |
| pst[(l=l+1)] <- paste0("* control data") | |
| # RSTFLE PESTMODE | |
| pst[(l=l+1)] <- paste0(" restart estimation") | |
| # NPAR NOBS NPARGP NPRIOR NOBSGP | |
| pst[(l=l+1)] <- paste(" ", NPAR, NOBS, NPARGP, NPRIOR, NOBSGP) | |
| # NTPLFLE NINSFLE PRECIS DPOINT NUMCOM JACFILE MESSFILE | |
| pst[(l=l+1)] <- paste0(" 1 1 single point 1 0 0") | |
| # RLAMBDA1 RLAMFAC PHIRATSUF PHIREDLAM NUMLAM | |
| pst[(l=l+1)] <- paste0(" 0.1 2.0 0.3 0.03 10") | |
| # RELPARMAX FACPARMAX FACORIG | |
| pst[(l=l+1)] <- paste0(" 5.0 5.0 0.001") | |
| # PHIREDSWH | |
| pst[(l=l+1)] <- paste0(" 0.1") | |
| # NOPTMAX PHIREDSTP NPHISTP NPHINORED RELPARSTP NRELPAR | |
| pst[(l=l+1)] <- paste0(" 30 .005 4 4 .005 4") | |
| # ICOV ICOR IEIG | |
| pst[(l=l+1)] <- paste0(" 1 1 1") | |
| # * singular value decomposition | |
| pst[(l=l+1)] <- paste0("* singular value decomposition") | |
| pst[(l=l+1)] <- paste0(" 1") | |
| pst[(l=l+1)] <- paste0(" 10 5.0000000E-07") | |
| pst[(l=l+1)] <- paste0(" 1") | |
| # * parameter groups | |
| pst[(l=l+1)] <- paste0("* parameter groups") | |
| # PARGPNME INCTYP DERINC DERINCLB FORCEN DERINCMUL DERMTHD | |
| # (one such line for each of the NPARGP parameter groups) | |
| pst[(l=l+1)] <- paste0(readLines("!!paragroup.dat"), collapse = "\n") | |
| # * parameter data | |
| pst[(l=l+1)] <- paste0("* parameter data") | |
| # PARNME PARTRANS PARCHGLIM PARVAL1 PARLBND PARUBND PARGP SCALE OFFSET DERCOM | |
| # (one such line for each of the NPAR parameters) | |
| pst[(l=l+1)] <- paste0(readLines("!!parameter.dat"), collapse = "\n") | |
| # PARNME PARTIED | |
| # (one such line for each tied parameter) | |
| # * observation groups | |
| pst[(l=l+1)] <- paste0("* observation groups") | |
| # OBGNME | |
| # (one such line for each observation group) | |
| pst[(l=l+1)] <- paste0(OBGNME, collapse = "\n") | |
| # * observation data | |
| pst[(l=l+1)] <- paste0("* observation data") | |
| # OBSNME OBSVAL WEIGHT OBGNME | |
| # (one such line for each of the NOBS observations) | |
| pst[(l=l+1)] <- paste0(readLines("!!OBData.txt"), collapse = "\n") | |
| # * model command line | |
| pst[(l=l+1)] <- paste0("* model command line") | |
| # write the command which PEST must use to run the model | |
| pst[(l=l+1)] <- paste0(" !!model_predict.bat") | |
| # * model input/output | |
| pst[(l=l+1)] <- paste0("* model input/output") | |
| # TEMPFLE INFLE | |
| # (one such line for each model input file containing parameters) | |
| pst[(l=l+1)] <- paste0(" !!par2par.tpl !!par2par.dat") | |
| # INSFLE OUTFLE | |
| # (one such line for each model output file containing observations) | |
| pst[(l=l+1)] <- paste0(" !!OB.ins output.rch") | |
| # * prior information | |
| pst[(l=l+1)] <- paste0("* prior information") | |
| # PILBL PIFAC * PARNME + PIFAC * log(PARNME) ... = PIVAL WEIGHT OBGNME | |
| # (one such line for each of the NPRIOR articles of prior information) | |
| writeLines(pst, "!!calib1.pst", sep = "\n") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment