In EURAC-Ecohydro/geotopOptim2: An R Plug-in to calibrate soil moisture output of the Distributed Hydrological Model GEOtop
Header row {data-height=15}
**Poster title**
**Author**
Abstract
Eco-hydrological models are increasingly used in the contexts of hydrology, ecology, precision agriculture for a detailed description of the water cycle at various scales: local scale, an hillslope or a watershed. However, with increasing computing power and observations available, bigger and bigger amount of raw data are produced. Therefore the need to develop flexible and user-oriented interfaces to visualize multiple outputs, perform sensitivity analyzes and compare against observations emerges.
This work presents two R open-source packages: geotopbricks and geotopOptim2. They offer an I/0 interface and R visualization tools the GEOtop hydrological distributed model (http://geotopmodel.github.io/geotop/), which solves water mass and energy budget equations
to describe water cycle in the Earth´s critical zone.
The package geotopbricks (https://github.com/ecor/geotopbricks and https://CRAN.R-project.org/package=geotopbricks ) is able to to read the GEOtop I/O data of the model. The package geotopOptim2 (https://github.com/EURAC-Ecohydro/geotopOptim2) calling the hydroPSO (https://CRAN.R-project.org/package=hydroPSO) package can be used for model calibration against observations. Further details and complete R package dependencies are listed in geotopOtim2 description file.
As a demonstration example, an analysis of modeled and observed soil moisture and evapotranspiration time series in some alpine agricultural sites (https://github.com/EURAC-Ecohydro/MonaLisa) are presented.
Body row {data-height=75}
GEOtop Hydrological Model
GEOtop (Endrizzi et al, 2014 and references therein) is a distributed model of the mass and energy balance of the hydrological cycle. GEOtop is applicable to simulations in continuum in small or relatively large montain catchments. GEOtop deals with the effects of topography on the interaction between energy balance and hydrological cycle (water, glacier and snow) with peculiar solutions. The source code of GEOtop 2.0 with detailed documentation is available for researchers and environmental/software engineers through the following links:
- https://github.com/geotopmodel/geotop
- https://github.com/se27xx/GEOtop/
heading
This section defines the first column of the body row of the poster.
You can include regular markdown text and R code chunks.
hist(volcano)
geotopbricks
Content for the second column.
To show the code in the chunks you have to ask for it with echo=TRUE.
filled.contour(volcano, color.palette = terrain.colors, asp = 1)
You can also include data frames as tables just by printing them
head(iris)
geotopOptim2
Introduction
The 'geotopOptim2' package has been built to perform a PSO Calibration through hydroPSO package of the the results of the physically-based hydrological model GEOtop.
The package makes use of the GEOtop-R plugin geotopbricks, which import/export GEOtop Data into R.
You need to have the following softwares installed: R, git and GEOtop
GEOtop Hydrological Model
GEOtop (Endrizzi et al, 2014 and references therein) is a distributed model of the mass and energy balance of the hydrological cycle. GEOtop is applicable to simulations in continuum in small or relatively large montain catchments. GEOtop deals with the effects of topography on the interaction between energy balance and hydrological cycle (water, glacier and snow) with peculiar solutions. The source code of GEOtop 2.0 with detailed documentation is available for researchers and environmental/software engineers through the following links:
- https://github.com/geotopmodel/geotop
- https://github.com/se27xx/GEOtop/
So far, goeotpOpintm2 has been usend only with the se27xx version under linux.
For Unix-like OS users, a C source code can be rapidly downloaded and built with the following instrunctions
(https://github.com/se27xx/GEOtop/):
- Open a console and go to the directory where to clone GEOtop source code;
- Clone the source code typing: "git clone https://github.com/se27xx/GEOtop";
- Enter GEOtop directory typing: "cd GEOtop";
- Create the subdirectory for the executable file "mkdir bin";
- Build and create GEOtop executale file: "make -f geotop.make"
GEOtop executable will be created in the subdirectory "bin".
For other versions of GEOtop, please see the related URLs.
It is ongoing the testing for the GEOtop2.1 version supported by a larger community.
Please read detailed information on the installation.
Overview of the main function
GEOtopOptim2 is a package which is able to enter a GEOtop simulation folder and calibarate some parameters (set by user) optimizing a Goodnes-of-Fit measure between point results and point and observations. eotopOptim2 is a set of functions that calls dependency function and the user can work with all or some of these functions as it in his or her preferences. Before entering geotopOptim2, here is a list of the main finctions:
- geotopExec : This function executes a simulation of GEOtop given a certain set of parameters:
library(geotopOptim2)
help(geotopExec,help_type="html")
help(get.geotop.keyword.inpts.value,help_type = "html")
- geotopLookUpTable: It read the lookup table that can be inserted within a GEOtop point simulation according Read point output from GEOtop for verification of the model results:
help(geotopLookUpTable,help_type="html")
To use this function, the observation point value must be set in the GEOtop simulation and inserted as additional keywords in the geotop.inpts file , see the following file for the package example simulation (look at the final OBSERVATION section):
wpath_sim <- system.file('geotop-simulation/B2site',package="geotopOptim2")
inpts.file.path <- paste(wpath_sim,"geotop.inpts",sep="/")
inpts <- readLines(inpts.file.path)
writeLines(inpts)
In the observation section, two now keywords are set (the symbols !>>! means that the keyword is not read by GEOtop!):
ObservationProfileFile, ObservationLookupTblFile. The first is similar to PointFile or SoilLiqContentProfileFile and refers to a set of csv files containing time series of observed variable values.
library(geotopbricks)
obs <- get.geotop.inpts.keyword.value("ObservationProfileFile",wpath=wpath_sim,data.frame=TRUE,date_field='Date12.DDMMYYYYhhmm.')
str(obs)
```
The second one refers to a txt file containing a table with the correspondance between the observation variables and the variables simulated by GEOtop:
(NOTE: here there is no formatter in the file neme and the seperator is *;* , so please set *formatter="",col_sep=";"* instead of default values):
```r
ltk <- get.geotop.inpts.keyword.value("ObservationLookupTblFile",wpath=wpath_sim,data.frame=TRUE,formatter="",col_sep=";")
head(ltk)
knitr::kable(ltk)
This table contains the following columns:
dl>
obs_var
name of the observed variables as reported in the csv file for observations;
geotop_what
GEOtop keyword corresponding for the oitput point file where the variable is written;
geotop_here
name of the corresponding variable simulated by GEOtop as reported in the output csv file;
unit
measurement unit.
The keywords in geotop_where or geotop_what can be more than one and can be operated through simple mathematical operations.
- geotoGOF : It is a wrapper for gof function, it calculates the goodness of fit index for the variables choosen in the data frame returned by geotopLookUpTable; it includes geotopExec and geotopLookukUpTale. Main arguments of geotopGOF are passed to geotopExec as main arguments.
help(geotoGOF,help_type="html")
help(gof,help_type="gof")
- geotopPSO : It is a wrapper for hydroPSO function, it calculates the goodness of fit index for the variables choosen in the data frame returned by geotopLookUpTable:
help(geotoGOF,help_type="html")
help(gof,help_type="gof")
These functions call other functions and methods of geotopOptim2 and its dependecies. Before entering geotopOptim2 in details, the readers' attentions must be focused on the main depency widely used in geotopOptim2:
-
geotopbricks: This package reads the configuration file 'geotop.inpts' of a GEOtop Hydrological Model simulation and get or set information on the date used or generated by the simulation (e.g. soil properties, weather time series, soil moisture and evapotranspiration, etc.);
-
hydroGOF: It is a set of S3 functions implementing both statistical and graphical goodness-of-fit measures between observed and simulated values, mainly oriented to be used during the calibration, validation, and application of hydrological models;
- hydroPSO : "This package implements a state-of-the-art version of the Particle Swarm Optimisation (PSO) algorithm (SPSO-2011 and SPSO-2007 capable). hydroPSO can be used as a replacement of the 'optim' R function for (global) optimization of non-smooth and non-linear functions. However, the main focus of hydroPSO is the calibration of environmental and other real-world models that need to be executed from the system console. hydroPSO is model-independent, allowing the user to easily interface any computer simulation model with the calibration engine (PSO). hydroPSO communicates with the model through the model's own input and output files, without requiring access to the model's source code. Several PSO variants and controlling options are included to fine-tune the performance of the calibration engine to different calibration problems. An advanced sensitivity analysis function together with user-friendly plotting summaries facilitate the interpretation and assessment of the calibration results. hydroPSO is parallel-capable, to alleviate the computational burden of complex models with "long" execution time. Bugs reports/comments/questions are very welcomed (in English, Spanish or Italian)." (Zambrano & Rojas, 2013)
A package can be drawn as a citation graph in whch functions depend on others. In the following the code lines can reproduce the dependency graph of geotopOptim2 package and the main functions belong to dependency packages.
#!/usr/bin/env Rscript
# file appendSmetData.R
#
# This script creates a graph of the package function and thair main external depencies
#
# author: Emanuele Cordano on 09-09-2015
#This program is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.
#
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with this program. If not, see <http://www.gnu.org/licenses/>.
###############################################################################
rm(list=ls())
library(geotopOptim2)
library(igraph)
set.seed(123)
list_envs <- list(environment(geotopExec),environment(get.geotop.inpts.keyword.value),environment(gof),environment(hydroPSO),environment(writeLines),environment(read.table),environment(terrain))
names(list_envs) <- c("geotopOptim2","geotopbricks","hydroGOF","hydroPSO","base","utils","raster")
color <- c("red","green","blue","orange","yellow","white","brown")
names(color) <- names(list_envs)
list_names <- lapply(X=list_envs,FUN=function(x){ls(env=x)})
list_df <- list()
for (it in names(list_envs)) {
list_df[[it]] <- data.frame(funx=list_names[[it]],env=it,color=color[it],stringsAsFactors=FALSE)
}
df <- do.call(what=rbind,args=list_df)
##### SEMPLIFICATE DF
onlyfun <- list(hydroGOF=c("gof"),hydroPSO=c("hydroPSO","lhoat"),geotopbricks=c("get.geotop.inpts.keyword.value"),
base=c("writeLines","readLines"),utils=c("read.table"),raster="raster") ### read.table was removed
for (it in names(onlyfun)) {
cond <- ((df$env==it) & (df$funx %in% onlyfun[[it]])) | (df$env!=it)
df <- df[cond,]
}
fun_names <- df$funx
names(fun_names) <- fun_names
########################################
########################################
########################################
########################################
########################################
lfunx <- lapply(X=fun_names,FUN=function(x,nx) {
o <- try(get(x),silent=TRUE)
if (class(o)=="try-error") {
o <- NA ### "It looks like a method!"
return(o)
}
o <- formals(o)
o <- lapply(X=o,FUN=as.character)
o <- unlist(o)
o <- o[o %in% nx]
src <- as.character(body(x))
src <- unlist(str_split(src, boundary("word")))
nx <- src[src %in% nx]
o <- c(o,nx)
o <- unique(o)
return(o)
},nx=fun_names)
for (it in names(lfunx)) {
temp <- lfunx[[it]]
ii <- which(temp!=it)
temp <- temp[ii]
nl <- length(temp)
lfunx[[it]] <- array(c(rep(it,nl),temp),c(nl,2))
}
#####edges
edges <- do.call(rbind,lfunx)
vertices <- unique(edges)
#####
env_base <- "base;utils"
df$env[df$env=="base"] <- env_base
df$env[df$env=="utils"] <- env_base
df$color[df$env==env_base] <- "white"
#####
color_ <- df$color
env_ <- df$env
names(color_) <- df$funx
names(env_) <- df$funx
######
gg <- graph_from_edgelist(edges)
vnames <- V(gg)$name
V(gg)$color <- color_[vnames]
vcodes <- sprintf("%02d",1:length(vnames))
names(vcodes) <- vnames
V(gg)$name <- vcodes
main <- "geeotopOptim2 Internal Functions"
plot(gg,main=main)
legend("bottomleft",legend=unique(env_),fill=unique(color_),ncol=2)
legend("topleft",legend=paste(vcodes,vnames,sep=" : "),ncol=3,cex=0.6)
The graph show
geotopExec
geotopGOF
geotopLook
geotopPSO
Note the various macros within the vignette section of the metadata block above. These are required in order to instruct R how to build the vignette. Note that you should change the title field and the \VignetteIndexEntry to match the title of your vignette.
Exercise
#!/usr/bin/env Rscript
# file pso_example_script.R
#
# This script is an examples of a GEOtop calibration via geotopOptim2
#
# author: Emanuele Cordano on 09-09-2015
#This program is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.
#
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with this program. If not, see <http://www.gnu.org/licenses/>.
###############################################################################
rm(list=ls())
library(zoo)
library(geotopOptim2)
set.seed(7988)
USE_RMPI <- FALSE
if (USE_RMPI==TRUE) {
library("parallel")
library(Rmpi)
require(snow)
if (mpi.comm.rank(0) > 0) {
sink(file="/dev/null")
#runMPIslave()
slaveLoop(makeMPImaster())
mpi.quit()
}
parallel <- "parallel"
npart <- 16
control <- list(maxit=5,npart=npart,parallel=parallel)
} else {
parellel <- "none"
npart <- 4
control <- list(maxit=2,npart=npart)
}
Here is a GEOtop simulation folder with observation:
tz <- "Etc/GMT-1"
wpath <- system.file('geotop-simulation/B2site',package="geotopOptim2")
Here is the full name (with complete path) of GEOtop executable built file:
bin <-'/home/ecor/local/geotop/GEOtop/bin/geotop-2.0.0'
In most cases, geotopOptim2 copies the simulation folder and run GEOtop with modified calibration input parameters in temporary directories:
## Local path where to write output for PSO
runpath <- "/home/lv70864/ecordano/temp/geotopOptim_tests"
runpath <- "/home/ecor/temp/geotopOptim_tests"
The function geotopExec running a GEOtop simulation can enter a vector of paratemers which are replaced with the ones formerly set in the original simulation directory wpath . Analogously the function geotopPSO performing a PSO calibration of GEOtop through geotopExec and geotopGOF needs two vectors: one for upper values, the other for lower values. A set of soil calibration parameters can be set though a csv file.
help(geotopExec,help_type="html")
help(geotopGOF,help_type="html")
help(geotopPSO,help_type="html")
Here is an example that can be downloaded by package examples:
geotop.soil.param.file <- system.file('examples_script/param/param_pso_c003.csv',package="geotopOptim2")
geotop.soil.param <- read.table(geotop.soil.param.file,header=TRUE,sep=",",stringsAsFactors=FALSE)
lower <- geotop.soil.param$lower
upper <- geotop.soil.param$upper
suggested <- geotop.soil.param$suggested
names(lower) <- geotop.soil.param$name
names(upper) <- geotop.soil.param$name
names(suggested) <- geotop.soil.param$name
knitr::kable(geotop.soil.param)
Here is the columns lower and upper stand for lower and upper values whereas the suggested column contains initial guess values for geotopPSO. Generally a good practical examples on how to set these parameters is given below:
prefix__name,lower,upperGrenarally
SOIL__N,1.45,1.89
SOIL__Alpha,0.00131,0.0131
SOIL__ThetaSat,0.41,0.53
SOIL__ThetaSat_bottomlayer,0.08,0.09
SOIL__ThetaRes,0.05,0.07
SOIL__LateralHydrConductivity,0.0923,0.1
SOIL__NormalHydrConductivity,0.0923,0.1
SOIL__LateralHydrConductivity_bottomlayer,0.00923,0.01
SOIL__NormalHydrConductivity_bottomlayer,0.00923,0.01
SOIL__SoilInitPresL0001,-10000,100
SOIL__SoilInitPresL0002,-10000,100
SOIL__SoilInitPresL0003,-10000,100
SOIL__SoilInitPresL0004,-10000,100
SOIL__SoilInitPresL0005,-10000,100
SOIL__SoilInitPresL0006,-10000,100
SOIL__SoilInitPresL0007,-10000,1000
SOIL__SoilInitPresL0008,-10000,1000
SOIL__SoilInitPresL0009,-10000,0
SOIL__PsiGamma,0.5,1
SOIL__SoilDepth,3000,30000
SOIL__NumberOfSoilLayers,9,20
VECTOR_1_LSAI,2,4
where "N", "Alpha", "ThetaSat", "LateralHydrConductivity", "NormalHydrConductivity", "ThetaRes",,.. are GEOtop keyword referred to the respective soil parameters. By default, "geotopPSO" considers soil parameters uniformly distributed within the soil profile unless they are repated in the CSV file with some suffixes, like "_bottomlayer" or "_ALL". In "_bottomlayer" case, the parameter (upper and lower) values are referred to the first (near surface) layer and the last (bottom) layer and the values of internal layers are exponentially interpolated. In this case a decrease of hydraulic conductivity or soil porosity can be modeled. In "_ALL" case, soil parameter is taken as variables with soil layers. So the reported values refer to a range for so many soil parameter of the same type how many the soil layers are. If the keywords contains the suffix "_V_L%04d" with the decimal formatter, the soil parameter is calibrated only for the soil specific layer.
The keyword "SoilInitPres" refers to the initial soil water pressure head. If the formatter "L%04d" is appended as a suffix, the value is referred to the indicated layer. The keyword "PsiGamma" refers to the soil water pressure gradient along the terrain-normal downward direction and is applied to calculate initial soil water pressure head in the above layers assuming a continuous profile.
"SoilDepth" and "NumberOfSoilLayers" refer to the whole soil depth (which now corresponds to the whole soil column used as domain for balance equation integration) and the number of soil layers in which the soil column is divided. Soil layer thickness increase with depth following a geometric progression.
Then,the target variables are here set :
var <- 'soil_moisture_content_50'
x <- (upper+lower)/2
pso <- geotopPSO(par=suggested,run.geotop=TRUE,bin=bin,
simpath=wpath,runpath=runpath,clean=TRUE,data.frame=TRUE,
level=1,intern=TRUE,target=var,gof.mes="RMSE",lower=lower,upper=upper,control=control)
if (USE_RMPI==TRUE) mpi.finalize()
Note the various macros within the vignette section of the metadata block above. These are required in order to instruct R how to build the vignette. Note that you should change the title field and the \VignetteIndexEntry to match the title of your vignette.
Examples:
Examples of simulations to optimize could be found here https://github.com/EURAC-Ecohydro/geotopOptim2/tree/master/inst/examples_script
Examples of scripts to launch optimization in background on a cluster for the EURAC Monalisa dataset could be found here:
https://github.com/EURAC-Ecohydro/MonaLisa/tree/master/Rscript
Visulatization tools:
R-based graphical visualization tool of the optimized simulations could be found here:
https://github.com/EURAC-Ecohydro/geotopOptim2/tree/master/inst/examples_animation
Footer row {data-height=10}
{data-width=50}
{data-width=50 .small}
EURAC-Ecohydro/geotopOptim2 documentation built on March 3, 2021, 4:56 a.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.
Header row {data-height=15}
**Poster title** |
|
Abstract Eco-hydrological models are increasingly used in the contexts of hydrology, ecology, precision agriculture for a detailed description of the water cycle at various scales: local scale, an hillslope or a watershed. However, with increasing computing power and observations available, bigger and bigger amount of raw data are produced. Therefore the need to develop flexible and user-oriented interfaces to visualize multiple outputs, perform sensitivity analyzes and compare against observations emerges. This work presents two R open-source packages: geotopbricks and geotopOptim2. They offer an I/0 interface and R visualization tools the GEOtop hydrological distributed model (http://geotopmodel.github.io/geotop/), which solves water mass and energy budget equations to describe water cycle in the Earth´s critical zone. The package geotopbricks (https://github.com/ecor/geotopbricks and https://CRAN.R-project.org/package=geotopbricks ) is able to to read the GEOtop I/O data of the model. The package geotopOptim2 (https://github.com/EURAC-Ecohydro/geotopOptim2) calling the hydroPSO (https://CRAN.R-project.org/package=hydroPSO) package can be used for model calibration against observations. Further details and complete R package dependencies are listed in geotopOtim2 description file. As a demonstration example, an analysis of modeled and observed soil moisture and evapotranspiration time series in some alpine agricultural sites (https://github.com/EURAC-Ecohydro/MonaLisa) are presented.
Body row {data-height=75}
GEOtop Hydrological Model
GEOtop (Endrizzi et al, 2014 and references therein) is a distributed model of the mass and energy balance of the hydrological cycle. GEOtop is applicable to simulations in continuum in small or relatively large montain catchments. GEOtop deals with the effects of topography on the interaction between energy balance and hydrological cycle (water, glacier and snow) with peculiar solutions. The source code of GEOtop 2.0 with detailed documentation is available for researchers and environmental/software engineers through the following links:
- https://github.com/geotopmodel/geotop
- https://github.com/se27xx/GEOtop/
heading
This section defines the first column of the body row of the poster.
You can include regular markdown text and R code chunks.
hist(volcano)
geotopbricks
Content for the second column.
To show the code in the chunks you have to ask for it with echo=TRUE.
filled.contour(volcano, color.palette = terrain.colors, asp = 1)
You can also include data frames as tables just by printing them
head(iris)
geotopOptim2
Introduction
The 'geotopOptim2' package has been built to perform a PSO Calibration through hydroPSO package of the the results of the physically-based hydrological model GEOtop. The package makes use of the GEOtop-R plugin geotopbricks, which import/export GEOtop Data into R.
You need to have the following softwares installed: R, git and GEOtop
GEOtop Hydrological Model
GEOtop (Endrizzi et al, 2014 and references therein) is a distributed model of the mass and energy balance of the hydrological cycle. GEOtop is applicable to simulations in continuum in small or relatively large montain catchments. GEOtop deals with the effects of topography on the interaction between energy balance and hydrological cycle (water, glacier and snow) with peculiar solutions. The source code of GEOtop 2.0 with detailed documentation is available for researchers and environmental/software engineers through the following links:
- https://github.com/geotopmodel/geotop
- https://github.com/se27xx/GEOtop/
So far, goeotpOpintm2 has been usend only with the se27xx version under linux.
For Unix-like OS users, a C source code can be rapidly downloaded and built with the following instrunctions (https://github.com/se27xx/GEOtop/):
- Open a console and go to the directory where to clone GEOtop source code;
- Clone the source code typing: "git clone https://github.com/se27xx/GEOtop";
- Enter GEOtop directory typing: "cd GEOtop";
- Create the subdirectory for the executable file "mkdir bin";
- Build and create GEOtop executale file: "make -f geotop.make"
GEOtop executable will be created in the subdirectory "bin". For other versions of GEOtop, please see the related URLs.
It is ongoing the testing for the GEOtop2.1 version supported by a larger community. Please read detailed information on the installation.
Overview of the main function
GEOtopOptim2 is a package which is able to enter a GEOtop simulation folder and calibarate some parameters (set by user) optimizing a Goodnes-of-Fit measure between point results and point and observations. eotopOptim2 is a set of functions that calls dependency function and the user can work with all or some of these functions as it in his or her preferences. Before entering geotopOptim2, here is a list of the main finctions:
- geotopExec : This function executes a simulation of GEOtop given a certain set of parameters:
library(geotopOptim2) help(geotopExec,help_type="html") help(get.geotop.keyword.inpts.value,help_type = "html")
- geotopLookUpTable: It read the lookup table that can be inserted within a GEOtop point simulation according Read point output from GEOtop for verification of the model results:
help(geotopLookUpTable,help_type="html")
To use this function, the observation point value must be set in the GEOtop simulation and inserted as additional keywords in the geotop.inpts file , see the following file for the package example simulation (look at the final OBSERVATION section):
wpath_sim <- system.file('geotop-simulation/B2site',package="geotopOptim2") inpts.file.path <- paste(wpath_sim,"geotop.inpts",sep="/") inpts <- readLines(inpts.file.path) writeLines(inpts)
In the observation section, two now keywords are set (the symbols !>>! means that the keyword is not read by GEOtop!): ObservationProfileFile, ObservationLookupTblFile. The first is similar to PointFile or SoilLiqContentProfileFile and refers to a set of csv files containing time series of observed variable values.
library(geotopbricks) obs <- get.geotop.inpts.keyword.value("ObservationProfileFile",wpath=wpath_sim,data.frame=TRUE,date_field='Date12.DDMMYYYYhhmm.') str(obs) ``` The second one refers to a txt file containing a table with the correspondance between the observation variables and the variables simulated by GEOtop: (NOTE: here there is no formatter in the file neme and the seperator is *;* , so please set *formatter="",col_sep=";"* instead of default values): ```r ltk <- get.geotop.inpts.keyword.value("ObservationLookupTblFile",wpath=wpath_sim,data.frame=TRUE,formatter="",col_sep=";") head(ltk) knitr::kable(ltk)
This table contains the following columns: dl>
The keywords in geotop_where or geotop_what can be more than one and can be operated through simple mathematical operations.
- geotoGOF : It is a wrapper for gof function, it calculates the goodness of fit index for the variables choosen in the data frame returned by geotopLookUpTable; it includes geotopExec and geotopLookukUpTale. Main arguments of geotopGOF are passed to geotopExec as main arguments.
help(geotoGOF,help_type="html") help(gof,help_type="gof")
- geotopPSO : It is a wrapper for hydroPSO function, it calculates the goodness of fit index for the variables choosen in the data frame returned by geotopLookUpTable:
help(geotoGOF,help_type="html") help(gof,help_type="gof")
These functions call other functions and methods of geotopOptim2 and its dependecies. Before entering geotopOptim2 in details, the readers' attentions must be focused on the main depency widely used in geotopOptim2:
-
geotopbricks: This package reads the configuration file 'geotop.inpts' of a GEOtop Hydrological Model simulation and get or set information on the date used or generated by the simulation (e.g. soil properties, weather time series, soil moisture and evapotranspiration, etc.);
-
hydroGOF: It is a set of S3 functions implementing both statistical and graphical goodness-of-fit measures between observed and simulated values, mainly oriented to be used during the calibration, validation, and application of hydrological models;
- hydroPSO : "This package implements a state-of-the-art version of the Particle Swarm Optimisation (PSO) algorithm (SPSO-2011 and SPSO-2007 capable). hydroPSO can be used as a replacement of the 'optim' R function for (global) optimization of non-smooth and non-linear functions. However, the main focus of hydroPSO is the calibration of environmental and other real-world models that need to be executed from the system console. hydroPSO is model-independent, allowing the user to easily interface any computer simulation model with the calibration engine (PSO). hydroPSO communicates with the model through the model's own input and output files, without requiring access to the model's source code. Several PSO variants and controlling options are included to fine-tune the performance of the calibration engine to different calibration problems. An advanced sensitivity analysis function together with user-friendly plotting summaries facilitate the interpretation and assessment of the calibration results. hydroPSO is parallel-capable, to alleviate the computational burden of complex models with "long" execution time. Bugs reports/comments/questions are very welcomed (in English, Spanish or Italian)." (Zambrano & Rojas, 2013)
A package can be drawn as a citation graph in whch functions depend on others. In the following the code lines can reproduce the dependency graph of geotopOptim2 package and the main functions belong to dependency packages.
#!/usr/bin/env Rscript # file appendSmetData.R # # This script creates a graph of the package function and thair main external depencies # # author: Emanuele Cordano on 09-09-2015 #This program is free software: you can redistribute it and/or modify #it under the terms of the GNU General Public License as published by #the Free Software Foundation, either version 3 of the License, or #(at your option) any later version. # #This program is distributed in the hope that it will be useful, #but WITHOUT ANY WARRANTY; without even the implied warranty of #MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #GNU General Public License for more details. # #You should have received a copy of the GNU General Public License #along with this program. If not, see <http://www.gnu.org/licenses/>. ############################################################################### rm(list=ls()) library(geotopOptim2) library(igraph) set.seed(123) list_envs <- list(environment(geotopExec),environment(get.geotop.inpts.keyword.value),environment(gof),environment(hydroPSO),environment(writeLines),environment(read.table),environment(terrain)) names(list_envs) <- c("geotopOptim2","geotopbricks","hydroGOF","hydroPSO","base","utils","raster") color <- c("red","green","blue","orange","yellow","white","brown") names(color) <- names(list_envs) list_names <- lapply(X=list_envs,FUN=function(x){ls(env=x)}) list_df <- list() for (it in names(list_envs)) { list_df[[it]] <- data.frame(funx=list_names[[it]],env=it,color=color[it],stringsAsFactors=FALSE) } df <- do.call(what=rbind,args=list_df) ##### SEMPLIFICATE DF onlyfun <- list(hydroGOF=c("gof"),hydroPSO=c("hydroPSO","lhoat"),geotopbricks=c("get.geotop.inpts.keyword.value"), base=c("writeLines","readLines"),utils=c("read.table"),raster="raster") ### read.table was removed for (it in names(onlyfun)) { cond <- ((df$env==it) & (df$funx %in% onlyfun[[it]])) | (df$env!=it) df <- df[cond,] } fun_names <- df$funx names(fun_names) <- fun_names ######################################## ######################################## ######################################## ######################################## ######################################## lfunx <- lapply(X=fun_names,FUN=function(x,nx) { o <- try(get(x),silent=TRUE) if (class(o)=="try-error") { o <- NA ### "It looks like a method!" return(o) } o <- formals(o) o <- lapply(X=o,FUN=as.character) o <- unlist(o) o <- o[o %in% nx] src <- as.character(body(x)) src <- unlist(str_split(src, boundary("word"))) nx <- src[src %in% nx] o <- c(o,nx) o <- unique(o) return(o) },nx=fun_names) for (it in names(lfunx)) { temp <- lfunx[[it]] ii <- which(temp!=it) temp <- temp[ii] nl <- length(temp) lfunx[[it]] <- array(c(rep(it,nl),temp),c(nl,2)) } #####edges edges <- do.call(rbind,lfunx) vertices <- unique(edges) ##### env_base <- "base;utils" df$env[df$env=="base"] <- env_base df$env[df$env=="utils"] <- env_base df$color[df$env==env_base] <- "white" ##### color_ <- df$color env_ <- df$env names(color_) <- df$funx names(env_) <- df$funx ###### gg <- graph_from_edgelist(edges) vnames <- V(gg)$name V(gg)$color <- color_[vnames] vcodes <- sprintf("%02d",1:length(vnames)) names(vcodes) <- vnames V(gg)$name <- vcodes main <- "geeotopOptim2 Internal Functions" plot(gg,main=main) legend("bottomleft",legend=unique(env_),fill=unique(color_),ncol=2) legend("topleft",legend=paste(vcodes,vnames,sep=" : "),ncol=3,cex=0.6)
The graph show
geotopExec geotopGOF geotopLook geotopPSO
Note the various macros within the vignette section of the metadata block above. These are required in order to instruct R how to build the vignette. Note that you should change the title field and the \VignetteIndexEntry to match the title of your vignette.
Exercise
#!/usr/bin/env Rscript # file pso_example_script.R # # This script is an examples of a GEOtop calibration via geotopOptim2 # # author: Emanuele Cordano on 09-09-2015 #This program is free software: you can redistribute it and/or modify #it under the terms of the GNU General Public License as published by #the Free Software Foundation, either version 3 of the License, or #(at your option) any later version. # #This program is distributed in the hope that it will be useful, #but WITHOUT ANY WARRANTY; without even the implied warranty of #MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #GNU General Public License for more details. # #You should have received a copy of the GNU General Public License #along with this program. If not, see <http://www.gnu.org/licenses/>. ############################################################################### rm(list=ls()) library(zoo) library(geotopOptim2) set.seed(7988) USE_RMPI <- FALSE if (USE_RMPI==TRUE) { library("parallel") library(Rmpi) require(snow) if (mpi.comm.rank(0) > 0) { sink(file="/dev/null") #runMPIslave() slaveLoop(makeMPImaster()) mpi.quit() } parallel <- "parallel" npart <- 16 control <- list(maxit=5,npart=npart,parallel=parallel) } else { parellel <- "none" npart <- 4 control <- list(maxit=2,npart=npart) }
Here is a GEOtop simulation folder with observation:
tz <- "Etc/GMT-1" wpath <- system.file('geotop-simulation/B2site',package="geotopOptim2")
Here is the full name (with complete path) of GEOtop executable built file:
bin <-'/home/ecor/local/geotop/GEOtop/bin/geotop-2.0.0'
In most cases, geotopOptim2 copies the simulation folder and run GEOtop with modified calibration input parameters in temporary directories:
## Local path where to write output for PSO runpath <- "/home/lv70864/ecordano/temp/geotopOptim_tests" runpath <- "/home/ecor/temp/geotopOptim_tests"
The function geotopExec running a GEOtop simulation can enter a vector of paratemers which are replaced with the ones formerly set in the original simulation directory wpath . Analogously the function geotopPSO performing a PSO calibration of GEOtop through geotopExec and geotopGOF needs two vectors: one for upper values, the other for lower values. A set of soil calibration parameters can be set though a csv file.
help(geotopExec,help_type="html") help(geotopGOF,help_type="html") help(geotopPSO,help_type="html")
Here is an example that can be downloaded by package examples:
geotop.soil.param.file <- system.file('examples_script/param/param_pso_c003.csv',package="geotopOptim2") geotop.soil.param <- read.table(geotop.soil.param.file,header=TRUE,sep=",",stringsAsFactors=FALSE) lower <- geotop.soil.param$lower upper <- geotop.soil.param$upper suggested <- geotop.soil.param$suggested names(lower) <- geotop.soil.param$name names(upper) <- geotop.soil.param$name names(suggested) <- geotop.soil.param$name knitr::kable(geotop.soil.param)
Here is the columns lower and upper stand for lower and upper values whereas the suggested column contains initial guess values for geotopPSO. Generally a good practical examples on how to set these parameters is given below:
prefix__name,lower,upperGrenarally SOIL__N,1.45,1.89 SOIL__Alpha,0.00131,0.0131 SOIL__ThetaSat,0.41,0.53 SOIL__ThetaSat_bottomlayer,0.08,0.09 SOIL__ThetaRes,0.05,0.07 SOIL__LateralHydrConductivity,0.0923,0.1 SOIL__NormalHydrConductivity,0.0923,0.1 SOIL__LateralHydrConductivity_bottomlayer,0.00923,0.01 SOIL__NormalHydrConductivity_bottomlayer,0.00923,0.01 SOIL__SoilInitPresL0001,-10000,100 SOIL__SoilInitPresL0002,-10000,100 SOIL__SoilInitPresL0003,-10000,100 SOIL__SoilInitPresL0004,-10000,100 SOIL__SoilInitPresL0005,-10000,100 SOIL__SoilInitPresL0006,-10000,100 SOIL__SoilInitPresL0007,-10000,1000 SOIL__SoilInitPresL0008,-10000,1000 SOIL__SoilInitPresL0009,-10000,0 SOIL__PsiGamma,0.5,1 SOIL__SoilDepth,3000,30000 SOIL__NumberOfSoilLayers,9,20 VECTOR_1_LSAI,2,4
where "N", "Alpha", "ThetaSat", "LateralHydrConductivity", "NormalHydrConductivity", "ThetaRes",,.. are GEOtop keyword referred to the respective soil parameters. By default, "geotopPSO" considers soil parameters uniformly distributed within the soil profile unless they are repated in the CSV file with some suffixes, like "_bottomlayer" or "_ALL". In "_bottomlayer" case, the parameter (upper and lower) values are referred to the first (near surface) layer and the last (bottom) layer and the values of internal layers are exponentially interpolated. In this case a decrease of hydraulic conductivity or soil porosity can be modeled. In "_ALL" case, soil parameter is taken as variables with soil layers. So the reported values refer to a range for so many soil parameter of the same type how many the soil layers are. If the keywords contains the suffix "_V_L%04d" with the decimal formatter, the soil parameter is calibrated only for the soil specific layer. The keyword "SoilInitPres" refers to the initial soil water pressure head. If the formatter "L%04d" is appended as a suffix, the value is referred to the indicated layer. The keyword "PsiGamma" refers to the soil water pressure gradient along the terrain-normal downward direction and is applied to calculate initial soil water pressure head in the above layers assuming a continuous profile. "SoilDepth" and "NumberOfSoilLayers" refer to the whole soil depth (which now corresponds to the whole soil column used as domain for balance equation integration) and the number of soil layers in which the soil column is divided. Soil layer thickness increase with depth following a geometric progression.
Then,the target variables are here set :
var <- 'soil_moisture_content_50' x <- (upper+lower)/2
pso <- geotopPSO(par=suggested,run.geotop=TRUE,bin=bin, simpath=wpath,runpath=runpath,clean=TRUE,data.frame=TRUE, level=1,intern=TRUE,target=var,gof.mes="RMSE",lower=lower,upper=upper,control=control) if (USE_RMPI==TRUE) mpi.finalize()
Note the various macros within the vignette section of the metadata block above. These are required in order to instruct R how to build the vignette. Note that you should change the title field and the \VignetteIndexEntry to match the title of your vignette.
Examples:
Examples of simulations to optimize could be found here https://github.com/EURAC-Ecohydro/geotopOptim2/tree/master/inst/examples_script
Examples of scripts to launch optimization in background on a cluster for the EURAC Monalisa dataset could be found here: https://github.com/EURAC-Ecohydro/MonaLisa/tree/master/Rscript
Visulatization tools: R-based graphical visualization tool of the optimized simulations could be found here: https://github.com/EURAC-Ecohydro/geotopOptim2/tree/master/inst/examples_animation
Footer row {data-height=10}
{data-width=50}
{data-width=50 .small}
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.