# lumpR/db_prepare_musle.R
# Copyright (C) 2019 Till Francke
#
# 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/>.
#' Prepare various input related to use of (M)USLE module (erosion)
#'
#' Function to write parameter values relevant for modelling application with the
#' WASA hydrological model into an existing database, preferably created with
#' \code{\link[lumpR]{db_create}}.
#'
#' @param dbname Name of the data source (DSN) registered at ODBC. See \code{Details} of
#' \code{\link[lumpR]{db_create}}.
#'
#' @param compute_K Compute (M)USLE K-factor (soil erodibility) from properties of topmost soil horizon, write it to table "soil" (use \code{copy_from_other_tables} to transfer it to table "soil_veg_components"). \code{c(TRUE, FALSE)}
#'
#' @param copy_from_other_tables Copy values to table "soil_veg_components": MUSLE-C and Manning-n from table "vegetation", MUSLE-K from "soils". Array of strings with possible values \code{c("MUSLE-C", "Manning-n" ,"MUSLE-K")} of
#' \code{\link[lumpR]{db_create}}.
#'
#' @param setP Set entire column of (M)USLE P-factor (protection measures) to a single value (usually 1). Default NULL (do not set).
#' @param verbose \code{logical}. Should detailed information during execution be
#' printed? Default: \code{TRUE}.
#'
#' @details
#' after Williams (1995). [explanation to be elaborated]
#'
#'
#' @references
#' \bold{lumpR paper}\cr
#' Pilz, T.; Francke, T.; Bronstert, A. (2017): lumpR 2.0.0: an R package facilitating
#' landscape discretisation for hillslope-based hydrological models.
#' \emph{Geosci. Model Dev.}, 10, 3001-3023, doi: 10.5194/gmd-10-3001-2017
#'
#'
#'
#'
#' @author
#' Till Francke \email{francke@@uni-potsdam.de}
#'
db_prepare_musle <- function(
dbname,
compute_K=FALSE,
copy_from_other_tables=NULL,
setP=NULL,
verbose=TRUE
) {
#this is a quick-and-dirty conversion from legacy PHP-code. It could probaby be solved in a couple of lines in R. Sorry for not having had the time, though.
if(verbose) message("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
if(verbose) message("% START db_prepare_musle()")
if(verbose) message("%")
if(verbose) message("% Connecting to and checking database ...")
# connect to ODBC registered database
con <- connect_db(dbname)
#modify error handler to gracefully close ODBC-connection before aborting (otherwise, ODBC-handles are left open)
org_error_handler = getOption("error") #original error handler
closeODBC_and_stop = function()
{
odbcCloseAll()
options(error=org_error_handler) #restore original handler
}
options(error=closeODBC_and_stop) #modify error handler
# ensure MySQL/MariaDB uses ANSI quotation (double quotes instead of back ticks)
if(grepl("MariaDB", odbcGetInfo(con)["DBMS_Name"], ignore.case=T))
sqlQuery(con, "SET sql_mode='ANSI';")
# check database (all 'tables' have to exist)
tables=c("soils", "horizons", "soil_veg_components", "vegetation")
check_tbl <- tables %in% sqlTables(con)$TABLE_NAME
if (any(!check_tbl))
stop(paste0("Table(s) ", paste(tables[which(!check_tbl)], sep=", "), " could not be found in database '", dbname, "'."))
if(verbose) message("% OK")
if (compute_K)
{
#these are the particle size classes that are converted to (USDA-soil classification), used by Williams (1995)
clay_upper_limit=0.002
silt_upper_limit=0.05
sand_upper_limit=2.0
res = sqlQuery2(con, "select * from r_soil_contains_particles", info="check <r_soil_contains_particles>")
if(nrow(res)==0)
stop("Table <r_soil_contains_particles> is empty. Cannot compute K-factor.")
res = sqlQuery2(con, "select * from soils", info="check <soils>")
if(nrow(res)==0)
stop("Table <soils> is empty. Cannot compute K-factor.")
res = sqlQuery2(con, "select * from soil_veg_components", info="check <soil_veg_components>")
if(nrow(res)==0)
stop("Table <soil_veg_components> is empty. Cannot compute K-factor.")
res = sqlQuery2(con, "select * from particle_classes", info="check <particle_classes>")
if(nrow(res)==0)
stop("Table <particle_classes> is empty. Cannot compute K-factor.")
#### compute clay fraction according to USDA-----------------------------------------------------
statement = paste0("select class_id, upper_limit from particle_classes where upper_limit>=",clay_upper_limit, " order by upper_limit")
#statement = paste0("select * from particle_classes;")
#get ID of class that is at or just over usda clay
res = sqlQuery2(con, statement, info="get class-ID")
if(nrow(res)==0)
stop("In table <particle_classes> all specified particle classes are finer than USDA-clay. Cannot compute K-factor.")
class_above_usda_clay = res[1, "class_id" ] #get the id of the user-defined class that sits just above USDA-clay
class_above_usda_clay_limit = res[1, "upper_limit"] #get uppper limit of respective class
statement = "DROP TABLE t_cum_above"
res = sqlQuery(con, query = statement, errors = FALSE) #delete any existing table
statement = paste0("SELECT soil_id, sum(fraction) AS a_cum_above INTO t_cum_above
FROM r_soil_contains_particles WHERE class_id<=", class_above_usda_clay, " GROUP BY soil_id")
res = sqlQuery2(con, statement, info="create temporary table <t_cum_above>") #produce a table containing the cumulative fractions up to class_above_usda_clay for each soil
statement = "DROP TABLE t_cum_below"
res = sqlQuery(con, query = statement, errors = FALSE) #delete any existing table
class_below_usda_clay=class_above_usda_clay-1 #the class below USDA clay
if (class_below_usda_clay==0) #no lower classes that can be used as a point for interpolation
{
class_below_usda_clay_limit=0 #interpolation starts at 0
statement="SELECT soil_id, 0 AS a_cum_below INTO t_cum_below
FROM r_soil_contains_particles GROUP BY soil_id" #statement-statement that produces a table containing zeros for each soil
} else #use the nearest lower class as a point for interpolation
{
statement = paste0("select upper_limit from particle_classes where class_id<=",class_below_usda_clay)
res = sqlQuery2(con, statement, info="get class limit")
if (nrow(res) == 0)
stop("USER-Class below USDA-clay not found.")
class_below_usda_clay_limit=res[1,"upper_limit"] #get uppper limit of respective class
statement=paste0("SELECT soil_id, sum(fraction) AS a_cum_below INTO t_cum_below",
" FROM r_soil_contains_particles WHERE class_id<=",class_below_usda_clay,
" GROUP BY soil_id") #statement-statement that produces a table containing the cumulative fractions up to class_below_usda_clay for each soil
}
res = sqlQuery2(con, statement, info="create temporary table <t_cum_below>") #produce a table containing the cumulative fractions up to class_below_usda_clay for each soil
statement = "select count(*) from soils"
res = sqlQuery2(con, statement, info="checking entries in <soil>")
if(res[1,1]==0)
stop("No rows found in table <soils>.")
statement = "UPDATE soils SET a_clay=0,a_silt=0,a_sand=0"
res = sqlQuery2(con, statement, info="initialise USDA-fractions columns in <soil>")
print("computing USDA-clay-content...")
statement = paste0("UPDATE (soils LEFT JOIN t_cum_above ON soils.pid=t_cum_above.soil_id) LEFT JOIN t_cum_below ON soils.pid=t_cum_below.soil_id",
" SET a_clay = ((a_cum_above-a_cum_below)*(",clay_upper_limit-class_below_usda_clay_limit,"))/(",class_above_usda_clay_limit-class_below_usda_clay_limit,")")
res = sqlQuery2(con, statement, info="compute USDA-clay content")
print("OK")
####compute silt fraction according to USDA-----------------------------------------------------
statement = paste0("select class_id, upper_limit from particle_classes where upper_limit>=", silt_upper_limit," order by upper_limit")
#get ID of class that is at or just over usda silt
res = sqlQuery2(con, statement, info="ID of class that is at or just over usda silt")
if (nrow(res) == 0)
stop("All specified particle classes are finer than USDA-silt, cannot compute K-factor.")
class_above_usda_silt =res[1,"class_id"] #get the user-class that sits just above USDA-silt
class_above_usda_silt_limit=res[1,"upper_limit"] #get uppper limit of respective class
statement = "DROP TABLE t_cum_above"
res = sqlQuery(con, query = statement, errors = FALSE) #delete any existing table
statement = paste0("SELECT soil_id, sum(fraction) AS a_cum_above INTO t_cum_above",
" FROM r_soil_contains_particles WHERE class_id<=",class_above_usda_silt," GROUP BY soil_id")
res = sqlQuery2(con, statement, info="t_cum_above for silt") #produce a table containing the cumulative fractions up to class_above_usda_silt for each soil
statement = "DROP TABLE t_cum_below"
res = sqlQuery(con, query = statement, errors = FALSE) #delete any existing table
class_below_usda_silt=class_above_usda_silt-1 #the class below USDA silt
if (class_below_usda_silt==0) #no lower classes that can be used as a point for interpolation
{
class_below_usda_silt_limit=0 #interpolation starts at 0
statement=paste0("SELECT soil_id, 0 AS a_cum_below INTO t_cum_below",
" FROM r_soil_contains_particles GROUP BY soil_id") #statement-statement that produces a table containing zeros for each soil
} else #use the nearest lower class as a point for interpolation
{
statement = paste0("select upper_limit from particle_classes where class_id=",class_below_usda_silt) #?
res = sqlQuery2(con, statement, info="USER-Class below USDA-silt not found.")
class_below_usda_silt_limit=res[1,"upper_limit"] #get uppper limit of respective class
statement=paste0("SELECT soil_id, sum(fraction) AS a_cum_below INTO t_cum_below",
" FROM r_soil_contains_particles WHERE class_id<=",class_below_usda_silt,
" GROUP BY soil_id") #statement-statement that produces a table containing the cumulative fractions up to class_below_usda_silt for each soil
}
res = sqlQuery2(con, statement, info="cumulative fractions up to class_below_usda_silt") #produce a table containing the for each soil
# print("computing USDA-silt-content...")
# statement = "UPDATE (soils LEFT JOIN t_cum_above ON soils.pid=t_cum_above.soil_id) LEFT JOIN t_cum_below ON soils.pid=t_cum_below.soil_id".
# " SET a_silt = ((a_cum_above-a_cum_below)*(silt_upper_limit-class_below_usda_silt_limit))/(class_above_usda_silt_limit-class_below_usda_silt_limit)"
# res = sqlQuery2(con, statement, info="")
# if(!res)
# die("Could not compute USDA-silt content (sql_err_msg).")
#
print("OK")
#### compute sand fraction according to USDA ####
statement = paste0("select class_id, upper_limit from particle_classes where upper_limit>=",sand_upper_limit," order by upper_limit")
#get ID of class that is at or just over usda sand
res = sqlQuery2(con, statement, info="")
if (nrow(res) ==0)
stop("All specified particle classes are finer than USDA-sand, cannot compute K-factor.")
class_above_usda_sand =res[1, "class_id"] #get the user-class that sits just above USDA-sand
class_above_usda_sand_limit=res[1, "upper_limit"] #get uppper limit of respective class
statement = "DROP TABLE t_cum_above"
res = sqlQuery(con, query = statement, errors = FALSE) #delete any existing table
statement = paste0("SELECT soil_id, sum(fraction) AS a_cum_above INTO t_cum_above",
" FROM r_soil_contains_particles WHERE class_id<=", class_above_usda_sand," GROUP BY soil_id")
res = sqlQuery2(con, statement, info="cumulative fractions up to class_above_usda_sand") #produce a table containing the cumulative fractions up to class_above_usda_sand for each soil
statement = "DROP TABLE t_cum_below"
res = sqlQuery(con, query = statement, errors = FALSE) #delete any existing table
class_below_usda_sand=class_above_usda_sand-1 #the class below USDA sand
if (class_below_usda_sand==0) #no lower classes that can be used as a point for interpolation
{
class_below_usda_sand_limit=0 #interpolation starts at 0
statement=paste0("SELECT soil_id, 0 AS a_cum_below INTO t_cum_below",
" FROM r_soil_contains_particles GROUP BY soil_id") #statement-statement that produces a table containing zeros for each soil
} else #use the nearest lower class as a point for interpolation
{
statement = paste0("select upper_limit from particle_classes where class_id<=", class_below_usda_sand)
res = sqlQuery2(con, statement, info="compute sand")
if (nrow(res)==0)
stop("USER-Class below USDA-sand not found.")
class_below_usda_sand_limit=res[1, "upper_limit"] #get uppper limit of respective class
statement=paste0("SELECT soil_id, sum(fraction) AS a_cum_below INTO t_cum_below",
" FROM r_soil_contains_particles WHERE class_id<=", class_below_usda_sand,
" GROUP BY soil_id") #statement-statement that produces a table containing the cumulative fractions up to class_below_usda_sand for each soil
}
res = sqlQuery2(con, statement, info="fractions up to class_below_usda_sand") #produce a table containing the cumulative fractions up to class_below_usda_sand for each soil
print("computing USDA-sand-content...")
statement = paste0("UPDATE (soils LEFT JOIN t_cum_above ON soils.pid=t_cum_above.soil_id) LEFT JOIN t_cum_below ON soils.pid=t_cum_below.soil_id",
" SET a_sand = ((a_cum_above-a_cum_below)*(",sand_upper_limit-class_below_usda_sand_limit, "))/(",class_above_usda_sand_limit-class_below_usda_sand_limit,")")
res = sqlQuery2(con, statement, info="compute USDA-sand")
print("OK")
statement = paste0("UPDATE (soils LEFT JOIN t_cum_above ON soils.pid=t_cum_above.soil_id) LEFT JOIN t_cum_below ON soils.pid=t_cum_below.soil_id",
" SET a_silt=1-a_sand-a_clay")
res = sqlQuery2(con, statement, info="set USDA-silt")
#### compute subfactors ####
print("computing f_cl_si factor...")
statement = "UPDATE soils SET a_f_cl_si=1 WHERE a_silt=0" #formula doesn't work for zero silt content - factor is set to 1
res = sqlQuery(con, statement, errors = FALSE) #raises an error, when no row affected, so suppress error message
statement = "UPDATE soils SET a_f_cl_si = exp(0.3*log(a_silt/(a_clay+a_silt))) WHERE a_silt>0"
res = sqlQuery(con, statement, errors=FALSE)
print("OK")
print("computing f_csand factor...")
#statement = "UPDATE soils SET a_f_csand = (0.2+(0.3*exp(-0.256*a_sand*100*(1-a_silt))))" #as cited in SWAT manual
statement = "UPDATE soils SET a_f_csand = (0.2+(0.3*exp(-0.0256*a_sand*100*(1-a_silt))))" #as in Williams, 1995
res = sqlQuery2(con, statement, info="computing f_csand factor")
print("OK")
print("computing f_hisand factor...")
statement = "UPDATE soils SET a_f_hisand = (1-(0.7*(1-a_sand)/((1-a_sand)+exp(-5.51+22.9*(1-a_sand)))))"
res = sqlQuery2(con, statement, info="f_hisand factor")
print("OK")
print("computing f_orgc factor...")
statement = "UPDATE soils SET a_f_orgc = (1-(0.25*b_om/1.72*100/(b_om/1.72*100+exp(3.72-2.95*b_om/1.72*100))))"
res = sqlQuery2(con, statement, info="computing f_orgc factor")
print("OK")
#### compute K and insert results into table soil
print("computing MUSLE-K and inserting into <soils>...")
statement = "UPDATE soils SET a_musle_k = a_f_csand*a_f_cl_si*a_f_orgc*a_f_hisand"
res = sqlQuery2(con, statement, info="computing MUSLE-K")
statement = "DROP TABLE t_cum_above"
res = sqlQuery(con, query = statement, errors = FALSE) #delete any existing table
statement = "DROP TABLE t_cum_below"
res = sqlQuery(con, query = statement, errors = FALSE) #delete any existing table
} #end computation musle k
#copy from other tables ####
if ("MUSLE-K" %in% toupper(copy_from_other_tables))
{
statement = "UPDATE soil_veg_components INNER JOIN soils ON soils.pid=soil_veg_components.soil_id SET musle_k = a_musle_k;"
res = sqlQuery2(con, statement, info="copy MUSLE-K")
}
if ("MANNING-N" %in% toupper(copy_from_other_tables))
{
statement = "UPDATE soil_veg_components INNER JOIN vegetation ON vegetation.pid=soil_veg_components.veg_id SET manning_n = c_manning_n;"
res = sqlQuery2(con, statement, info="copy Manning-n")
}
if ("MUSLE-C" %in% toupper(copy_from_other_tables))
{
statement = "UPDATE soil_veg_components INNER JOIN vegetation ON vegetation.pid=soil_veg_components.veg_id SET musle_c = c_musle_c;"
res = sqlQuery2(con, statement, info="copy MUSLE-C")
}
if (is.finite(setP))
{
if (setP < 0 | setP>1)
warning("P-factor (setP) should be between 0 and 1, ignored.")
else
{
statement = paste0("UPDATE soil_veg_components set musle_p = ", setP,";")
res = sqlQuery2(con, statement, info="copy MUSLE-C")
}
}
# update table meta_info
write_metainfo(con,
"db_prepare_musle()",
c("soils","soil_veg_components"), "soils.a_*, soil_veg_components.musle_k",
"Computation of MUSLE K-factor ",
FALSE)
if(verbose) message("% OK")
if(verbose) message("%")
if(verbose) message("% All data written successfully. Close ODBC connection.")
tryCatch(odbcClose(con), error=function(e){})
if(verbose) message("%")
if(verbose) message("% DONE!")
if(verbose) message("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
} # EOF
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.