Nothing
#' Calculate attributes of the sites.
#'
#' For each site (observation or prediction) the total catchment area is
#' calculated ('H2OArea'). Additionally, other attributes (predictor variables)
#' can be derived based on given raster or vector maps. This function calculates
#' exact values for catchments derived with
#' \href{https://grass.osgeo.org/grass78/manuals/addons/r.stream.basins.html}{r.stream.basins}
#' and can take considerable time if there are many sites.
#' Catchment raster maps can optionally be stored as "sitename_catchm_X" (X = locID).
#' @import progress
#'
#' @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_raster character vector (optional); name of additional raster
#' maps to calculate attributes from.
#' @param stat_rast character vector (optional); statistics to be calculated, one of:
#' "min", "max", "mean", "stddev", "variance", "sum", "median", "percent", "area" or "percentile_X" (where X
#' gives the desired percentile e.g. 25 for the first). Must be provided if
#' \code{input_raster} are given.
#' @param attr_name_rast character vector (optional); column name for the attributes
#' to be calculated. Attribute names must not be longer than 10 characters.
#' Must be provided if \code{input_raster} are given.
#' @param input_vector character string vector (optional); name of additional vector
#' maps to calculate attributes from.
#' @param stat_vect character string vector (optional); statistics to be calculated,
#' one of: "percent", "area" or "count." Must be provided if \code{input_vector} is given.
#' @param attr_name_vect character string vector (optional); column name(s) in
#' the vector file provided to calculate the attributes from (if \code{input_vector}
#' is a polygon map and stat_vect is "percent") or giving the new name attributes
#' to calculate (if \code{input_vector} is a point map and stat_vect is "count".
#' Must be provided if \code{input_vector} is given.
#' @param round_dig integer; number of digits to round results to. Can be a vector
#' of different values or just one value for all attributes.
#' @param calc_basin_area boolean; shall the catchment area be calculated? (Useful
#' to set to FALSE if the function has been called before with \code{keep_basins = TRUE}.)
#' @param keep_basins boolean; shall raster and vector maps of all the watersheds be kept?
#' Defaults to FALSE.
#' @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{'H2OArea':} {Total watershed area of the watershed upstream of each site.}
#' \item{attr_name_rast:} {Additional optional attributes calculated based on input_raster maps.}
#' \item{attributes form vector maps:}{Additional optional attributes
#' calculated based on input_vector maps. The column names are based on the unique entries
#' of the column(s) given in \code{attr_name_vect}.}
#' }
#' Please note that for sampling points that lie in the same DEM raster cell
#' along a stream identical values are calculated because identical watersheds
#' are derived.
#'
#' @note \code{\link{import_data}}, \code{\link{derive_streams}},
#' \code{\link{calc_edges}} and \code{\link{calc_sites}} or
#' \code{\link{calc_prediction_sites}} must be run before.
#'
#' If \code{calc_basin_area = F} but there are no raster maps called 'sitename_catchm_x'
#' with x = locID of all sites the catchments (and their area) are derived.
#'
#' @author Mira Kattwinkel, \email{mira.kattwinkel@@gmx.net},
#' Eduard Szoecs, \email{eduardszoecs@@gmail.com}
#' @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
#' dem_path <- system.file("extdata", "nc", "elev_ned_30m.tif", package = "openSTARS")
#' sites_path <- system.file("extdata", "nc", "sites_nc.shp", package = "openSTARS")
#' import_data(dem = dem_path, sites = sites_path)
#'
#' # Derive streams from DEM
#' derive_streams(burn = 0, accum_threshold = 700, condition = TRUE, clean = TRUE)
#'
#' # Prepare edges
#' calc_edges()
#'
#' # caluclate slope as potential predictor
#' 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")
#'
#' # Prepare sites
#' calc_sites()
#' calc_attributes_sites_approx(input_attr_name = "maxSlo", output_attr_name = "maxSloA", stat = "max")
#' calc_attributes_sites_exact(input_raster = "slope", attr_name_rast = "maxSloE", stat_rast = "max")
#'
#' # Plot data
#' library(sp)
#' dem <- readRAST('dem', ignore.stderr = TRUE, plugin = FALSE)
#' edges <- readVECT('edges', ignore.stderr = TRUE)
#' sites <- readVECT('sites', ignore.stderr = TRUE)
#' plot(dem, col = gray(seq(0,1,length.out=20)))
#' mm <- range(c(edges$maxSlo_e, sites$maxSloA, sites$maxSloE))
#' b <- seq(from = mm[1], to = mm[2] + diff(mm) * 0.01, length.out = 10)
#' c_ramp <- colorRampPalette(c("white", "blue", "orange", "red"))
#' cols <- c_ramp(length(b))[as.numeric(cut(edges$maxSlo_e, breaks = b, right = FALSE))]
#' # plot stream edges, color depending on maxSlope of the edge
#' plot(edges,col = cols, lwd = 2, add = TRUE)
#' cols <- c_ramp(length(b))[as.numeric(cut(sites$maxSloA,breaks = b,right = FALSE))]
#' # plot sites as points with color corresponding to maxSlop approximate
#' plot(sites, pch = 19, col = cols, cex = 2, add = TRUE)
#' cols <- c_ramp(length(b))[as.numeric(cut(sites$maxSloE,breaks = b,right = FALSE))]
#' #' # plot sites as ring around points with color corresponding to maxSlop exact
#' plot(sites, pch = 21, bg = cols, cex = 1.1, add = TRUE)
#' # Some points in the lower centre of the map indicate a difference in max slope between
#' # approximate and exact calculation (different colors for inner and outer points). However,
#' # for most points the values are similar
#' }
calc_attributes_sites_exact <- function(sites_map = "sites",
input_raster = NULL,
stat_rast = NULL,
attr_name_rast = NULL,
input_vector = NULL,
stat_vect = NULL,
attr_name_vect = NULL,
round_dig = 2,
calc_basin_area = TRUE,
keep_basins = FALSE,
overwrite = FALSE){
if(length(input_raster) != length(stat_rast) | length(input_raster) != length(attr_name_rast) | length(attr_name_rast) != length(stat_rast))
stop(paste0("There must be the same number of input raster files (",length(input_raster), "), statistics to calculate (",
length(stat_rast), ") and attribute names (", length(attr_name_rast),")."))
if(length(input_vector) != length(stat_vect))
stop(paste0("There must be the same number of input vector files (",length(input_vector), ")
and statistics to calculate (", length(stat_vect), ")."))
if(!is.null(stat_rast) & any(stat_rast %in% c("min","max", "mean", "stddev","variance","sum","median", "percent","area") + grepl("percentile", stat_rast)) == 0) # TRUE = 1, FALSE = 0
stop('Statistisc to calculate from raster data must be one of "min","max", "mean", "stddev","variance","sum", "median", "percentile_X", "percent" or "area".')
if(any(!stat_vect %in% c("percent", "count", "area")))
stop('statistics to calculate from vector data must be one of "count", "percent" or "area".')
# 1 for area, 2 for points
vtype <- rep(1, length(stat_vect))
if(!is.null(stat_vect)){
for(i in 1:length(stat_vect)){
a <- execGRASS("v.info", flags = "t",
parameters = list(
map = input_vector[i]
), intern = TRUE)
a <- do.call(rbind,strsplit(a, "="))
k <- which(a[,1] == "points")
if(as.numeric(a[k,2]) != 0){
vtype[i] <- 2
if(stat_vect[i] != "count"){
stop('If an input vector is of type point the statistic to calculate must be "count".')
}
}
}
}
if(is.null(stat_rast) & is.null(stat_vect) & !calc_basin_area)
stop("Either the catchment areas are calculated or a statistic to calculate must be provided.")
temp_dir <- tempdir()
if(length(round_dig) == 1)
round_dig <- rep(round_dig, length(stat_rast)+1)
if(length(round_dig) == length(stat_rast))
round_dig <- c(round_dig[1], round_dig)
rast <- execGRASS("g.list",
parameters = list(
type = "rast"
),
intern = TRUE)
if ("MASK" %in% rast)
execGRASS("r.mask",flags = c("r", "quiet"))
if(!all(input_raster %in% rast)){
if(length(input_raster)>1){
i <- which(input_raster %in% rast)
mes <- input_raster[-i]
} else {
mes <- input_raster
}
stop(paste0("Missing input raster data ", paste0("'", mes, "'", collapse = ", "),
". Please give valid raster names. \nAvailable raster are ",
paste0("'", rast, "'", collapse = ", ")))
}
vect <- execGRASS("g.list",
parameters = list(
type = "vector"
),
intern = TRUE)
if(!all(input_vector %in% vect)){
if(length(input_vector)>1){
i <- which(input_vector %in% vect)
mes <- input_vector[-i]
} else {
mes <- input_vector
}
stop(paste0("Missing input vector data ", paste0("'", mes, "'", collapse = ", "),
". Please give valid vector file names. \nAvailable vector files are ",
paste0("'", vect, "'", collapse = ", ")))
}
d.sites <- readVECT(sites_map, ignore.stderr = TRUE)
if(!all(paste0(sites_map, "_catchm_", d.sites@data$locID) %in% rast)){
calc_basin_area <- TRUE
}
if(any(d.sites@data$ratio == 0) & calc_basin_area){
d.edges <- readVECT("edges", ignore.stderr = TRUE)
dt.edges <- setDT(d.edges@data)
dt.edges[, colnames(dt.edges)[-which(colnames(dt.edges) %in% c("cat", "stream","prev_str01","prev_str02","rid","H2OArea"))] := NULL]
rm(d.edges)
}
locID <- unique(d.sites@data$locID)
dat.h2o <- dat.rast <- dat.vect <- NULL
if(calc_basin_area){
dat.h2o <- matrix(nrow = length(locID), ncol = 1, data = 0)
colnames(dat.h2o) <- c("H2OArea")
}
if(!is.null(input_raster)){
vals.rast <- lapply(1:length(input_raster), function(x) ifelse(stat_rast[x] == "percent", list(get_all_raster_values(input_raster[x])), 1))
ndat <- length(unlist(vals.rast))
dat.rast <- matrix(nrow = length(locID), ncol = ndat, data = 0)
cnames.rast <- lapply(1:length(attr_name_rast), function(x) ifelse(stat_rast[x] %in% c("percent", "area"),
list(paste(attr_name_rast[x], unlist(vals.rast[[x]]), sep = "_")),
attr_name_rast[x]))
colnames(dat.rast) <- unlist(cnames.rast)
}
if(!is.null(input_vector)){
attribute_cats <- NULL
for(i in 1:length(input_vector)){
if(vtype[i] == 1){
cnames <- execGRASS("db.columns", flags = "quiet",
parameters = list(
table = input_vector[i]
), intern = TRUE)
if(!attr_name_vect[i] %in% cnames)
stop(paste0("Invalid vector attribute name ", paste0("'", attr_name_vect[i], "'", collapse = ", "),
". Please give a valid column name. \nAvailable columns are ",
paste0("'", cnames, "'", collapse = ", ")))
attribute_cats <- c(attribute_cats,
paste0(unique(execGRASS("db.select", flags = c("c"),
parameters = list(
sql = paste0("select ", attr_name_vect[i], " from ",input_vector[i])
), ignore.stderr = TRUE, intern = TRUE)), substr(stat_vect[i], 1, 1))
)
} else {
attribute_cats <- c(attribute_cats, attr_name_vect[i])
}
if(length(round_dig) < length(attribute_cats) + length(input_raster) + calc_basin_area)
round_dig <- c(round_dig, rep(round_dig[1], length(attribute_cats)))
}
dat.vect <- matrix(ncol = length(attribute_cats), nrow = length(locID), data = 0)
colnames(dat.vect) <- c(attribute_cats)
}
cnames_new <- c()
if(exists("dat.h2o"))
cnames_new <- c(cnames_new, colnames(dat.h2o))
if(exists("dat.rast"))
cnames_new <- c(cnames_new, colnames(dat.rast))
if(exists("dat.vect"))
cnames_new <- c(cnames_new, colnames(dat.vect))
cellsize <- NULL
if(any(stat_rast == "area")){
cellsize <- execGRASS("g.region", flags = "p",intern=T)
cellsize <- as.numeric(do.call(rbind,strsplit(cellsize[grep("res",cellsize)],split=":"))[,2])
cellsize <- prod(cellsize)
}
cnames_sites <- execGRASS("db.columns", flags = "quiet",
parameters = list(
table = "sites"
), intern = TRUE)
if(any(cnames_new %in% cnames_sites) & overwrite == FALSE){
stop(paste0("Column(s) ", paste0(cnames_new[which(cnames_new %in% cnames_sites)], collapse = ", "), " already exist; use overwrite = TRUE to overwrite them."))
}
message("Intersecting maps for ", nrow(d.sites@data)," sites ...")
# progress bar
pb <- progress::progress_bar$new(total = nrow(d.sites@data))
for (i in seq_along(locID)) {
#message(i)
id <- locID[i]
#dat[i,"locID"] <- locID
# get first entry from d.sites@data with locID
ii <- which(d.sites@data$locID == id)[1]
# calculate drainage area in km^2
if(calc_basin_area){
# If the distance ratio is 0, the site lies within the outflow cell of
# the edge; then, r.stream.basins will extract a too large basin including
# the second tributary of the confluence
if(d.sites@data$ratio[ii] == 0){
j <- which(dt.edges[, "rid"] == d.sites@data$rid[ii])
# use $ here to get single value from data.table to match dimensions of dat
dat.h2o[i,"H2OArea"] <- round(dt.edges$H2OArea[j], round_dig[1])
streams <- get_streams_edges_in_catchment(dt.edges, dt.edges$stream[j])
## new, simpler approach
m <- paste0("if(", paste0("rca == ", streams, collapse = " || "), ", 1, 0)")
execGRASS("r.mapcalc", flags = c("overwrite","quiet"),
parameters = list(
expression = paste0(sites_map, "_catchm_", id, " = ", m)
), intern = TRUE, ignore.stderr = TRUE)
## old approach
# ## more than 106 categories in r.mask crashes (i.e. no MASK is created)
# ## workaround:
# n <- length(cats) %/% 100
# if(n == 0){
# execGRASS("r.mask", flags = c("overwrite","quiet"),
# parameters = list(
# raster = "rca",
# maskcats = paste0(cats, collapse = " "))
# )
# execGRASS("r.mapcalc", flags = c("overwrite","quiet"),
# parameters = list(
# expression = paste0("'", paste0(sites_map, "_catchm_",locID), "=MASK'")
# ))
# execGRASS("r.mask", flags = c("r","quiet")
# )
# } else {
# count <- 0
# for(j in 1:n){
# execGRASS("r.mask", flags = c("overwrite","quiet"),
# parameters = list(
# raster = "rca",
# maskcats = paste0(cats[((j-1)*100+1):(j*100)], collapse = " "))
# )
# execGRASS("r.mapcalc", flags = c("overwrite","quiet"),
# parameters = list(
# expression = paste0("'temp", j, "=MASK'")
# ))
# execGRASS("r.mask", flags = c("r","quiet")
# )
# count <- count + 1
# }
# if(length(cats) %% 100 != 0){
# execGRASS("r.mask", flags = c("overwrite","quiet"),
# parameters = list(
# raster = "rca",
# maskcats = paste0(cats[(n*100+1):length(cats)], collapse = " "))
# )
# execGRASS("r.mapcalc", flags = c("overwrite","quiet"),
# parameters = list(
# expression = paste0("'temp", n+1, "=MASK'")
# ))
# execGRASS("r.mask", flags = c("r","quiet")
# )
# count <- count + 1
# }
# m <- paste0("temp",1:count)
# m <- paste0("if(", paste0(m, collapse = "|||"),",1)")
# execGRASS("r.mapcalc", flags = c("overwrite","quiet"),
# parameters = list(
# expression = paste0('\"', paste0(sites_map, '_catchm_',locID), '=', m,'\"')
# ))
# execGRASS("g.remove", flags = c("f","quiet"),
# parameters = list(
# type = "raster",
# name = paste0("temp",1:count)
# ))
#}
} else {
coord <- sp::coordinates(d.sites[ii,])
execGRASS("r.stream.basins",
flags = c("overwrite", "l", "quiet"),
parameters = list(direction = "dirs",
coordinates = coord,
basins = paste0(sites_map, "_catchm_",id)))
dat.h2o[i,"H2OArea"] <- round(as.numeric(as.character(strsplit(
execGRASS("r.stats",
flags = c("a", "quiet"),
parameters = list(input = paste0(sites_map, "_catchm_",id)),
intern = TRUE)[1], split = ' ')[[1]][[2]]))/1000000,round_dig[1])
}
}
##########################
# raster data
# calculate unviriate statistics per watershed
if(length(stat_rast) > 0){
# set mask to the current basin
execGRASS("r.mask",
flags = c("overwrite", "quiet"),
parameters = list(
raster = paste0(sites_map, "_catchm_",id)))
st <- execGRASS("r.univar",
flags = c("overwrite", "quiet", "t"),
parameters = list(
map = "MASK",
separator = ","),
ignore.stderr = TRUE,
intern = TRUE)
st <- data.frame(do.call(rbind,strsplit(st,",")), stringsAsFactors = FALSE)
non_null_cells <- as.numeric(st[2,which(st[1,] == "non_null_cells")])
for(j in 1:length(stat_rast)){
if(stat_rast[j] == "percent" | stat_rast[j] == "area"){
st <- execGRASS("r.stats", flags = c("c"),
parameters = list(
input = input_raster[j], #paste0("rca,", input_raster[j]),
separator = "comma"
),
ignore.stderr = TRUE,
intern = TRUE)
st <- data.frame(do.call(rbind,strsplit(st,",")), stringsAsFactors = FALSE)
colnames(st) <- c("class", "cellcount")
st$cellcount <- as.numeric(st$cellcount)
#suppressWarnings(st$zone <- as.numeric(st$zone))
setDT(st)
#st[, `:=`(names(st), lapply(.SD, as.numeric))]
#st <- dcast(st, "zone ~ class", value.var = "cellcount")
#st[, all_cells := rowSums(.SD, na.rm = TRUE), .SDcols = 2:ncol(st)]
#st[, "*" := NULL]
#st[, 2:(ncol(st)-1) := lapply(.SD, function(x) x / st$all_cells), .SDcols = 2:(ncol(st)-1)]
#st[, all_cells := NULL]
#colnames(st)[-1] <- paste(attr_name_rast[j], colnames(st)[-1], sep = "_")
#rca_cell_count <- merge(rca_cell_count, st, by = "zone", all.x = TRUE)
#stat_r <- c(stat_r, rep("percent_class", ncol(st) - 1))
if(length(ind <- which(st$class == "*")) > 0)
st <- st[-ind,, drop = FALSE]
if(nrow(st) > 0){
if(stat_rast[j] == "percent"){
dat.rast[i, paste(attr_name_rast[j], st$class, sep = "_")] <- st$cellcount / non_null_cells
} else { #"area
dat.rast[i, paste(attr_name_rast[j], st$class, sep = "_")] <- st$cellcount * cellsize
}
}
} else {
if(stat_rast[j] == "median"){
st <- execGRASS("r.univar",
flags = c("overwrite", "quiet","g", "e"),
parameters = list(
map = input_raster[j]), intern = TRUE)
} else if(grepl("percentile",stat_rast[j])){
p <- as.numeric(as.character(unlist(strsplit(stat_rast[j],"_"))[2]))
st <- execGRASS("r.univar",
flags = c("overwrite", "quiet","g", "e"),
parameters = list(
map = input_raster[j],
percentile = p), intern = TRUE)
} else {
st <- execGRASS("r.univar",
flags = c("overwrite", "quiet","g"),
parameters = list(
map = input_raster[j]), intern = TRUE)
}
st <- setDT(data.frame(do.call(rbind,strsplit(st,"="))))
if(nrow(st) > 0){
st[,X2 := as.numeric(as.character(X2))]
if(grepl("percent", stat_rast[j])){
if(st[X1 == "variance", X2] == 0){ # if coded as something and NA, null(), no data valu
dat.rast[i, cnames.rast[[j]]] <- round(st[X1 == "n", X2] /non_null_cells, round_dig[j + 1])
} else{ # if coded as 1 and 0, "mean" gives ratio
dat.rast[i, cnames.rast[[j]]] <- round(st[X1 == "mean", X2], round_dig[j + 1])
}
} else
dat.rast[i, cnames.rast[[j]]] <- round(st[X1 == stat_rast[j], X2], round_dig[j + 1])
} else
dat.rast[i, cnames.rast[[j]]] <- 0
}
}
# Remove the mask!
execGRASS("r.mask",
flags = c("r", "quiet"))
}
###################################################
# vector data
if(length(input_vector) > 0){
# convert raster catchment to vector
rname <- paste0(sites_map, "_catchm_",id)
vname <- paste0(sites_map, "_catchm_",id, "_v")
execGRASS("r.to.vect", flags = c("overwrite","v", "quiet"),
parameters = list(
input = rname,
output = vname,
type = "area"
))
# calculate area
## GRASS version below 7.8
## v.to.db needs the column to be pobulated to exist; from 7.8 onward this column is created and existing ones are not automatically overwritten
# if(compareVersion(strsplit(system2("grass", "--version", stdout = TRUE, stderr = TRUE)[1], " ")[[1]][3], "7.8") < 0){
# execGRASS("v.db.addcolumn",
# parameters = list(
# map = vname,
# columns = "area double precision"
# ))
# }
# execGRASS("v.to.db",flags = "quiet",
# parameters = list(
# map = vname,
# option = "area",
# columns = "area"
# ), ignore.stderr = TRUE)
grass_v.to.db(map = vname, option = "area", columns = "area", format = "double precision")
carea <- sum(as.numeric(execGRASS("v.db.select",flags = "quiet",
parameters = list(
map = vname,
columns = "area"
), intern = TRUE)[-1]))
j.count <- 1 + calc_basin_area * length(input_raster)
for(j in 1:length(input_vector)){
# if this is no point vector
if(vtype[j] == 1){
# intersect with catchment
execGRASS("v.overlay", flags = c("overwrite","quiet"),
parameters = list(
ainput = input_vector[j],
binput = vname,
operator = "and",
output = "intersect_out",
olayer = "1,0,0"
), ignore.stderr = TRUE, intern = TRUE)
# calculate area of all features
## GRASS version below 7.8
## v.to.db needs the column to be pobulated to exist; from 7.8 onward this column is created and existing ones are not automatically overwritten
# if(compareVersion(strsplit(system2("grass", "--version", stdout = TRUE, stderr = TRUE)[1], " ")[[1]][3], "7.8") < 0){
# execGRASS("v.db.addcolumn",
# parameters = list(
# map = "intersect_out",
# columns = "area double precision"
# ))
# }
# execGRASS("v.to.db",flags = "quiet",
# parameters = list(
# map = "intersect_out",
# option = "area",
# columns = "area"
# ), ignore.stderr = TRUE)
grass_v.to.db(map = "intersect_out", option = "area", columns = "area", format = "double precision")
# get the areas per value of the attribute
a <- execGRASS("db.select",flags = "c",
parameters = list(
sql = paste0("select a_", attr_name_vect[j],",sum(area) from intersect_out group by a_", attr_name_vect[j])
), intern = TRUE)
if(length(a) > 0){
a <- do.call(rbind,strsplit(a,split = '\\|'))
a <- data.frame(a, stringsAsFactors = F)
if(stat_vect[j] == "percent"){ # do not divide by area if calculating total area
a[,2] <- round(as.numeric(a[,2]) / carea, round_dig[j.count])
} else {
a[,2] <- round(as.numeric(a[,2]), round_dig[j.count])
}
a[,1] <- paste0(a[,1], substr(stat_vect[j], 1, 1))
dat.vect[i, a[,1]] <- a[,2]
}
} else { # if this is a point vector
dat.vect[i, attr_name_vect[j]] <- as.numeric(unlist(strsplit(
execGRASS("v.vect.stats", flags = c("p", "quiet"),
parameters = list(
points = input_vector[j],
areas = vname,
separator = ","
), intern = TRUE)[2], split = ","))[2])
}
j.count <- j.count + 1
}
}
###################################################
# Delete watershed raster
if (!keep_basins) {
rast <- execGRASS("g.list",
parameters = list(
type = "rast"
),
intern = TRUE)
if(paste0(sites_map, "_catchm_",id) %in% rast){
execGRASS("g.remove",
flags = c("quiet", "f"),
parameters = list(
type = "raster",
name = paste0(sites_map, "_catchm_",id)
))
}
vect <- execGRASS("g.list",
parameters = list(
type = "vector"
),
intern = TRUE)
if(paste0(sites_map, "_catchm_",id, "_v") %in% vect){
execGRASS("g.remove",
flags = c("quiet", "f"),
parameters = list(
type = "vector",
name = paste0(sites_map, "_catchm_",id, "_v")
), ignore.stderr = TRUE)
}
}
pb$tick()
}
# change column names if there is just one raster class (and NULL)
if(!is.null(input_raster)){
cnames.rast <- lapply(1:length(attr_name_rast), function(x) ifelse(stat_rast[x] %in% c("percent", "area") & length(unlist(vals.rast[[x]])) > 1,
list(paste(attr_name_rast[x], unlist(vals.rast[[x]]), sep = "_")),
attr_name_rast[x]))
colnames(dat.rast) <- unlist(cnames.rast)
}
dat <- cbind(dat.h2o, dat.rast, dat.vect, locID)
# truncate column names
n <- nchar(colnames(dat))
nn <- which(n > 10)
if(length(nn) > 0){
colnames(dat)[nn] <- paste0(substr(colnames(dat)[nn], start = 1, stop = 7), substr(colnames(dat)[nn], start = nchar(colnames(dat)[nn]), stop = nchar(colnames(dat)[nn])))
}
# Join attributes to sites attribute table
message("Joining new attributes to attribute table ...")
utils::write.csv(dat, file.path(temp_dir,"sites_attributes_exact.csv"),row.names = F)
write.table(t(gsub("numeric","Real",apply(dat,2,class))),file.path(temp_dir,"sites_attributes_exact.csvt"),quote=T,sep=",",row.names = F,col.names = F)
execGRASS("db.in.ogr", flags = c("overwrite","quiet"),
parameters = list(
input = file.path(temp_dir,"sites_attributes_exact.csv"),
output = "sites_attributes_exact"
),ignore.stderr = T)
execGRASS("v.db.join", flags = "quiet",
parameters = list(
map = sites_map,
column = "locID",
other_table = "sites_attributes_exact",
other_column = "locID"
))
execGRASS("db.droptable", flags = c("quiet","f"),
parameters = list(
table = "sites_attributes_exact"
))
invisible(file.remove(file.path(temp_dir,"sites_attributes_exact.csv")))
invisible(file.remove(file.path(temp_dir,"sites_attributes_exact.csvt")))
# cnames_sites2 <- execGRASS("db.columns", flags = "quiet",
# parameters = list(
# table = sites_map
# ), intern = TRUE)
# cnames_sites2 <- cnames_sites2[-(which(cnames_sites2 %in% cnames_sites))]
#message(writeLines(strwrap(paste0("\nNew attributes values are stored as ", paste("'", cnames_sites2, "'", sep = "", collapse = ", "), " in 'sites'."),
# width = 80)))
message(paste0("\nNew attributes values are stored as ", paste("'", cnames_new, "'", sep = "", collapse = ", "), " in 'sites'."))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.