Nothing
## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## -----------------------------------------------------------------------------
library(HBV.IANIGLA)
data("semi_distributed_hbv")
str(semi_distributed_hbv)
## ----definition, echo = TRUE--------------------------------------------------
## brief arguments description
# basin: data frame with the same structure of the data("semi_distributed_hbv) (colnames included).
# tair: numeric matrix with air temperature inputs.
# precip: numeric matrix with precipitation inputs.
# pet: numeric matrix with potential eavapotranspiration inputs.
# param_snow: numeric vector with snow module parameters.
# param_soil: numeric vector with soil moisture parameters.
# param_routing: numeric vector with the routing parameters.
# param_tf: numeric vector with the transfer function parameter.
# init_snow: numeric value with initial snow water equivalent. Default value being 20 mm.
# init_soil: numeric value with initial soil moisture content. Default value being 100 mm.
# init_routing: numeric vector with bucket water initial values. Default values are 0 mm.
## output
# simulated streamflow series.
hydrological_hbv <- function(basin,
tair,
precip,
pet,
param_snow,
param_soil,
param_route,
param_tf,
init_snow = 20,
init_soil = 0,
init_routing = c(0, 0, 0)
){
n_it <- nrow(basin)
# create output lists
snow_module <- list()
soil_module <- list()
route_module <- list()
tf_module <- list()
# snow and soil module in every elevation band
for(i in 1:n_it){
snow_module[[ i ]] <-
SnowGlacier_HBV(model = 1, inputData = cbind(tair[ , i], precip[ , i]),
initCond = c(init_snow, 2), param = param_snow)
soil_module[[ i ]] <-
Soil_HBV(model = 1, inputData = cbind(snow_module[[i]][ , 5] , pet[ , i]),
initCond = c(init_soil, basin[i, 'rel_area']), param = param_soil )
} # end for
# get total soil discharge
soil_disch <- lapply(X = 1:n_it, FUN = function(x){
out <- soil_module[[x]][ , 1]
})
soil_disch <- Reduce(f = `+`, x = soil_disch)
# route module
route_module <- Routing_HBV(model = 1, lake = F, inputData = as.matrix(soil_disch),
initCond = init_routing, param = param_route )
# transfer function
tf_module <- round(
UH(model = 1, Qg = route_module[ , 1], param = param_tf), 4
)
return(tf_module)
}# end fun
## -----------------------------------------------------------------------------
# first we are going to create set the parameter range
# snow module
snow_range <- rbind(
sfcf = c(0, 1.5),
tr = c(-1, 1),
tt = c(0, 3),
fm = c(1.5, 4)
)
# soil module
soil_range <- rbind(
fc = c(100, 200),
lp = c(0.5, 1),
beta = c(1, 3)
)
# routing module (here I will give you the correct values)
routing_range <- rbind(
k0 = c(0.09, 0.09),
k1 = c(0.07, 0.07),
k2 = c(0.05, 0.05),
uzl = c(5, 5),
perc = c(2, 2)
)
# transfer function module (I will give the correct value)
tf_range <- rbind(
bmax = c(2.25, 2.25)
)
## -----------------------------------------------------------------------------
param_range <-
rbind(
snow_range,
soil_range,
routing_range,
tf_range
)
head(param_range)
## -----------------------------------------------------------------------------
# set the number of model runs that you want to try
n_run <- 1000
# build the matrix
n_it <- nrow(param_range)
param_sets <- matrix(NA_real_, nrow = n_run, ncol = n_it)
colnames(param_sets) <- rownames(param_range)
for(i in 1:n_it){
param_sets[ , i] <- runif(n = n_run,
min = param_range[i, 1],
max = param_range[i, 2]
)
}
head(param_sets)
## -----------------------------------------------------------------------------
# goodness of fit vector
gof <- c()
# make a loop
for(i in 1:n_run){
streamflow <- hydrological_hbv(
basin = semi_distributed_hbv$basin,
tair = semi_distributed_hbv$tair,
precip = semi_distributed_hbv$prec,
pet = semi_distributed_hbv$pet,
param_snow = param_sets[i, rownames(snow_range) ],
param_soil = param_sets[i, rownames(soil_range)],
param_route = param_sets[i, rownames(routing_range)],
param_tf = param_sets[i, rownames(tf_range)]
)
gof[i] <- cor(x = streamflow, y = semi_distributed_hbv$qout)
}
param_sets <- cbind(param_sets, gof)
head(param_sets)
## -----------------------------------------------------------------------------
# get the row index
max_gof <- which.max(param_sets[ , "gof"])
# extract the parameter set
param_opt <- param_sets[max_gof, ]
param_opt
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.