R/lump_grass_post.R
In tpilz/lumpR: Landscape Unit Mapping Program for R
Defines functions
flow_velocity
channel_depth
channel_width
channel_slope
channel_length
sub_route
centroid
lump_grass_post
Documented in
lump_grass_post
# lumpR/lump_grass_post.R
# Copyright (C) 2014-2018 Tobias Pilz, 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/>.
#' Post-processing of Landscape Unit deviation using GRASS GIS
#'
#' Creates raster map of Landscape Units and files containing information of
#' and parameter estimation for Subbasins and Landscape Units in the catchment
#' using outputs of \code{\link[lumpR]{lump_grass_prep}} and \code{\link[lumpR]{prof_class}}.
#'
#' @param mask Name of mask raster map masking the study area. E.g. output \code{basin_out}
#' of \code{\link[lumpR]{calc_subbas}}.
#' @param dem DEM raster map in GRASS location as used in \code{\link[lumpR]{lump_grass_prep}}.
#' @param subbasin Subbasin raster map in GRASS location as used in \code{\link[lumpR]{lump_grass_prep}}
#' and/or created by \code{\link[lumpR]{calc_subbas}}.
#' @param recl_lu Name of GRASS reclassification file: EHA -> LU. Output of
#' \code{\link[lumpR]{prof_class}}. If omitted, the existing raster layer \code{lu} is used.
#' @param lu Input or Output: Name of Landscape Units (LU) raster map already existing in
#' GRASS location or to be generated using \code{recl_lu}.
#' @param eha Name of Elementary Hillslope Areas (EHA) raster map in GRASS
#' location. Output of \code{\link[lumpR]{lump_grass_prep}}.
#' @param fill_holes TRUE: fill any holes in map \code{eha} (e.g. as result of skipped/non-classified
#' EHAs) by growing
#' @param flowacc Name of flow accumulation raster map in GRASS location. Can
#' be created with \code{\link[lumpR]{lump_grass_prep}}.
#' @param flowdir Name of flow direction raster map in GRASS location. Can
#' be created with \code{\link[lumpR]{lump_grass_prep}}.
#' @param stream_horton Name of Horton stream order raster map in GRASS location. Can
#' be created with \code{\link[lumpR]{lump_grass_prep}}. If left empty, the channel length,
#' slope and retention times are set to \code{NA}.
#' @param soil_depth Name of soil depth [mm] raster map in GRASS location. If \code{NULL}
#' (default), \code{na_val} is used.
#' @param sdr Name of sediment delivery ratio [-] raster map in GRASS location. If empty,
#' this optional column is omitted.
#' @param dir_out Character string specifying output directory (will be created;
#' any overwriting will be prompted).
#' @param sub_ofile Output: Name of subbasin statistics file containing subbasin
#' parameters. See \code{Details} below. If \code{NULL} (default) this file
#' will not be generated.
#' @param lu_ofile Output: Name of file containing subbasins and the corresponding
#' LUs with their fraction of area in the subbasin. If \code{NULL} (default)
#' this file will not be generated.
#' @param lupar_ofile Output: Name of file containing LUs and related parameters.
#' See \code{Details} below. If \code{NULL} (default) this file will not be generated.
#' @param groundwater Flag: 1: respect groundwater and infer parameters. 0 (default):
#' Ignore groundwater and associated parameters.
#' @param na_val Value used to indicate \code{NA} values in output files. Default: \code{NA}.
#' @param keep_temp \code{logical}. Set to \code{TRUE} if temporary files shall be kept
#' in the GRASS location, e.g. for debugging or further analyses. Default: \code{FALSE}.
#' @param overwrite \code{logical}. Shall output of previous calls of this function be
#' deleted? If \code{FALSE} the function returns an error if output already exists.
#' Default: \code{FALSE}.
#' @param silent \code{logical}. Shall the function be silent (also suppressing warnings
#' of internally used GRASS functions)? Default: \code{FALSE}.
#'
#' @return Function returns nothing. Output files (\code{sub_ofile, lupar_ofile, lu_ofile})
#' are written into output directory and raster map (\code{lu}) exported into GRASS
#' location as specified in arguments.
#'
#' @note Prepare GRASS location and necessary raster files in advance (e.g. using
#' \code{\link[lumpR]{lump_grass_prep}}) and start GRASS session in R using
#' \code{\link[rgrass7]{initGRASS}}.
#'
#' \bold{IMPORTANT:} Herein, when specifying the GRASS input maps, please do
#' explicitly refer to the mapset if it is different from the mapset given in
#' initGRASS() (even PERMANENT!), as otherwise internally used readRAST() command
#' resulted in errors under Windows.
#'
#' TODO:\cr
#' - check empirical formulas for channel width and channel depth\cr
#' - LU parameter estimation\cr
#' - include options to add parameters manually in case data are available\cr
#' - include option to make function more efficient regarding RAM usage (e.g. by
#' wrting/reading temporary raster data to/from disk) at the cost of higher computational burden
#'
#'
#' @details
#' \bold{Subbasin parameters}\cr
#' Subbasin parameter estimation given in \code{sub_ofile} contains:
#'
#' \emph{pid}\cr
#' Subbasin identifier.
#'
#' \emph{description}\cr
#' Subbasin description you can fill in manually if needed.
#'
#' \emph{lat}\cr
#' Latitude of subbasin centroid in \emph{decimal degrees} (negative values for southern hemisphere).
#'
#' \emph{lon}\cr
#' Longitude of subbasin centroid in \emph{decimal degrees west of Greenwhich}, e.g.
#' Greenwich: 0°, New York: 75°, Berlin: 345°.
#'
#' \emph{elev}\cr
#' Average elevation above sea level of subbasin \emph{m}.
#'
#' \emph{area}\cr
#' Subbasin area in \emph{km^2}.
#'
#' \emph{x}\cr
#' X-coordinate of subbasin centroid in units of GRASS location.
#'
#' \emph{y}\cr
#' Y-coordinate of subbasin centroid in units of GRASS location.
#'
#' \emph{drains_to}\cr
#' Pid of subbasin the current subbasin drains to. It is determined by identifying
#' the cell with the highest flow accumulation value (= subbasin outlet). By
#' deriving the drainage direction value of that raster cell the next downstream
#' cell is identified along with corresponding subbasin number.
#'
#' \emph{lag_time}\cr
#' Time in \emph{days} a runoff signal in the current subbasin needs to be directed
#' from the subbasin input to the outlet. Estimated from channel geometry (see below).
#'
#' \emph{retention}\cr
#' Maximum time period in \emph{days} over which a runoff signal is distributed by
#' the routing process. Estimated from channel geometry (see below).
#'
#' \emph{chan_len}\cr
#' Length of the main channel of the respective subbasin in \emph{m}. Estimated
#' from channel geometry (see below).
#'
#' \emph{channel geometry}\cr
#' Main channel length: For each subbasin the main channel is determined from Horton
#' stream raster map. Its length is then calculated depending on raster resolution
#' and flow direction.\cr
#' Channel slope: Minimum (= inflow) and maximum (= outflow) flow accumulation and
#' elevation of the corresponding raster cells are determined. The difference in
#' elevation is divided by main channel length.\cr
#' Channel width: Maximum flow accumulation is determined and the corresponding
#' drainage area calculated by resolution of raster cells. Channel width is then
#' calculated from the empirical formula: width[m] = 1.29 * darea[km2] ^ 0.6.\cr
#' Channel depth: Empirical formula based on drainage area: depth[m] = 0.13 * darea[km2] ^ 0.4.\cr
#' Flow velocity: Is calculated using Mannings equation with an n-value of 0.075.
#' Flow velocities are calculated to derive flow travel times (velocity / channel
#' length) for bankfull (= high flow condition), 2/3 (= average conditions) and
#' 1/10 (= low flow conditions) water levels to derive lag time (travel time for
#' average conditions) and retention time (max - min travel time).
#'
#'
#' \bold{Landscape Units in Subbasins}\cr
#' \code{lu_ofile} contains:
#'
#' \emph{subbas_id}\cr
#' Subbasin identifier.
#'
#' \emph{lu_id}\cr
#' Landscape Unit identifier.
#'
#' \emph{fraction}\cr
#' Areal fraction of Landscape Unit within corresponding Subbasin.
#'
#'
#' \bold{Landscape Unit parameters}\cr
#' Landscape Unit parameter estimation given in \code{lupar_ofile} contains:
#'
#' \emph{pid}\cr
#' Landscape Unit identifier.
#'
#' \emph{description}\cr
#' Description for this Landscape Unit. Can be adjusted manually if you want.
#' Generally set to \code{na_val}.
#'
#' \emph{kf_bedrock}\cr
#' Hydraulic conductivity of bedrock. Fill in values manually. Generally set
#' to \code{na_val}. Use in WASA model is optional (see WASA documentation
#' -> 'kfsu' and notes on 'bedrock').
#'
#' \emph{slopelength}\cr
#' Slope length of Landscape Unit [m]. Value can be obtained from output of
#' \code{\link[lumpR]{prof_class}}. Herein set to \code{na_val}.
#'
#' \emph{soil_depth}\cr
#' Soil depth in \emph{mm} averaged over respective landscape unit. Use in
#' WASA model is optional (see WASA documentation -> 'meandep' and notes on
#' 'bedrock'). Herein set to \code{na_val}.
#'
#' \emph{allu_depth}\cr
#' Depth of alluvial soils in \emph{mm}. Use in WASA model is optional (see
#' WASA documentation -> 'maxdep' and notes on 'bedrock'). Herein set to \code{na_val}.
#'
#' \emph{riverbed_depth}\cr
#' Depth of river bed below Terrain Component in \emph{mm}. Use in WASA model
#' is optional (see WASA documentation -> 'riverbed' and notes on 'bedrock').
#' Herein set to \code{na_val}.
#'
#' \emph{gw_flag}\cr
#' Groundwater flag: 0: no groundwater in this LU. 1: LU contains groundwater.
#' At the moment set to 0 by default for every LU. Use in WASA model is optional
#' (see WASA documentation -> 'gw_flag' and notes on groundwater).
#'
#' \emph{gw_dist}\cr
#' Initial depth of groundwater below surface in \emph{mm}. At the moment set to
#' 1000 by default for every LU if \code{groundwater = 1}. Use in WASA model is optional
#' (see WASA documentation -> 'gw_dist' and notes on groundwater).
#'
#' \emph{frgw_delay}\cr
#' Storage coefficient for groundwater outflow in \emph{days}. At the moment set to
#' 200 by default for every LU if \code{groundwater = 1}. Use in WASA model is optional
#' (see WASA documentation -> 'frgw_delay' and notes on groundwater).
#'
#' \emph{sdr_lu} (optional)\cr
#' sediment delivery ratio from raster map \code{sdr}, if specified
#'
#'
#' @references Source code based on \code{SHELL} and \code{MATLAB} scripts of Till Francke.
#'
#' lumpR package introduction with literature study and sensitivity analysis:\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
#'
#' Theory of LUMP:\cr
#' Francke, T.; Guentner, A.; Mamede, G.; Mueller, E. N. and Bronstert, A (2008):
#' Automated catena-based discretization of landscapes for the derivation of
#' hydrological modelling units. \emph{International Journal of Geographical
#' Information Science, Informa UK Limited}, 22(2), 111-132, DOI: 10.1080/13658810701300873
#'
#' Subbasin Parameters:\cr
#' Bronstert, A., Guentner, A., Jaeger, A., Krol, M. & Krywkow, J. (1999): Grossraeumige
#' hydrologische Parametrisierung und Modellierung als Teil der integrierten
#' Modellierung. In: Fohrer, N. & Doell, P. (Eds.): Modellierung des Wasser- und
#' Stofftransports in grossen Einzugsgebieten. \emph{Kassel University Press}, Kassel,
#' Germany, 31-40.
#'
#' Guentner, A. (2002): Large-scale hydrological modelling in the semi-arid
#' North-East of Brazil. \emph{PIK Report 77}, Potsdam Institute for Climate
#' Impact Research, Potsdam, Germany.
#'
#'
#' @author Tobias Pilz \email{tpilz@@uni-potsdam.de}
lump_grass_post <- function(
### INPUT ###
mask=NULL,
dem=NULL,
recl_lu=NULL,
lu=NULL,
subbasin=NULL,
eha=NULL,
flowacc=NULL,
flowdir=NULL,
stream_horton=NULL,
soil_depth=NULL,
sdr=NULL,
### OUTPUT ###
dir_out="./",
sub_ofile=NULL,
lu_ofile=NULL,
lupar_ofile=NULL,
### OPTIONS ###
fill_holes=T,
groundwater=0,
na_val=NA,
keep_temp=F,
overwrite=F,
silent=F
) {
### PREPROCESSING ###----------------------------------------------------------
if(!silent) message("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
if(!silent) message("% START lump_grass_post()")
if(!silent) message("%")
if(!silent) message("% Initialise function...")
# checks #---------------------------------------------------------------------
tryCatch(gmeta(), error = function(e) stop("Cannot execute GRASS commands. Maybe you forgot to run initGRASS()?"))
if(is.null(mask))
stop("The name of a raster within the mapset of your initialised GRASS session to be used as catchment MASK in GRASS has to be given!")
if(!is.character(mask))
stop("The name of a raster to be used as catchment MASK in GRASS must be a string (argument 'mask')!")
if(is.null(dem))
stop("The name of a DEM within the mapset of your initialised GRASS session has to be given!")
if(is.null(subbasin))
stop("A name for needed subbasin raster map within the mapset of your initialised GRASS session has to be given!")
if(is.null(eha))
stop("A name for needed EHA raster map within the mapset of your initialised GRASS session has to be given!")
if(is.null(flowdir))
stop("A name for needed flow direction raster map within the mapset of your initialised GRASS session has to be given!")
if(is.null(flowacc))
stop("A name for needed flow accumulation raster map within the mapset of your initialised GRASS session has to be given!")
if(is.null(lu))
stop("A name for the landscape units raster map to be generated (recl_lu = NULL) or already within the mapset of your initialised GRASS session has to be given")
if(!is.null(recl_lu) && !file.exists(recl_lu))
stop(paste0("Could not find file '",recl_lu,"' (argument 'recl_lu')"))
#check existence of raster maps
check_raster(mask,"mask")
check_raster(dem,"dem")
if(is.null(recl_lu)) check_raster(lu,"lu")
check_raster(subbasin,"subbasin")
check_raster(eha,"eha")
check_raster(flowacc,"flowacc")
check_raster(flowdir,"flowdir")
check_raster(stream_horton,"stream_horton")
if(!is.null(soil_depth)) check_raster(soil_depth,"soil_depth")
if(!is.null(sdr)) check_raster(sdr,"sdr")
# suppress annoying GRASS outputs
tmp_file <- file(tempfile(), open="wt")
sink(tmp_file, type="output")
# also supress warnings in silent mode
if(silent){
tmp_file2 <- file(tempfile(), open="wt")
sink(tmp_file2, type="message")
oldw <- getOption("warn")
options(warn = -1)
}
# create output dir
dir.create(dir_out, recursive=T, showWarnings=F)
# check output directory
if (!overwrite & (file.exists(paste(dir_out,sub_ofile,sep="/")) |
file.exists(paste(dir_out,lu_ofile,sep="/")) |
file.exists(paste(dir_out,lupar_ofile,sep="/"))) )
stop(paste0("Output file(s) ", sub_ofile, ", ",lu_ofile, ", and/or ", lupar_ofile, " already exist(s) in ", dir_out, "!"))
# remove mask if there is any (and ignore error in case there is no mask)
# set mask to make sure calculations are done exactly within expected area
cmd_out <- execGRASS("g.copy", raster=paste0(mask,",MASK"), flags = "overwrite", intern = T)
# remove output of previous function calls if overwrite=T
if (overwrite) {
cmd_out <- execGRASS("g.remove", type="raster", pattern=paste("*_t,", lu, sep=","), flags=c("f", "b"), intern=T)
} else {
# remove temporary maps in any case
cmd_out <- execGRASS("g.remove", type="raster", pattern="*_t", flags=c("f", "b"), intern=T)
}
if(!silent) message("% OK")
### CALCULATIONS ###-----------------------------------------------------------
# prepare output directory, read GRASS data #-----------------------------
tryCatch({
if(!silent) message("%")
if(!silent) message("% Prepare output and import GRASS data...")
# load rasters into R
dem_rast <- read_raster(dem)
accum_rast <- read_raster(flowacc)
dir_rast <- read_raster(flowdir)
cur_mapset = execGRASS("g.mapset", flags="p", intern=TRUE) #determine name of current mapset
sub_rast <- read_raster(subbasin)
horton_rast <- read_raster(stream_horton)
# ... raster values as matrix
dem_mat <- getValues(dem_rast, format="matrix")
rm(dem_rast)
sub_mat <- getValues(sub_rast, format="matrix")
coords <- xyFromCell(sub_rast, 1:length(sub_rast))
rm(sub_rast)
accum_mat <- getValues(accum_rast, format="matrix")
rm(accum_rast)
dir_mat <- getValues(dir_rast, format="matrix")
rm(dir_rast)
horton_mat <- getValues(horton_rast, format="matrix")
rm(horton_rast)
gc(); gc()
# resolution
g_meta <- gmeta()
resol <- mean(g_meta$nsres, g_meta$ewres)
if(!silent) message("% OK")
}, error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
stop(paste(e))
})
# generate lu raster #---------------------------------------------------------
tryCatch({
if(!silent) message("%")
if (is.null(recl_lu)) #use existing LU-map
{
if(!silent) message(paste0("% NOTE: Use existing raster map ", lu, " as specified."))
} else {
if(!silent) message(paste0("% Generate landscape units raster map '", lu, "'..."))
#reclass EHA according to reclass file generated by prof_class to get LU
cmd_out=execGRASS("r.reclass", input=eha, output=lu, rules=recl_lu, intern = TRUE)
if(!is.null(attr(cmd_out, "status")))
stop(paste0("Could not reclass '",eha,"' into '",lu,"'"))
# growing radius (parameter for r.grow)
GROWRAD <- 20
if (fill_holes)
while (TRUE)
{
# look for empty patches
cmd_out=execGRASS("r.mapcalc", expression=paste0("grow_eval_t = if(isnull(", lu, "), 1, 0)"), flags = "overwrite", intern=TRUE)
na_eval <- execGRASS("r.stats", input=paste0("grow_eval_t"), flags=c("n", "p"), intern=TRUE, ignore.stderr = T)
if (length(na_eval)==0 ) break #no empty patches found
# if (grepl(pattern="[0-9]+.*[\b]+",x=tail(na_eval, n=1)))
# na_eval = na_eval[-length(na_eval)] #last line contains progress indicator, remove
na_line=grepl(na_eval, pattern = "^1 .*\\%$")
if (!any(na_line)) break #no empty patches found
na_res <- strsplit(x=na_eval[na_line], split=" +")
if(!silent) message(paste0("% -> ", tail(unlist(na_res),1), " NAs in LU-map '", lu,"'"))
# otherwise grow LU map to fill gaps
if(!silent) message(paste0("% -> Filling empty patches in ",lu,", iterating..."))
cmd_out=execGRASS("r.grow", input=lu, output="lu_t", radius=GROWRAD, flags="overwrite", intern=TRUE)
if(!is.null(attr(cmd_out, "status")))
stop(paste0("Could not grow '",lu,"'"))
# r.grow converts type CELL to type DCELL; convert back to CELL
cmd_out <- execGRASS("r.mapcalc", expression=paste0(lu," = int(lu_t)"), flags = "overwrite", intern = T)
cmd_out <- execGRASS("g.remove", type="raster", pattern="lu_t", flags=c("f", "b"), intern=T)
# x=execGRASS("r.mapcalc", expression = paste0(lu, " = ", mask, " * tt_eha_grown"), flags="overwrite", intern=TRUE)
#
# x=execGRASS("g.remove", type="raster", name="tt_eha_grown", intern=TRUE)
} # while true (grow-loop)
cmd_out=execGRASS("g.remove", type="raster", name="grow_eval_t", intern=TRUE)
if(!silent) message("% OK")
} # recl_lu was given
},
error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
# remove output so far
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern=paste0("grow_eval_t,", lu), flags=c("f", "b"), intern=TRUE)
stop(paste(e))
}
)
# subbasin statistics #--------------------------------------------------------
tryCatch({
if(!is.null(sub_ofile)) {
if(!silent) message("%")
if(!silent) message("% Calculate subbasin statistics...")
sub_stats <- execGRASS("r.stats", input=subbasin, flags=c("a", "n"), intern=TRUE, ignore.stderr = T)
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(sub_stats, n=1)))
sub_stats <- sub_stats[-length(sub_stats)] #last line contains progress indicator, remove
sub_stats <- matrix(as.numeric(gsub("%", "", unlist(strsplit(sub_stats, split=" +")))), ncol=2, byrow=T,
dimnames=list(NULL,c("pid", "area")))
# convert m2 to km2
sub_stats[,"area"] <- round(sub_stats[,"area"]/1e6, 2)
# calculate stats of LUs in each subbasin and subbasin drainage ("drains_to")
sub_stats <- cbind(sub_stats, na_val, na_val, na_val, na_val, na_val, na_val, na_val, na_val, na_val, na_val, na_val)
colnames(sub_stats)[c(3:13)] <- c("x", "y", "lon", "lat", "elev", "drains_to", "lag_time", "retention", "description", "a_stream_order", "chan_len")
sub_lu_stats <- NULL
s_row <- 0
for (SUB in sub_stats[,1]) {
s_row <- s_row +1
# bounding box of current subbasin
sub_arrind <- which(sub_mat == SUB, arr.ind = T)
rmin <- min(sub_arrind[,"row"])
rmax <- max(sub_arrind[,"row"])
cmin <- min(sub_arrind[,"col"])
cmax <- max(sub_arrind[,"col"])
# crop raster matrices
dem_crop <- dem_mat[rmin:rmax, cmin:cmax]
horton_crop <- horton_mat[rmin:rmax, cmin:cmax]
sub_crop <- sub_mat[rmin:rmax, cmin:cmax]
dir_crop <- dir_mat[rmin:rmax, cmin:cmax]
accum_crop <- accum_mat[rmin:rmax, cmin:cmax]
# get coordinates of cells
coords_crop <- coords[which(t(sub_mat) == SUB),]
# set all values not overlapping with the current subbasin to NA
na_cells <- which(sub_crop != SUB)
sub_crop[na_cells] <- NA
dem_crop[na_cells] <- NA
horton_crop[na_cells] <- NA
dir_crop[na_cells] <- NA
accum_crop[na_cells] <- NA
# COORDINATES OF SUBBASIN centroids in GRASS units #
sub_centr <- centroid(coords_crop)
sub_stats[s_row, "x"] <- round(sub_centr[,"x"])
sub_stats[s_row, "y"] <- round(sub_centr[,"y"])
# convert to decimal degree
sub_centr <- data.frame(x = sub_centr[,"x"], y = sub_centr[,"y"])
coordinates(sub_centr) <- c("x", "y")
proj4string(sub_centr) <- getLocationProj()
longlat <- spTransform(sub_centr, "+proj=longlat")
sub_stats[s_row, "lon"] <- round(coordinates(longlat)[,"x"], 4)
sub_stats[s_row, "lat"] <- round(coordinates(longlat)[,"y"], 4)
# longitude as decimal degrees west of Greenwich
if(sub_stats[s_row, "lon"] < 0)
sub_stats[s_row, "lon"] <- -1*sub_stats[s_row, "lon"]
else
sub_stats[s_row, "lon"] <- 360 - sub_stats[s_row, "lon"]
# AVERAGE SUBBASIN ELEVATION #
sub_stats[s_row, "elev"] <- round(mean(dem_crop, na.rm=T))
# SUBBASIN drainage #
sub_stats[s_row,"drains_to"] <- sub_route(SUB,sub_mat,accum_mat,dir_mat) # internal function, see below
# SUBBASIN PARAMETERS #
# calc main channel length
chan_len <- channel_length(horton_crop, dir_crop, resol)
# calc main channel average slope
chan_slope <- channel_slope(horton_crop, accum_crop, dem_crop, chan_len)
# calc main channel width
chan_width <- channel_width(accum_crop, resol)
# calc main channel depth
chan_depth <- channel_depth(accum_crop, resol)
# calc flow velocites for bankful, 2/3 and 1/10 filling
# Manning's n = 0.075 (very weedy reaches, deep pools, or floodways with heavy stand of timber and underbrush)
# http://www.fsl.orst.edu/geowater/FX3/help/8_Hydraulic_Reference/Mannings_n_Tables.htm
chan_velo_full <- flow_velocity(chan_width, chan_depth, chan_slope, n=0.075)
chan_velo_med <- flow_velocity(chan_width, 2/3 * chan_depth, chan_slope, n=0.075)
chan_velo_low <- flow_velocity(chan_width, 1/10 * chan_depth, chan_slope, n=0.075)
# calc travel times (translation and retention) [d]
# approach according to Bronstert et al. (1999), Guentner (2002)
flowtime_min <- chan_len / chan_velo_full / 86400
flowtime_med <- chan_len / chan_velo_med / 86400
flowtime_max <- chan_len / chan_velo_low / 86400
# check results
retention <- flowtime_max - flowtime_min
if(!is.finite(flowtime_med) || !is.finite(retention))
warning(paste0("Could not calculate finite subbasin parameters for subbasin ", SUB))
# save
sub_stats[s_row, "lag_time"] <- round(flowtime_med, 3)
sub_stats[s_row, "retention"] <- round(retention, 3)
sub_stats[s_row, "chan_len"] <- round(chan_len)
} # subbasin loop
# write output file
write.table(sub_stats, paste(dir_out, sub_ofile, sep="/"), quote=F, row.names=F, sep="\t")
if(!silent) message("% OK")
} # sub_ofile given?
}, error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
# remove output so far
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern=paste0("*_t,", lu), flags=c("f", "b"), intern=TRUE)
stop(paste(e))
})
# lu statistics #--------------------------------------------------------------
tryCatch({
if(!is.null(lu_ofile)) {
if(!silent) message("%")
if(!silent) message("% Calculate LU statistics...")
sub_lu_stats_t <- execGRASS("r.stats", input=paste0(subbasin,",",lu), flags=c("n", "c"), intern=TRUE, ignore.stderr = T)
if(!is.null(attr(sub_lu_stats_t, "status")))
stop(paste0("Could not access lu-map '",lu,"'."))
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(sub_lu_stats_t, n=1))) #check if last line contains progress indicator, remove
sub_lu_stats_t = sub_lu_stats_t[-length(sub_lu_stats_t)]
sub_lu_stats_t = sub_lu_stats_t[grepl(x=sub_lu_stats_t, pattern = "^[0-9]")] #remove non-number lines
sub_lu_stats_t = as.numeric(gsub("%", "", unlist(strsplit(sub_lu_stats_t, split=" +")))) #convert to numbers and remove %
sub_lu_stats_t2 =data.frame(matrix(sub_lu_stats_t, ncol=3, byrow=T))
names(sub_lu_stats_t2)=c("subbas_id","lu_id","cells")
subbasin_cells=aggregate(x = sub_lu_stats_t2$cells, by = list(subbas_id=sub_lu_stats_t2$subbas_id), FUN=sum)
sub_lu_stats_t2=merge(sub_lu_stats_t2,subbasin_cells)
sub_lu_stats_t2$fraction=sub_lu_stats_t2$cells / sub_lu_stats_t2$x #compute fraction
write.table(sub_lu_stats_t2[,c("subbas_id", "lu_id", "fraction")], paste(dir_out, lu_ofile, sep="/"), quote=F, row.names=F, sep="\t")
if(!silent) message("% OK")
} # lu_ofile given?
}, error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
# remove output so far
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern=paste0("*_t,", lu), flags=c("f", "b"), intern=TRUE)
stop(paste(e))
})
# lu parameters #--------------------------------------------------------------
tryCatch({
if(!is.null(lupar_ofile)) {
if(!silent) message("%")
if(!silent) message("% Calculate LU parameters")
# identify LUs
lu_ids <- execGRASS("r.stats", input=lu, flags=c("n"), intern=TRUE, ignore.stderr = T)
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(lu_ids, n=1)))
lu_ids <- as.numeric(lu_ids[-length(lu_ids)]) #last line contains progress indicator, remove
# initialise LU parameter output object
lu_par <- matrix(na_val, ncol=11, nrow=length(lu_ids),
dimnames=list(NULL, c("pid", "description", "kf_bedrock", "slopelength",
"soil_depth", "allu_depth", "riverbed_depth",
"gw_flag", "gw_dist", "frgw_delay", "sdr_lu")))
# calculate mean soil depth for every LU
if (!is.null(soil_depth) && soil_depth!="")
{
cmd_out <- execGRASS("r.univar", zones=lu, map=soil_depth, fs=",", flags=c("t"),intern=T, ignore.stderr = T)
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(cmd_out, n=1)))
cmd_out <- cmd_out[-length(cmd_out)] #last line contains progress indicator, remove
cmd_out <- strsplit(cmd_out, ",")
cmd_out <- matrix(unlist(cmd_out[-1]), ncol=length(cmd_out[[1]]), byrow=T,
dimnames=list(NULL, cmd_out[[1]]))
lu_depth <- as.numeric(cmd_out[,"mean"])
lu_ids = as.numeric(cmd_out[,"zone"])
# quick check of soil depths
if(any(lu_depth > 10000))
warning("There are average LU soil depths of more than 10 m which is a large (but not impossible) value. Check your input data and units ('soil_depth' in [mm])!")
if(any(lu_depth < 100))
warning("There are average LU soil depths of less than 10 cm which is a small (but not impossible) value. Check your input data and units ('soil_depth' in [mm])!")
} else lu_depth=na_val
lu_par[,"pid"] <- lu_ids
lu_par[,"soil_depth"] <- lu_depth
if (!is.null(sdr) && sdr!="")
{
cmd_out <- execGRASS("r.univar", zones=lu, map=sdr, fs=",", flags=c("t"),intern=T, ignore.stderr = T)
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(cmd_out, n=1)))
cmd_out <- cmd_out[-length(cmd_out)] #last line contains progress indicator, remove
cmd_out <- strsplit(cmd_out, ",")
cmd_out <- matrix(unlist(cmd_out[-1]), ncol=length(cmd_out[[1]]), byrow=T,
dimnames=list(NULL, cmd_out[[1]]))
sdr_vals <- as.numeric(cmd_out[,"mean"])
lu_ids = as.numeric(cmd_out[,"zone"])
sdr_lu <- sdr_vals[match(lu_par[,"pid"], lu_ids)]
lu_par[,"sdr_lu"] <- sdr_lu
}
# groundwater parameters (so far only default values)
# TODO: alternative approaches?
# groundwater for every LU
if(groundwater) {
lu_par[,"gw_flag"] <- rep(1, length(lu_depth))
# initial depth of groundwater below surface: 1000 mm
lu_par[,"gw_dist"] <- rep(1000, length(lu_depth))
# storage coefficient for groundwater outflow [days]; TODO: estimate from baseflow analysis?!
lu_par[,"frgw_delay"] <- rep(200, length(lu_depth))
} else {
lu_par[,"gw_flag"] <- rep(0, length(lu_depth))
# initial depth of groundwater below surface: 1000 mm
lu_par[,"gw_dist"] <- rep(na_val, length(lu_depth))
# storage coefficient for groundwater outflow [days]; TODO: estimate from baseflow analysis?!
lu_par[,"frgw_delay"] <- rep(na_val, length(lu_depth))
}
# write output
write.table(lu_par, paste(dir_out, lupar_ofile, sep="/"), quote=F, row.names=F, sep="\t")
if(!silent) message("% OK")
} # lupar_ofile fiven?
}, error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
# remove output so far
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern=paste0("*_t,", lu), flags=c("f", "b"), intern=TRUE)
stop(paste(e))
})
# remove temp files
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern="*_t", flags=c("f", "b"), intern=TRUE)
if(!silent) message("%")
if(!silent) message("% DONE!")
if(!silent) message("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
} # EOF
### INTERNAL FUNCTIONS ###-----------------------------------------------------
# subbasin centroids #---------------------------------------------------------
# returns a matrix with x and y value of the centroid of the object given in coords_sub
centroid <- function(coords_sub) {
# convert so sp class
coords_sub <- as.data.frame(coords_sub)
coordinates(coords_sub) <- c("x", "y")
# calculate centroid
centr <- gCentroid(coords_sub)
centr <- centr@coords
return(centr)
}
# subbasin routing #-----------------------------------------------------------
# returns ID of downstream subbasin for the current subbasin 'sub_no'
# determined from flow accumulation and flow direction map
sub_route <- function(sub_no,sub_mat,accum_mat,dir_mat) {
# extract highest flowacc in subbasin sub_no
sub_ids <- which(sub_mat == sub_no)
accum_sub <- accum_mat[sub_ids]
accum_sub_max <- which.max(accum_sub)
# extract corresp. flowdir
dir_sub_out <- dir_mat[sub_ids[accum_sub_max]]
# determine row and col no in sub_mat
sub_rowcol <- arrayInd(sub_ids[accum_sub_max], dim(sub_mat))
# determine value of sub_mat (subbasin no.) of neighbour cell according to dir_sub_out
if (dir_sub_out < 0) { # outlet of catchment
return(9999)
} else {
if (dir_sub_out == 1) { # NE
sub_rowcol_out <- sub_rowcol + c(-1,1)
} else if (dir_sub_out == 2) { # N
sub_rowcol_out <- sub_rowcol + c(-1,0)
} else if (dir_sub_out == 3) { # NW
sub_rowcol_out <- sub_rowcol + c(-1,-1)
} else if (dir_sub_out == 4) { # W
sub_rowcol_out <- sub_rowcol + c(0,-1)
} else if (dir_sub_out == 5) { # SW
sub_rowcol_out <- sub_rowcol + c(+1,-1)
} else if (dir_sub_out == 6) { # S
sub_rowcol_out <- sub_rowcol + c(+1,0)
} else if (dir_sub_out == 7) { # SE
sub_rowcol_out <- sub_rowcol + c(+1,+1)
} else if (dir_sub_out == 8) { # E
sub_rowcol_out <- sub_rowcol + c(0,+1)
} else {
stop(paste("While determining subbasin drainage: Determined flow direction at outlet
of subbasin ", sub_no, " has value ", dir_sub_out, " but should be one of {1,2,3,4,5,6,7,8} or a negative number.", sep=""))
}
# determine subbasin the current subbasin drains to
sub_out <- sub_mat[sub_rowcol_out]
# update stats table
if (is.na(sub_out)) { # catchment outlet
return(9999)
} else {
if (sub_out == sub_no)
stop(paste0("In subbasin no. ", sub_no, " the determined downstream subbasin is this subbasin!"))
return(sub_out)
}
}
} # EOF
# main channel length #--------------------------------------------------------
# returns length of the main channel (largest value in Horton order) in [m]
# computes main stream temporary raster used in further calculations
channel_length <- function(horton, flowdir, resol) {
# determine main channel (largest value in Horton stream order)
main_chan <- na.exclude(unique(c(horton)))
# if there is no main channel (e.g. in very small reservoir subbasins) assume one cell of main stream
if (length(main_chan) == 0) {
chan_len <- resol
# accmax <- max(as.numeric(execGRASS("r.stats", input=flowacc, flags=c("n", "1"), intern=T, ignore.stderr = T))) # max flowacc
# execGRASS("r.mapcalculator", amap=flowacc, outfile="stream_main_t",
# formula=paste0("if(A == ", sprintf(accmax, fmt="%d"), ", 1, null())"))
warning(paste("Subbasin ", sub_no, " has no main channel. Assume at least one raster cell.", sep=""))
} else {
# diagonal raster cell length
dia <- sqrt(2*resol^2)
# determine raster cells of main channel
max_val <- max(main_chan)
r_main <- which(horton==max_val)
# calculate length of the main channel
chanlen_t <- flowdir[r_main] %% 2 == 0 # horizontal/vertical stream
chan_len <- length(which(chanlen_t))*resol + length(which(!chanlen_t))*dia
}
return(chan_len)
} # EOF
# main channel slope #---------------------------------------------------------
# returns average slope of main channel [m/m]
channel_slope <- function(horton, accum, dem, chan_len) {
if (is.na(chan_len) ) return(NA)
# determine main channel (largest value in Horton stream order)
main_chan <- na.exclude(unique(c(horton)))
# determine raster cells of main channel
max_val <- max(main_chan)
r_main <- which(horton==max_val)
# min and max values of flow accumulation within main channel to determin flow direction
acc_min <- which.min(accum[r_main])
acc_max <- which.max(accum[r_main])
# get dem values for respective raster cells of min and max flowacc
dem_vals <- dem[r_main][c(acc_max, acc_min)]
# compute slope [m/m]
chan_slope <- diff(dem_vals) / chan_len
# slope zero not allowed -> at least 0.00001
chan_slope <- max(chan_slope, 0.00001)
return(chan_slope)
}
# main channel width #---------------------------------------------------------
# calculate main channel width based on empirical formula: width[m] = 1.29 * darea[km2] ^ 0.6
# darea = drainage area determined from maximum flow accumulation and resolution
channel_width <- function(accum, resol) {
# determine maximum flow accumulation
max_acc <- max(accum, na.rm = T)
# calculate drainage area [km^2]
drain_area <- max_acc * resol^2 / 1e6
# calculate width according to empirical function
chan_width <- 1.29 * drain_area^(0.6)
if(chan_width > 3000)
stop(paste("Calculated channel width is ", chan_width, "m which seems unrealistic!", sep=""))
return(chan_width)
}
# main channel depth #---------------------------------------------------------
# calculate main channel width based on empirical formula: width[m] = 0.13 * darea[km2] ^ 0.4
# darea = drainage area determined from maximum flow accumulation and resolution
channel_depth <- function(accum, resol) {
# determine maximum flow accumulation
max_acc <- max(accum, na.rm = T)
# calculate drainage area [km^2]
drain_area <- max_acc * resol^2 / 1e6
# calculate depth according to empirical function
chan_depth <- 0.13 * drain_area^(0.4)
if(chan_depth > 50)
stop(paste("Calculated channel depth is ", chan_depth, "m which seems unrealistic!", sep=""))
return(chan_depth)
}
# flow velocity #--------------------------------------------------------------
# function returns flow velocity of main channel [m/s]
# calc by Gauckler-Manning-Strickler
# assume simple rectangle cross section as width >> depth for main channels
flow_velocity <- function(chan_width, chan_depth, chan_slope, n) {
# calculate cross section area (assume simple rectangle as width >> depth for main channels)
cross_area <- chan_width * chan_depth
# calculate hydraulic radius (wetted perimeter = width + 2 * depth)
hyd_radius <- cross_area / (chan_width + 2 * chan_depth)
# calculate flow velocity
chan_velo <- 1/n * hyd_radius^(2/3) * chan_slope^(1/2)
return(chan_velo)
}
tpilz/lumpR documentation built on Aug. 5, 2023, 1:31 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.
Defines functions flow_velocity channel_depth channel_width channel_slope channel_length sub_route centroid lump_grass_post
Documented in lump_grass_post
# lumpR/lump_grass_post.R
# Copyright (C) 2014-2018 Tobias Pilz, 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/>.
#' Post-processing of Landscape Unit deviation using GRASS GIS
#'
#' Creates raster map of Landscape Units and files containing information of
#' and parameter estimation for Subbasins and Landscape Units in the catchment
#' using outputs of \code{\link[lumpR]{lump_grass_prep}} and \code{\link[lumpR]{prof_class}}.
#'
#' @param mask Name of mask raster map masking the study area. E.g. output \code{basin_out}
#' of \code{\link[lumpR]{calc_subbas}}.
#' @param dem DEM raster map in GRASS location as used in \code{\link[lumpR]{lump_grass_prep}}.
#' @param subbasin Subbasin raster map in GRASS location as used in \code{\link[lumpR]{lump_grass_prep}}
#' and/or created by \code{\link[lumpR]{calc_subbas}}.
#' @param recl_lu Name of GRASS reclassification file: EHA -> LU. Output of
#' \code{\link[lumpR]{prof_class}}. If omitted, the existing raster layer \code{lu} is used.
#' @param lu Input or Output: Name of Landscape Units (LU) raster map already existing in
#' GRASS location or to be generated using \code{recl_lu}.
#' @param eha Name of Elementary Hillslope Areas (EHA) raster map in GRASS
#' location. Output of \code{\link[lumpR]{lump_grass_prep}}.
#' @param fill_holes TRUE: fill any holes in map \code{eha} (e.g. as result of skipped/non-classified
#' EHAs) by growing
#' @param flowacc Name of flow accumulation raster map in GRASS location. Can
#' be created with \code{\link[lumpR]{lump_grass_prep}}.
#' @param flowdir Name of flow direction raster map in GRASS location. Can
#' be created with \code{\link[lumpR]{lump_grass_prep}}.
#' @param stream_horton Name of Horton stream order raster map in GRASS location. Can
#' be created with \code{\link[lumpR]{lump_grass_prep}}. If left empty, the channel length,
#' slope and retention times are set to \code{NA}.
#' @param soil_depth Name of soil depth [mm] raster map in GRASS location. If \code{NULL}
#' (default), \code{na_val} is used.
#' @param sdr Name of sediment delivery ratio [-] raster map in GRASS location. If empty,
#' this optional column is omitted.
#' @param dir_out Character string specifying output directory (will be created;
#' any overwriting will be prompted).
#' @param sub_ofile Output: Name of subbasin statistics file containing subbasin
#' parameters. See \code{Details} below. If \code{NULL} (default) this file
#' will not be generated.
#' @param lu_ofile Output: Name of file containing subbasins and the corresponding
#' LUs with their fraction of area in the subbasin. If \code{NULL} (default)
#' this file will not be generated.
#' @param lupar_ofile Output: Name of file containing LUs and related parameters.
#' See \code{Details} below. If \code{NULL} (default) this file will not be generated.
#' @param groundwater Flag: 1: respect groundwater and infer parameters. 0 (default):
#' Ignore groundwater and associated parameters.
#' @param na_val Value used to indicate \code{NA} values in output files. Default: \code{NA}.
#' @param keep_temp \code{logical}. Set to \code{TRUE} if temporary files shall be kept
#' in the GRASS location, e.g. for debugging or further analyses. Default: \code{FALSE}.
#' @param overwrite \code{logical}. Shall output of previous calls of this function be
#' deleted? If \code{FALSE} the function returns an error if output already exists.
#' Default: \code{FALSE}.
#' @param silent \code{logical}. Shall the function be silent (also suppressing warnings
#' of internally used GRASS functions)? Default: \code{FALSE}.
#'
#' @return Function returns nothing. Output files (\code{sub_ofile, lupar_ofile, lu_ofile})
#' are written into output directory and raster map (\code{lu}) exported into GRASS
#' location as specified in arguments.
#'
#' @note Prepare GRASS location and necessary raster files in advance (e.g. using
#' \code{\link[lumpR]{lump_grass_prep}}) and start GRASS session in R using
#' \code{\link[rgrass7]{initGRASS}}.
#'
#' \bold{IMPORTANT:} Herein, when specifying the GRASS input maps, please do
#' explicitly refer to the mapset if it is different from the mapset given in
#' initGRASS() (even PERMANENT!), as otherwise internally used readRAST() command
#' resulted in errors under Windows.
#'
#' TODO:\cr
#' - check empirical formulas for channel width and channel depth\cr
#' - LU parameter estimation\cr
#' - include options to add parameters manually in case data are available\cr
#' - include option to make function more efficient regarding RAM usage (e.g. by
#' wrting/reading temporary raster data to/from disk) at the cost of higher computational burden
#'
#'
#' @details
#' \bold{Subbasin parameters}\cr
#' Subbasin parameter estimation given in \code{sub_ofile} contains:
#'
#' \emph{pid}\cr
#' Subbasin identifier.
#'
#' \emph{description}\cr
#' Subbasin description you can fill in manually if needed.
#'
#' \emph{lat}\cr
#' Latitude of subbasin centroid in \emph{decimal degrees} (negative values for southern hemisphere).
#'
#' \emph{lon}\cr
#' Longitude of subbasin centroid in \emph{decimal degrees west of Greenwhich}, e.g.
#' Greenwich: 0°, New York: 75°, Berlin: 345°.
#'
#' \emph{elev}\cr
#' Average elevation above sea level of subbasin \emph{m}.
#'
#' \emph{area}\cr
#' Subbasin area in \emph{km^2}.
#'
#' \emph{x}\cr
#' X-coordinate of subbasin centroid in units of GRASS location.
#'
#' \emph{y}\cr
#' Y-coordinate of subbasin centroid in units of GRASS location.
#'
#' \emph{drains_to}\cr
#' Pid of subbasin the current subbasin drains to. It is determined by identifying
#' the cell with the highest flow accumulation value (= subbasin outlet). By
#' deriving the drainage direction value of that raster cell the next downstream
#' cell is identified along with corresponding subbasin number.
#'
#' \emph{lag_time}\cr
#' Time in \emph{days} a runoff signal in the current subbasin needs to be directed
#' from the subbasin input to the outlet. Estimated from channel geometry (see below).
#'
#' \emph{retention}\cr
#' Maximum time period in \emph{days} over which a runoff signal is distributed by
#' the routing process. Estimated from channel geometry (see below).
#'
#' \emph{chan_len}\cr
#' Length of the main channel of the respective subbasin in \emph{m}. Estimated
#' from channel geometry (see below).
#'
#' \emph{channel geometry}\cr
#' Main channel length: For each subbasin the main channel is determined from Horton
#' stream raster map. Its length is then calculated depending on raster resolution
#' and flow direction.\cr
#' Channel slope: Minimum (= inflow) and maximum (= outflow) flow accumulation and
#' elevation of the corresponding raster cells are determined. The difference in
#' elevation is divided by main channel length.\cr
#' Channel width: Maximum flow accumulation is determined and the corresponding
#' drainage area calculated by resolution of raster cells. Channel width is then
#' calculated from the empirical formula: width[m] = 1.29 * darea[km2] ^ 0.6.\cr
#' Channel depth: Empirical formula based on drainage area: depth[m] = 0.13 * darea[km2] ^ 0.4.\cr
#' Flow velocity: Is calculated using Mannings equation with an n-value of 0.075.
#' Flow velocities are calculated to derive flow travel times (velocity / channel
#' length) for bankfull (= high flow condition), 2/3 (= average conditions) and
#' 1/10 (= low flow conditions) water levels to derive lag time (travel time for
#' average conditions) and retention time (max - min travel time).
#'
#'
#' \bold{Landscape Units in Subbasins}\cr
#' \code{lu_ofile} contains:
#'
#' \emph{subbas_id}\cr
#' Subbasin identifier.
#'
#' \emph{lu_id}\cr
#' Landscape Unit identifier.
#'
#' \emph{fraction}\cr
#' Areal fraction of Landscape Unit within corresponding Subbasin.
#'
#'
#' \bold{Landscape Unit parameters}\cr
#' Landscape Unit parameter estimation given in \code{lupar_ofile} contains:
#'
#' \emph{pid}\cr
#' Landscape Unit identifier.
#'
#' \emph{description}\cr
#' Description for this Landscape Unit. Can be adjusted manually if you want.
#' Generally set to \code{na_val}.
#'
#' \emph{kf_bedrock}\cr
#' Hydraulic conductivity of bedrock. Fill in values manually. Generally set
#' to \code{na_val}. Use in WASA model is optional (see WASA documentation
#' -> 'kfsu' and notes on 'bedrock').
#'
#' \emph{slopelength}\cr
#' Slope length of Landscape Unit [m]. Value can be obtained from output of
#' \code{\link[lumpR]{prof_class}}. Herein set to \code{na_val}.
#'
#' \emph{soil_depth}\cr
#' Soil depth in \emph{mm} averaged over respective landscape unit. Use in
#' WASA model is optional (see WASA documentation -> 'meandep' and notes on
#' 'bedrock'). Herein set to \code{na_val}.
#'
#' \emph{allu_depth}\cr
#' Depth of alluvial soils in \emph{mm}. Use in WASA model is optional (see
#' WASA documentation -> 'maxdep' and notes on 'bedrock'). Herein set to \code{na_val}.
#'
#' \emph{riverbed_depth}\cr
#' Depth of river bed below Terrain Component in \emph{mm}. Use in WASA model
#' is optional (see WASA documentation -> 'riverbed' and notes on 'bedrock').
#' Herein set to \code{na_val}.
#'
#' \emph{gw_flag}\cr
#' Groundwater flag: 0: no groundwater in this LU. 1: LU contains groundwater.
#' At the moment set to 0 by default for every LU. Use in WASA model is optional
#' (see WASA documentation -> 'gw_flag' and notes on groundwater).
#'
#' \emph{gw_dist}\cr
#' Initial depth of groundwater below surface in \emph{mm}. At the moment set to
#' 1000 by default for every LU if \code{groundwater = 1}. Use in WASA model is optional
#' (see WASA documentation -> 'gw_dist' and notes on groundwater).
#'
#' \emph{frgw_delay}\cr
#' Storage coefficient for groundwater outflow in \emph{days}. At the moment set to
#' 200 by default for every LU if \code{groundwater = 1}. Use in WASA model is optional
#' (see WASA documentation -> 'frgw_delay' and notes on groundwater).
#'
#' \emph{sdr_lu} (optional)\cr
#' sediment delivery ratio from raster map \code{sdr}, if specified
#'
#'
#' @references Source code based on \code{SHELL} and \code{MATLAB} scripts of Till Francke.
#'
#' lumpR package introduction with literature study and sensitivity analysis:\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
#'
#' Theory of LUMP:\cr
#' Francke, T.; Guentner, A.; Mamede, G.; Mueller, E. N. and Bronstert, A (2008):
#' Automated catena-based discretization of landscapes for the derivation of
#' hydrological modelling units. \emph{International Journal of Geographical
#' Information Science, Informa UK Limited}, 22(2), 111-132, DOI: 10.1080/13658810701300873
#'
#' Subbasin Parameters:\cr
#' Bronstert, A., Guentner, A., Jaeger, A., Krol, M. & Krywkow, J. (1999): Grossraeumige
#' hydrologische Parametrisierung und Modellierung als Teil der integrierten
#' Modellierung. In: Fohrer, N. & Doell, P. (Eds.): Modellierung des Wasser- und
#' Stofftransports in grossen Einzugsgebieten. \emph{Kassel University Press}, Kassel,
#' Germany, 31-40.
#'
#' Guentner, A. (2002): Large-scale hydrological modelling in the semi-arid
#' North-East of Brazil. \emph{PIK Report 77}, Potsdam Institute for Climate
#' Impact Research, Potsdam, Germany.
#'
#'
#' @author Tobias Pilz \email{tpilz@@uni-potsdam.de}
lump_grass_post <- function(
### INPUT ###
mask=NULL,
dem=NULL,
recl_lu=NULL,
lu=NULL,
subbasin=NULL,
eha=NULL,
flowacc=NULL,
flowdir=NULL,
stream_horton=NULL,
soil_depth=NULL,
sdr=NULL,
### OUTPUT ###
dir_out="./",
sub_ofile=NULL,
lu_ofile=NULL,
lupar_ofile=NULL,
### OPTIONS ###
fill_holes=T,
groundwater=0,
na_val=NA,
keep_temp=F,
overwrite=F,
silent=F
) {
### PREPROCESSING ###----------------------------------------------------------
if(!silent) message("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
if(!silent) message("% START lump_grass_post()")
if(!silent) message("%")
if(!silent) message("% Initialise function...")
# checks #---------------------------------------------------------------------
tryCatch(gmeta(), error = function(e) stop("Cannot execute GRASS commands. Maybe you forgot to run initGRASS()?"))
if(is.null(mask))
stop("The name of a raster within the mapset of your initialised GRASS session to be used as catchment MASK in GRASS has to be given!")
if(!is.character(mask))
stop("The name of a raster to be used as catchment MASK in GRASS must be a string (argument 'mask')!")
if(is.null(dem))
stop("The name of a DEM within the mapset of your initialised GRASS session has to be given!")
if(is.null(subbasin))
stop("A name for needed subbasin raster map within the mapset of your initialised GRASS session has to be given!")
if(is.null(eha))
stop("A name for needed EHA raster map within the mapset of your initialised GRASS session has to be given!")
if(is.null(flowdir))
stop("A name for needed flow direction raster map within the mapset of your initialised GRASS session has to be given!")
if(is.null(flowacc))
stop("A name for needed flow accumulation raster map within the mapset of your initialised GRASS session has to be given!")
if(is.null(lu))
stop("A name for the landscape units raster map to be generated (recl_lu = NULL) or already within the mapset of your initialised GRASS session has to be given")
if(!is.null(recl_lu) && !file.exists(recl_lu))
stop(paste0("Could not find file '",recl_lu,"' (argument 'recl_lu')"))
#check existence of raster maps
check_raster(mask,"mask")
check_raster(dem,"dem")
if(is.null(recl_lu)) check_raster(lu,"lu")
check_raster(subbasin,"subbasin")
check_raster(eha,"eha")
check_raster(flowacc,"flowacc")
check_raster(flowdir,"flowdir")
check_raster(stream_horton,"stream_horton")
if(!is.null(soil_depth)) check_raster(soil_depth,"soil_depth")
if(!is.null(sdr)) check_raster(sdr,"sdr")
# suppress annoying GRASS outputs
tmp_file <- file(tempfile(), open="wt")
sink(tmp_file, type="output")
# also supress warnings in silent mode
if(silent){
tmp_file2 <- file(tempfile(), open="wt")
sink(tmp_file2, type="message")
oldw <- getOption("warn")
options(warn = -1)
}
# create output dir
dir.create(dir_out, recursive=T, showWarnings=F)
# check output directory
if (!overwrite & (file.exists(paste(dir_out,sub_ofile,sep="/")) |
file.exists(paste(dir_out,lu_ofile,sep="/")) |
file.exists(paste(dir_out,lupar_ofile,sep="/"))) )
stop(paste0("Output file(s) ", sub_ofile, ", ",lu_ofile, ", and/or ", lupar_ofile, " already exist(s) in ", dir_out, "!"))
# remove mask if there is any (and ignore error in case there is no mask)
# set mask to make sure calculations are done exactly within expected area
cmd_out <- execGRASS("g.copy", raster=paste0(mask,",MASK"), flags = "overwrite", intern = T)
# remove output of previous function calls if overwrite=T
if (overwrite) {
cmd_out <- execGRASS("g.remove", type="raster", pattern=paste("*_t,", lu, sep=","), flags=c("f", "b"), intern=T)
} else {
# remove temporary maps in any case
cmd_out <- execGRASS("g.remove", type="raster", pattern="*_t", flags=c("f", "b"), intern=T)
}
if(!silent) message("% OK")
### CALCULATIONS ###-----------------------------------------------------------
# prepare output directory, read GRASS data #-----------------------------
tryCatch({
if(!silent) message("%")
if(!silent) message("% Prepare output and import GRASS data...")
# load rasters into R
dem_rast <- read_raster(dem)
accum_rast <- read_raster(flowacc)
dir_rast <- read_raster(flowdir)
cur_mapset = execGRASS("g.mapset", flags="p", intern=TRUE) #determine name of current mapset
sub_rast <- read_raster(subbasin)
horton_rast <- read_raster(stream_horton)
# ... raster values as matrix
dem_mat <- getValues(dem_rast, format="matrix")
rm(dem_rast)
sub_mat <- getValues(sub_rast, format="matrix")
coords <- xyFromCell(sub_rast, 1:length(sub_rast))
rm(sub_rast)
accum_mat <- getValues(accum_rast, format="matrix")
rm(accum_rast)
dir_mat <- getValues(dir_rast, format="matrix")
rm(dir_rast)
horton_mat <- getValues(horton_rast, format="matrix")
rm(horton_rast)
gc(); gc()
# resolution
g_meta <- gmeta()
resol <- mean(g_meta$nsres, g_meta$ewres)
if(!silent) message("% OK")
}, error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
stop(paste(e))
})
# generate lu raster #---------------------------------------------------------
tryCatch({
if(!silent) message("%")
if (is.null(recl_lu)) #use existing LU-map
{
if(!silent) message(paste0("% NOTE: Use existing raster map ", lu, " as specified."))
} else {
if(!silent) message(paste0("% Generate landscape units raster map '", lu, "'..."))
#reclass EHA according to reclass file generated by prof_class to get LU
cmd_out=execGRASS("r.reclass", input=eha, output=lu, rules=recl_lu, intern = TRUE)
if(!is.null(attr(cmd_out, "status")))
stop(paste0("Could not reclass '",eha,"' into '",lu,"'"))
# growing radius (parameter for r.grow)
GROWRAD <- 20
if (fill_holes)
while (TRUE)
{
# look for empty patches
cmd_out=execGRASS("r.mapcalc", expression=paste0("grow_eval_t = if(isnull(", lu, "), 1, 0)"), flags = "overwrite", intern=TRUE)
na_eval <- execGRASS("r.stats", input=paste0("grow_eval_t"), flags=c("n", "p"), intern=TRUE, ignore.stderr = T)
if (length(na_eval)==0 ) break #no empty patches found
# if (grepl(pattern="[0-9]+.*[\b]+",x=tail(na_eval, n=1)))
# na_eval = na_eval[-length(na_eval)] #last line contains progress indicator, remove
na_line=grepl(na_eval, pattern = "^1 .*\\%$")
if (!any(na_line)) break #no empty patches found
na_res <- strsplit(x=na_eval[na_line], split=" +")
if(!silent) message(paste0("% -> ", tail(unlist(na_res),1), " NAs in LU-map '", lu,"'"))
# otherwise grow LU map to fill gaps
if(!silent) message(paste0("% -> Filling empty patches in ",lu,", iterating..."))
cmd_out=execGRASS("r.grow", input=lu, output="lu_t", radius=GROWRAD, flags="overwrite", intern=TRUE)
if(!is.null(attr(cmd_out, "status")))
stop(paste0("Could not grow '",lu,"'"))
# r.grow converts type CELL to type DCELL; convert back to CELL
cmd_out <- execGRASS("r.mapcalc", expression=paste0(lu," = int(lu_t)"), flags = "overwrite", intern = T)
cmd_out <- execGRASS("g.remove", type="raster", pattern="lu_t", flags=c("f", "b"), intern=T)
# x=execGRASS("r.mapcalc", expression = paste0(lu, " = ", mask, " * tt_eha_grown"), flags="overwrite", intern=TRUE)
#
# x=execGRASS("g.remove", type="raster", name="tt_eha_grown", intern=TRUE)
} # while true (grow-loop)
cmd_out=execGRASS("g.remove", type="raster", name="grow_eval_t", intern=TRUE)
if(!silent) message("% OK")
} # recl_lu was given
},
error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
# remove output so far
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern=paste0("grow_eval_t,", lu), flags=c("f", "b"), intern=TRUE)
stop(paste(e))
}
)
# subbasin statistics #--------------------------------------------------------
tryCatch({
if(!is.null(sub_ofile)) {
if(!silent) message("%")
if(!silent) message("% Calculate subbasin statistics...")
sub_stats <- execGRASS("r.stats", input=subbasin, flags=c("a", "n"), intern=TRUE, ignore.stderr = T)
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(sub_stats, n=1)))
sub_stats <- sub_stats[-length(sub_stats)] #last line contains progress indicator, remove
sub_stats <- matrix(as.numeric(gsub("%", "", unlist(strsplit(sub_stats, split=" +")))), ncol=2, byrow=T,
dimnames=list(NULL,c("pid", "area")))
# convert m2 to km2
sub_stats[,"area"] <- round(sub_stats[,"area"]/1e6, 2)
# calculate stats of LUs in each subbasin and subbasin drainage ("drains_to")
sub_stats <- cbind(sub_stats, na_val, na_val, na_val, na_val, na_val, na_val, na_val, na_val, na_val, na_val, na_val)
colnames(sub_stats)[c(3:13)] <- c("x", "y", "lon", "lat", "elev", "drains_to", "lag_time", "retention", "description", "a_stream_order", "chan_len")
sub_lu_stats <- NULL
s_row <- 0
for (SUB in sub_stats[,1]) {
s_row <- s_row +1
# bounding box of current subbasin
sub_arrind <- which(sub_mat == SUB, arr.ind = T)
rmin <- min(sub_arrind[,"row"])
rmax <- max(sub_arrind[,"row"])
cmin <- min(sub_arrind[,"col"])
cmax <- max(sub_arrind[,"col"])
# crop raster matrices
dem_crop <- dem_mat[rmin:rmax, cmin:cmax]
horton_crop <- horton_mat[rmin:rmax, cmin:cmax]
sub_crop <- sub_mat[rmin:rmax, cmin:cmax]
dir_crop <- dir_mat[rmin:rmax, cmin:cmax]
accum_crop <- accum_mat[rmin:rmax, cmin:cmax]
# get coordinates of cells
coords_crop <- coords[which(t(sub_mat) == SUB),]
# set all values not overlapping with the current subbasin to NA
na_cells <- which(sub_crop != SUB)
sub_crop[na_cells] <- NA
dem_crop[na_cells] <- NA
horton_crop[na_cells] <- NA
dir_crop[na_cells] <- NA
accum_crop[na_cells] <- NA
# COORDINATES OF SUBBASIN centroids in GRASS units #
sub_centr <- centroid(coords_crop)
sub_stats[s_row, "x"] <- round(sub_centr[,"x"])
sub_stats[s_row, "y"] <- round(sub_centr[,"y"])
# convert to decimal degree
sub_centr <- data.frame(x = sub_centr[,"x"], y = sub_centr[,"y"])
coordinates(sub_centr) <- c("x", "y")
proj4string(sub_centr) <- getLocationProj()
longlat <- spTransform(sub_centr, "+proj=longlat")
sub_stats[s_row, "lon"] <- round(coordinates(longlat)[,"x"], 4)
sub_stats[s_row, "lat"] <- round(coordinates(longlat)[,"y"], 4)
# longitude as decimal degrees west of Greenwich
if(sub_stats[s_row, "lon"] < 0)
sub_stats[s_row, "lon"] <- -1*sub_stats[s_row, "lon"]
else
sub_stats[s_row, "lon"] <- 360 - sub_stats[s_row, "lon"]
# AVERAGE SUBBASIN ELEVATION #
sub_stats[s_row, "elev"] <- round(mean(dem_crop, na.rm=T))
# SUBBASIN drainage #
sub_stats[s_row,"drains_to"] <- sub_route(SUB,sub_mat,accum_mat,dir_mat) # internal function, see below
# SUBBASIN PARAMETERS #
# calc main channel length
chan_len <- channel_length(horton_crop, dir_crop, resol)
# calc main channel average slope
chan_slope <- channel_slope(horton_crop, accum_crop, dem_crop, chan_len)
# calc main channel width
chan_width <- channel_width(accum_crop, resol)
# calc main channel depth
chan_depth <- channel_depth(accum_crop, resol)
# calc flow velocites for bankful, 2/3 and 1/10 filling
# Manning's n = 0.075 (very weedy reaches, deep pools, or floodways with heavy stand of timber and underbrush)
# http://www.fsl.orst.edu/geowater/FX3/help/8_Hydraulic_Reference/Mannings_n_Tables.htm
chan_velo_full <- flow_velocity(chan_width, chan_depth, chan_slope, n=0.075)
chan_velo_med <- flow_velocity(chan_width, 2/3 * chan_depth, chan_slope, n=0.075)
chan_velo_low <- flow_velocity(chan_width, 1/10 * chan_depth, chan_slope, n=0.075)
# calc travel times (translation and retention) [d]
# approach according to Bronstert et al. (1999), Guentner (2002)
flowtime_min <- chan_len / chan_velo_full / 86400
flowtime_med <- chan_len / chan_velo_med / 86400
flowtime_max <- chan_len / chan_velo_low / 86400
# check results
retention <- flowtime_max - flowtime_min
if(!is.finite(flowtime_med) || !is.finite(retention))
warning(paste0("Could not calculate finite subbasin parameters for subbasin ", SUB))
# save
sub_stats[s_row, "lag_time"] <- round(flowtime_med, 3)
sub_stats[s_row, "retention"] <- round(retention, 3)
sub_stats[s_row, "chan_len"] <- round(chan_len)
} # subbasin loop
# write output file
write.table(sub_stats, paste(dir_out, sub_ofile, sep="/"), quote=F, row.names=F, sep="\t")
if(!silent) message("% OK")
} # sub_ofile given?
}, error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
# remove output so far
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern=paste0("*_t,", lu), flags=c("f", "b"), intern=TRUE)
stop(paste(e))
})
# lu statistics #--------------------------------------------------------------
tryCatch({
if(!is.null(lu_ofile)) {
if(!silent) message("%")
if(!silent) message("% Calculate LU statistics...")
sub_lu_stats_t <- execGRASS("r.stats", input=paste0(subbasin,",",lu), flags=c("n", "c"), intern=TRUE, ignore.stderr = T)
if(!is.null(attr(sub_lu_stats_t, "status")))
stop(paste0("Could not access lu-map '",lu,"'."))
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(sub_lu_stats_t, n=1))) #check if last line contains progress indicator, remove
sub_lu_stats_t = sub_lu_stats_t[-length(sub_lu_stats_t)]
sub_lu_stats_t = sub_lu_stats_t[grepl(x=sub_lu_stats_t, pattern = "^[0-9]")] #remove non-number lines
sub_lu_stats_t = as.numeric(gsub("%", "", unlist(strsplit(sub_lu_stats_t, split=" +")))) #convert to numbers and remove %
sub_lu_stats_t2 =data.frame(matrix(sub_lu_stats_t, ncol=3, byrow=T))
names(sub_lu_stats_t2)=c("subbas_id","lu_id","cells")
subbasin_cells=aggregate(x = sub_lu_stats_t2$cells, by = list(subbas_id=sub_lu_stats_t2$subbas_id), FUN=sum)
sub_lu_stats_t2=merge(sub_lu_stats_t2,subbasin_cells)
sub_lu_stats_t2$fraction=sub_lu_stats_t2$cells / sub_lu_stats_t2$x #compute fraction
write.table(sub_lu_stats_t2[,c("subbas_id", "lu_id", "fraction")], paste(dir_out, lu_ofile, sep="/"), quote=F, row.names=F, sep="\t")
if(!silent) message("% OK")
} # lu_ofile given?
}, error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
# remove output so far
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern=paste0("*_t,", lu), flags=c("f", "b"), intern=TRUE)
stop(paste(e))
})
# lu parameters #--------------------------------------------------------------
tryCatch({
if(!is.null(lupar_ofile)) {
if(!silent) message("%")
if(!silent) message("% Calculate LU parameters")
# identify LUs
lu_ids <- execGRASS("r.stats", input=lu, flags=c("n"), intern=TRUE, ignore.stderr = T)
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(lu_ids, n=1)))
lu_ids <- as.numeric(lu_ids[-length(lu_ids)]) #last line contains progress indicator, remove
# initialise LU parameter output object
lu_par <- matrix(na_val, ncol=11, nrow=length(lu_ids),
dimnames=list(NULL, c("pid", "description", "kf_bedrock", "slopelength",
"soil_depth", "allu_depth", "riverbed_depth",
"gw_flag", "gw_dist", "frgw_delay", "sdr_lu")))
# calculate mean soil depth for every LU
if (!is.null(soil_depth) && soil_depth!="")
{
cmd_out <- execGRASS("r.univar", zones=lu, map=soil_depth, fs=",", flags=c("t"),intern=T, ignore.stderr = T)
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(cmd_out, n=1)))
cmd_out <- cmd_out[-length(cmd_out)] #last line contains progress indicator, remove
cmd_out <- strsplit(cmd_out, ",")
cmd_out <- matrix(unlist(cmd_out[-1]), ncol=length(cmd_out[[1]]), byrow=T,
dimnames=list(NULL, cmd_out[[1]]))
lu_depth <- as.numeric(cmd_out[,"mean"])
lu_ids = as.numeric(cmd_out[,"zone"])
# quick check of soil depths
if(any(lu_depth > 10000))
warning("There are average LU soil depths of more than 10 m which is a large (but not impossible) value. Check your input data and units ('soil_depth' in [mm])!")
if(any(lu_depth < 100))
warning("There are average LU soil depths of less than 10 cm which is a small (but not impossible) value. Check your input data and units ('soil_depth' in [mm])!")
} else lu_depth=na_val
lu_par[,"pid"] <- lu_ids
lu_par[,"soil_depth"] <- lu_depth
if (!is.null(sdr) && sdr!="")
{
cmd_out <- execGRASS("r.univar", zones=lu, map=sdr, fs=",", flags=c("t"),intern=T, ignore.stderr = T)
if (grepl(pattern="[0-9]+.*[\b]+",x=tail(cmd_out, n=1)))
cmd_out <- cmd_out[-length(cmd_out)] #last line contains progress indicator, remove
cmd_out <- strsplit(cmd_out, ",")
cmd_out <- matrix(unlist(cmd_out[-1]), ncol=length(cmd_out[[1]]), byrow=T,
dimnames=list(NULL, cmd_out[[1]]))
sdr_vals <- as.numeric(cmd_out[,"mean"])
lu_ids = as.numeric(cmd_out[,"zone"])
sdr_lu <- sdr_vals[match(lu_par[,"pid"], lu_ids)]
lu_par[,"sdr_lu"] <- sdr_lu
}
# groundwater parameters (so far only default values)
# TODO: alternative approaches?
# groundwater for every LU
if(groundwater) {
lu_par[,"gw_flag"] <- rep(1, length(lu_depth))
# initial depth of groundwater below surface: 1000 mm
lu_par[,"gw_dist"] <- rep(1000, length(lu_depth))
# storage coefficient for groundwater outflow [days]; TODO: estimate from baseflow analysis?!
lu_par[,"frgw_delay"] <- rep(200, length(lu_depth))
} else {
lu_par[,"gw_flag"] <- rep(0, length(lu_depth))
# initial depth of groundwater below surface: 1000 mm
lu_par[,"gw_dist"] <- rep(na_val, length(lu_depth))
# storage coefficient for groundwater outflow [days]; TODO: estimate from baseflow analysis?!
lu_par[,"frgw_delay"] <- rep(na_val, length(lu_depth))
}
# write output
write.table(lu_par, paste(dir_out, lupar_ofile, sep="/"), quote=F, row.names=F, sep="\t")
if(!silent) message("% OK")
} # lupar_ofile fiven?
}, error = function(e) {
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
# remove mask
tryCatch(suppressWarnings(execGRASS("r.mask", flags=c("r"))), error=function(e){})
# remove output so far
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern=paste0("*_t,", lu), flags=c("f", "b"), intern=TRUE)
stop(paste(e))
})
# remove temp files
if(!keep_temp)
x=execGRASS("g.remove", type="raster", pattern="*_t", flags=c("f", "b"), intern=TRUE)
if(!silent) message("%")
if(!silent) message("% DONE!")
if(!silent) message("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
# stop sinking
closeAllConnections()
# restore original warning mode
if(silent)
options(warn = oldw)
} # EOF
### INTERNAL FUNCTIONS ###-----------------------------------------------------
# subbasin centroids #---------------------------------------------------------
# returns a matrix with x and y value of the centroid of the object given in coords_sub
centroid <- function(coords_sub) {
# convert so sp class
coords_sub <- as.data.frame(coords_sub)
coordinates(coords_sub) <- c("x", "y")
# calculate centroid
centr <- gCentroid(coords_sub)
centr <- centr@coords
return(centr)
}
# subbasin routing #-----------------------------------------------------------
# returns ID of downstream subbasin for the current subbasin 'sub_no'
# determined from flow accumulation and flow direction map
sub_route <- function(sub_no,sub_mat,accum_mat,dir_mat) {
# extract highest flowacc in subbasin sub_no
sub_ids <- which(sub_mat == sub_no)
accum_sub <- accum_mat[sub_ids]
accum_sub_max <- which.max(accum_sub)
# extract corresp. flowdir
dir_sub_out <- dir_mat[sub_ids[accum_sub_max]]
# determine row and col no in sub_mat
sub_rowcol <- arrayInd(sub_ids[accum_sub_max], dim(sub_mat))
# determine value of sub_mat (subbasin no.) of neighbour cell according to dir_sub_out
if (dir_sub_out < 0) { # outlet of catchment
return(9999)
} else {
if (dir_sub_out == 1) { # NE
sub_rowcol_out <- sub_rowcol + c(-1,1)
} else if (dir_sub_out == 2) { # N
sub_rowcol_out <- sub_rowcol + c(-1,0)
} else if (dir_sub_out == 3) { # NW
sub_rowcol_out <- sub_rowcol + c(-1,-1)
} else if (dir_sub_out == 4) { # W
sub_rowcol_out <- sub_rowcol + c(0,-1)
} else if (dir_sub_out == 5) { # SW
sub_rowcol_out <- sub_rowcol + c(+1,-1)
} else if (dir_sub_out == 6) { # S
sub_rowcol_out <- sub_rowcol + c(+1,0)
} else if (dir_sub_out == 7) { # SE
sub_rowcol_out <- sub_rowcol + c(+1,+1)
} else if (dir_sub_out == 8) { # E
sub_rowcol_out <- sub_rowcol + c(0,+1)
} else {
stop(paste("While determining subbasin drainage: Determined flow direction at outlet
of subbasin ", sub_no, " has value ", dir_sub_out, " but should be one of {1,2,3,4,5,6,7,8} or a negative number.", sep=""))
}
# determine subbasin the current subbasin drains to
sub_out <- sub_mat[sub_rowcol_out]
# update stats table
if (is.na(sub_out)) { # catchment outlet
return(9999)
} else {
if (sub_out == sub_no)
stop(paste0("In subbasin no. ", sub_no, " the determined downstream subbasin is this subbasin!"))
return(sub_out)
}
}
} # EOF
# main channel length #--------------------------------------------------------
# returns length of the main channel (largest value in Horton order) in [m]
# computes main stream temporary raster used in further calculations
channel_length <- function(horton, flowdir, resol) {
# determine main channel (largest value in Horton stream order)
main_chan <- na.exclude(unique(c(horton)))
# if there is no main channel (e.g. in very small reservoir subbasins) assume one cell of main stream
if (length(main_chan) == 0) {
chan_len <- resol
# accmax <- max(as.numeric(execGRASS("r.stats", input=flowacc, flags=c("n", "1"), intern=T, ignore.stderr = T))) # max flowacc
# execGRASS("r.mapcalculator", amap=flowacc, outfile="stream_main_t",
# formula=paste0("if(A == ", sprintf(accmax, fmt="%d"), ", 1, null())"))
warning(paste("Subbasin ", sub_no, " has no main channel. Assume at least one raster cell.", sep=""))
} else {
# diagonal raster cell length
dia <- sqrt(2*resol^2)
# determine raster cells of main channel
max_val <- max(main_chan)
r_main <- which(horton==max_val)
# calculate length of the main channel
chanlen_t <- flowdir[r_main] %% 2 == 0 # horizontal/vertical stream
chan_len <- length(which(chanlen_t))*resol + length(which(!chanlen_t))*dia
}
return(chan_len)
} # EOF
# main channel slope #---------------------------------------------------------
# returns average slope of main channel [m/m]
channel_slope <- function(horton, accum, dem, chan_len) {
if (is.na(chan_len) ) return(NA)
# determine main channel (largest value in Horton stream order)
main_chan <- na.exclude(unique(c(horton)))
# determine raster cells of main channel
max_val <- max(main_chan)
r_main <- which(horton==max_val)
# min and max values of flow accumulation within main channel to determin flow direction
acc_min <- which.min(accum[r_main])
acc_max <- which.max(accum[r_main])
# get dem values for respective raster cells of min and max flowacc
dem_vals <- dem[r_main][c(acc_max, acc_min)]
# compute slope [m/m]
chan_slope <- diff(dem_vals) / chan_len
# slope zero not allowed -> at least 0.00001
chan_slope <- max(chan_slope, 0.00001)
return(chan_slope)
}
# main channel width #---------------------------------------------------------
# calculate main channel width based on empirical formula: width[m] = 1.29 * darea[km2] ^ 0.6
# darea = drainage area determined from maximum flow accumulation and resolution
channel_width <- function(accum, resol) {
# determine maximum flow accumulation
max_acc <- max(accum, na.rm = T)
# calculate drainage area [km^2]
drain_area <- max_acc * resol^2 / 1e6
# calculate width according to empirical function
chan_width <- 1.29 * drain_area^(0.6)
if(chan_width > 3000)
stop(paste("Calculated channel width is ", chan_width, "m which seems unrealistic!", sep=""))
return(chan_width)
}
# main channel depth #---------------------------------------------------------
# calculate main channel width based on empirical formula: width[m] = 0.13 * darea[km2] ^ 0.4
# darea = drainage area determined from maximum flow accumulation and resolution
channel_depth <- function(accum, resol) {
# determine maximum flow accumulation
max_acc <- max(accum, na.rm = T)
# calculate drainage area [km^2]
drain_area <- max_acc * resol^2 / 1e6
# calculate depth according to empirical function
chan_depth <- 0.13 * drain_area^(0.4)
if(chan_depth > 50)
stop(paste("Calculated channel depth is ", chan_depth, "m which seems unrealistic!", sep=""))
return(chan_depth)
}
# flow velocity #--------------------------------------------------------------
# function returns flow velocity of main channel [m/s]
# calc by Gauckler-Manning-Strickler
# assume simple rectangle cross section as width >> depth for main channels
flow_velocity <- function(chan_width, chan_depth, chan_slope, n) {
# calculate cross section area (assume simple rectangle as width >> depth for main channels)
cross_area <- chan_width * chan_depth
# calculate hydraulic radius (wetted perimeter = width + 2 * depth)
hyd_radius <- cross_area / (chan_width + 2 * chan_depth)
# calculate flow velocity
chan_velo <- 1/n * hyd_radius^(2/3) * chan_slope^(1/2)
return(chan_velo)
}
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.