R/calc_attributes_sites_approx.R

In openSTARS: An Open Source Implementation of the 'ArcGIS' Toolbox 'STARS'

#' Calculate attributes of the sites.
#'
#' For each site (observations or predictions) attributes (potential predictor variables)
#' are derived based on the values calculated for the edge the site lies on.
#' This function calculates approximate values for site catchments as described
#' in Peterson & Ver Hoef, 2014: STARS: An ArcGIS Toolset Used to Calculate the
#' Spatial Information Needed to Fit Spatial Statistical Models to Stream
#' Network Data. J. Stat. Softw., 56 (2).
#'
#' @param sites_map character; name of the sites (observation or prediction) 
#'  attributes shall be calculated for. "sites" (default) refers to the observation sites.
#' @param input_attr_name character vector; input column name(s) in the
#'   attribute table of the vector map "edges".
#' @param output_attr_name character vector (optional); output column name(s)
#'   appended to the site attribute data table. If not provided it is set to
#'   \code{input_attr_name}. Attribute names must not be longer than 10
#'   characters.
#' @param stat name or character vector giving the statistics to be calculated.
#'   See details below.
#' @param round_dig integer; number of digits to round results to.
#' @param calc_basin_area boolean; shall the catchment area be calculated? (Useful
#'  to set to FALSE if the function has been called before.)
#' @param overwrite boolean; shall existing columns be overwritten; defaults to FALSE
#'
#' @return Nothing. The function appends new columns to the \code{sites_map}
#'   attribute table
#' \itemize{
#'  \item{'H2OAreaA':} {Total watershed area of the watershed upstream of each site.}
#'  \item{attr_name:} {Additional optional attributes calculated based on \code{input_attr_name}.}
#' }
#'
#' @details The approximate total catchment area (H2OAreaA) is always calculated
#'   if \code{calc_basin_area} is TRUE. If \code{stat} is one of 
#'   "min", "max", "mean" or "percent" the function assigns the value of the edge the site lies on. 
#'   Otherwise, the value is calculated as the sum of all edges upstream of the previous
#'   confluence and the proportional value of the edge the site lies on (based on
#'   the distance ratio 'ratio'); this is useful e.g. for counts of dams or waste water
#'   treatment plants or total catchment area.
#'   
#'   \code{input_attr_name} must give the column names of the edges attribute table
#'   for that the statistics should be calculated. 
#'
#' @note \code{\link{import_data}}, \code{\link{derive_streams}},
#'   \code{\link{calc_edges}}, \code{\link{calc_sites}} or
#'   \code{\link{calc_prediction_sites}} and \code{\link{calc_attributes_edges}}
#'   must be run before.
#'
#' @author Mira Kattwinkel, \email{mira.kattwinkel@@gmx.net}
#' @export
#' 
#' @examples 
#' \donttest{
#' # Initiate and setup GRASS
#' dem_path <- system.file("extdata", "nc", "elev_ned_30m.tif", package = "openSTARS")
#' if(.Platform$OS.type == "windows"){
#'   grass_program_path = "c:/Program Files/GRASS GIS 7.6"
#'   } else {
#'   grass_program_path = "/usr/lib/grass78/"
#'   }
#' 
#' setup_grass_environment(dem = dem_path, 
#'                         gisBase = grass_program_path,      
#'                         remove_GISRC = TRUE,
#'                         override = TRUE
#'                         )
#' gmeta()
#'                         
#' # Load files into GRASS
#' sites_path <- system.file("extdata", "nc", "sites_nc.shp", package = "openSTARS")
#' preds_v_path <- system.file("extdata", "nc", "geology.shp", package = "openSTARS")
#' import_data(dem = dem_path, sites = sites_path, predictor_vector = preds_v_path,)
#'
#' # Derive streams from DEM
#' derive_streams(burn = 0, accum_threshold = 700, condition = TRUE, clean = TRUE)
#'
#' # Check and correct complex confluences (there are no complex confluences in this
#' # example date set; set accum_threshold in derive_streams to a smaller value
#' # to create complex confluences)
#' cj <- check_compl_confluences()
#' if(cj){
#'   correct_compl_confluences()
#' }
#' 
#' # Prepare edges
#' calc_edges()
#' 
#' # Derive slope from the DEM as an example raster map to calculate attributes from
#' execGRASS("r.slope.aspect", flags = c("overwrite","quiet"),
#' parameters = list(
#'   elevation = "dem",
#'     slope = "slope"
#'     ))
#' calc_attributes_edges(input_raster = "slope", stat_rast = "max", attr_name_rast = "maxSlo",
#'                      input_vector = "geology", stat_vect = "percent", attr_name_vect = "GEO_NAME")
#' 
#' calc_sites() 
#'  
#' # approximate potential predictor variables for each site based on edge values
#' calc_attributes_sites_approx(
#' input_attr_name = c('maxSlo', 'CZamp', 'CZbgp', 'CZfgp', 'CZgp', 'CZigp', 'CZlgp', 'CZvep', 'Kmp'), 
#'   stat = c("max", rep("percent", 8)))
#' 
#' # plot share of a certain geology in the sampling point's catchment as
#' # point size
#' library(sp)
#' edges <- readVECT('edges', ignore.stderr = TRUE)
#' sites <- readVECT('sites', ignore.stderr = TRUE)
#' geo <- readVECT("geology", ignore.stderr = TRUE)
#' plot(geo, col = adjustcolor(1:8, alpha.f = 0.5)[as.factor(geo$GEO_NAME)])
#' plot(edges, col = "blue", add = TRUE)
#' plot(sites, col = 1, add = TRUE, pch = 19, cex = (sites$CZbgp + 0.15) * 2)
#' legend("left", col = adjustcolor(1:8, alpha.f = 0.5), bty = "n",
#' legend = unique(geo$GEO_NAME), pch = 15, title = "geology")
#' legend("right", col = 1, pch = 19, legend = seq(0, 1, 0.2), bty = "n",
#' title = "share CZbg\nin catchment", pt.cex =  (seq(0, 1, 0.2) + 0.15) * 2)
#'}
#'
calc_attributes_sites_approx <- function(sites_map = "sites",
                                         input_attr_name = NULL,
                                         output_attr_name = NULL,
                                         stat = NULL,
                                         round_dig = 2,
                                         calc_basin_area = TRUE,
                                         overwrite = FALSE){
  cnames <- execGRASS("db.columns", 
                      parameters = list(
                        table = sites_map), 
                      intern = TRUE)
  
  if(is.null(output_attr_name))
    output_attr_name <- input_attr_name
  if(calc_basin_area == TRUE){
    output_attr_name <- c("H2OAreaA", output_attr_name)
    input_attr_name <- c("H2OArea", input_attr_name)
    stat <- c("totalArea", stat)
    
    if("H2OAreaA" %in% cnames){
      execGRASS("v.db.dropcolumn", flags = "quiet",
                parameters = list(
                  map = sites_map,
                  columns = "H2OAreaA"
                ))
    }
  }
  if(calc_basin_area == FALSE & is.null(input_attr_name)){
    stop("Either input attribute name(s) must be given or calc_basin_area must be TRUE.")
  }
  
  if(any(output_attr_name %in% cnames) & overwrite == FALSE){
    stop(paste0("Column(s) ", paste0(output_attr_name[which(output_attr_name %in% cnames)], collapse = ", "), " already exist; use overwrite = TRUE to overwrite them."))
  }
  
  if(any(output_attr_name %in% cnames) & overwrite == TRUE){
    ii <- which(output_attr_name %in% cnames)
    jj <- which(output_attr_name[ii] == "H2OAreaA")
    if(length(jj) > 0)
      ii <- ii[-jj]
    for(i in ii){
      execGRASS("v.db.dropcolumn", flags = "quiet",
                parameters = list(
                  map = sites_map,
                  columns = output_attr_name[i]
                ))
    }
  }
  
  if(length(input_attr_name) != length(output_attr_name))
    stop("There must be the same number of input and output attribute names.")
  
  if(any(nchar(output_attr_name)) > 10)
    stop("Attribute names must not be longer than ten characters.")
  
  if(length(round_dig) == 1)
    round_dig <- rep(round_dig, length(output_attr_name))
  
  if(length(round_dig) < length(output_attr_name))
    round_dig <- c(max(round_dig), round_dig)
  
  if(length(unique(c(length(input_attr_name), length(stat),length(output_attr_name)))) > 1)
    stop(paste0("There must be the same number of input attribute names (",length(input_attr_name), "),
                output attribute names (", length(output_attr_name), ") and
                statistics to calculate  (", length(stat),")."))
  
  cnames <- execGRASS("db.columns", 
                      parameters = list(
                        table = sites_map), 
                      intern = TRUE)
  
  execGRASS("v.db.addcolumn",
            flags = c("quiet"),
            parameters = list(
              map = sites_map,
              columns = paste0(output_attr_name," double precision", collapse = ", ")
            ))
  
  for(i in seq_along(input_attr_name)){
    if(stat[i] %in% c("min", "max", "mean", "percent")){
      execGRASS("db.execute",
                parameters = list(
                  sql = paste0("UPDATE ", sites_map," SET ", output_attr_name[i], "=",
                               "(SELECT ", paste0(input_attr_name[i],"_c"),
                               " FROM edges WHERE edges.stream=", sites_map,".str_edge)")
                ))
    } else {
      # calculate site attribute as attribute of the two previous edges +
      # (1-ratio) * contribution of edge to total edge attribute
      # for H2O Area or e.g. for total numbers (no of WWTP per catchment)
      # e.g. calculated with stat = "sum" in calc_attributes_edges
      stream_prev1 <-  paste0("(SELECT prev_str01 FROM edges WHERE edges.stream=",sites_map,".str_edge)")
      stream_prev2 <-  paste0("(SELECT prev_str02 FROM edges WHERE edges.stream=",sites_map,".str_edge)")
      if(input_attr_name[i] == "H2OArea"){
        sql_str <-paste0("UPDATE ", sites_map," SET ",output_attr_name[i],
                         " = ROUND(((1-ratio)*",
                         "(SELECT rcaArea FROM edges WHERE ", sites_map,".str_edge = edges.stream) +",
                         "(SELECT H2OArea FROM edges WHERE edges.stream=",stream_prev1,") +",
                         "(SELECT H2OArea FROM edges WHERE edges.stream=",stream_prev2,")),",round_dig[i],")")
      } else {
        sql_str <-paste0("UPDATE ", sites_map," SET ",output_attr_name[i],
                         " = ROUND(((1-ratio)*",
                         "(SELECT ", paste0(input_attr_name[i],"_e"), " FROM edges WHERE ", sites_map,".str_edge = edges.stream) +",
                         "(SELECT ", paste0(input_attr_name[i],"_c"), " FROM edges WHERE edges.stream=",stream_prev1,") +",
                         "(SELECT ", paste0(input_attr_name[i],"_c"), " FROM edges WHERE edges.stream=",stream_prev2,")),",round_dig[i],")")
      }
      execGRASS("db.execute",
                parameters = list(
                  sql = sql_str
                ))
      # correct for those segments that do not have previous streams
      if(input_attr_name[i] == "H2OArea"){
        sql_str <- paste0("UPDATE ", sites_map," SET ",output_attr_name[i],
                          " = (1-ratio)*(SELECT rcaArea FROM edges WHERE ",
                          sites_map,".str_edge = edges.stream) WHERE str_edge IN ",
                          "(SELECT stream FROM edges WHERE prev_str01=0)")
      } else {
        sql_str <- paste0("UPDATE ", sites_map," SET ",output_attr_name[i],
                          " = (1-ratio)*(SELECT ", paste0(input_attr_name[i],"_e"),
                          " FROM edges WHERE ", sites_map,".str_edge = edges.stream) WHERE str_edge IN ",
                          "(SELECT stream FROM edges WHERE prev_str01=0)")
      }
      execGRASS("db.execute",
                parameters = list(
                  sql = sql_str
                ))
      # ROUND does not work with WHERE ... IN ...
      execGRASS("db.execute",
                parameters = list(
                  sql = paste0("UPDATE ",sites_map, " SET ",output_attr_name[i],
                               "= ROUND(",output_attr_name[i],",",round_dig[i],")")
                ))
    }
  }
  cnames2 <- execGRASS("db.columns", flags = "quiet",
                       parameters = list(
                         table = sites_map
                       ), intern = TRUE)
  cnames2 <- cnames2[-(which(cnames2 %in% cnames))]
  #message(writeLines(strwrap(paste0("\nNew attributes values are stored as ", paste("'", cnames2, "'", sep = "", collapse = ", "), " in 'sites'."),
  #        width = 80)))
  message(paste0("\nNew attributes values are stored as ", paste("'", cnames2, "'", sep = "", collapse = ", "), " in 'sites'."))
  
}