Nothing
#'Calculate prediction sites for 'SSN' object.
#'
#'@description A vector (points) map of prediction sites is created and several
#'attributes are assigned.
#'
#'@param predictions string giving the name for the prediction sites map.
#'@param dist number giving the distance between the points to create in map
#' units.
#'@param nsites integer giving the approximate number of sites to create
#'@param netIDs integer (optional): create prediction sites only on streams with
#' these netID(s).
#'
#'@details Either \code{dist} or \code{nsites} must be provided. If \code{dist}
#'is NULL, it is estimated by dividing the total stream length in the map by
#'\code{nsites}; the number of sites actually derived might therefore be a bit
#'smaller than \code{nsites}.
#'
#'Steps include:
#'\itemize{
#'\item{Place points on edges with given distance from each other}
#'\item{Save the point coordinates in NEAR_X and NEAR_Y.}
#'\item{Assign unique identifiers (needed by the 'SSN' package) 'pid'
#'and 'locID'.}
#'\item{Get 'rid' and 'netID' of the stream segment the site
#'intersects with (from map 'edges').}
#'\item{Calculate upstream distance for
#'each point ('upDist').}
#'\item{Calculate distance ratio ('distRatio') between
#'position of the site on the edge (= distance traveled from lower end of the
#'edge to the site) and the total length of the edge.} }
#'
#''pid' and 'locID' are identical, unique numbers. 'upDist' is calculated using
#'\href{https://grass.osgeo.org/grass78/manuals/addons/r.stream.distance.html}{r.stream.distance}.
#'Points are created using
#'\href{https://grass.osgeo.org/grass78/manuals/v.segment.html}{v.segment}.
#'
#'@note \code{\link{import_data}}, \code{\link{derive_streams}} and
#' \code{\link{calc_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
#' 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)
#'
#' check_compl_confluences()
#' calc_edges()
#' calc_sites()
#' calc_prediction_sites(predictions = "preds", dist = 2500)
#'
#' library(sp)
#' dem <- readRAST('dem', ignore.stderr = TRUE, plugin = FALSE)
#' sites <- readVECT('sites', ignore.stderr = TRUE)
#' preds <- readVECT('preds', ignore.stderr = TRUE)
#' edges <- readVECT('edges', ignore.stderr = TRUE)
#' plot(dem, col = terrain.colors(20))
#' lines(edges, col = 'blue', lwd = 2)
#' points(sites, pch = 4)
#' points(preds, pch = 19, col = "darkred")
#' }
calc_prediction_sites <- function(predictions, dist = NULL, nsites = 10,
netIDs = NULL) {
# MiKatt 20200717
# WARNING: Values in column <cat> will be overwritten
vect <- execGRASS("g.list",
parameters = list(
type = "vect"
),
intern = TRUE)
if (!"edges" %in% vect)
stop("Edges not found. Did you run calc_edges()?")
if(predictions %in% vect)
execGRASS("g.remove",
flags = c("quiet", "f"),
parameters = list(
type = "vector",
name = predictions
))
if(all(is.null(c(dist, nsites))))
stop("Either the distance between prediction sites (dist) or the number of
prediction sites (nsites) must be given.")
temp_dir <- tempdir()
dt.streams <- do.call(rbind,strsplit(
execGRASS("db.select",
parameters = list(
sql = "select cat, stream, next_str, prev_str01,prev_str02,netID,Length from edges"
),intern = T),
split = '\\|'))
colnames(dt.streams) <- dt.streams[1,]
dt.streams <- data.table(dt.streams[-1,,drop = FALSE])
dt.streams[, names(dt.streams) := lapply(.SD, as.numeric)]
dt.streams[, offset := 0]
# omit all segements that do not belong to the netIDs given
if(!is.null(netIDs)){
dt.streams[!(netID %in% netIDs), offset := NA]
dt.streams <- stats::na.omit(dt.streams, cols = "offset")
}
if(nrow(dt.streams) == 0)
stop("No streams to place prediction points on. Please check netIDs.")
if(is.null(dist))
dist <- ceiling(sum(dt.streams[,Length]) / nsites)
message("Calculating point positions ...")
outlets <- dt.streams[next_str == -1, stream]
for(i in outlets){
calc_offset(dt.streams, id=i, offs = 0, dist)
}
pt <- 1
str1 <- NULL
for(i in 1:nrow(dt.streams)){
offs <- dt.streams[i, "offset", with = FALSE]
while(offs > 0){
str1 <- paste0(str1, "\n", paste("P", pt, dt.streams[i, "cat", with = FALSE], offs, sep=" "))
offs <- offs - dist
pt <- pt + 1
}
}
str1 <- substring(str1, 2)
write(str1, file.path(temp_dir,"pt.txt"))
execGRASS("v.segment", flags = c("overwrite", "quiet"),
parameters = list(
input = "edges",
output = predictions,
rules = file.path(temp_dir,"pt.txt")
))
# MiKatt: No line break in long strings on Windows!
message("Creating attribute table ...")
execGRASS("v.db.addtable", flags = c("quiet"),
parameters = list(
map = predictions,
columns = "cat_edge int,str_edge int,dist double precision,nx double precision,ny double precision,pid int,loc int,net int,rid int,out_dist double,distalong double precision,ratio double precision"
), ignore.stderr = TRUE, intern = TRUE)
# MiKatt: Necessary to get upper and lower case column names
execGRASS("v.db.renamecolumn", flags = "quiet",
parameters = list(
map = predictions,
column = "loc,locID"
))
execGRASS("v.db.renamecolumn", flags = "quiet",
parameters = list(
map = predictions,
column = "net,netID"
))
execGRASS("v.db.renamecolumn", flags = "quiet",
parameters = list(
map = predictions,
column = "out_dist,upDist"
))
execGRASS("v.db.renamecolumn", flags = "quiet",
parameters = list(
map = predictions,
column = "nx,NEAR_X"
))
execGRASS("v.db.renamecolumn", flags = "quiet",
parameters = list(
map = predictions,
column = "ny,NEAR_Y"
))
message("Setting cat_edge ...")
# MiKatt: additionally get x and y coordinate
execGRASS("v.distance",
flags = c("overwrite", "quiet"),
parameters = list(from = predictions,
to = "edges",
upload = "cat,dist,to_x,to_y",
column = "cat_edge,dist,NEAR_X,NEAR_Y"))
message("Setting pid and locID ...")
execGRASS("v.db.update",
parameters = list(map = predictions,
column = "pid",
value = "cat"))
execGRASS("v.db.update",
parameters = list(map = predictions,
column = "locID",
value = "pid"))
# Set netID and rid from network ---------
message("Assigning netID and rid ...")
sql_str<- paste0("UPDATE ", predictions, " SET rid=(SELECT rid FROM edges WHERE ",
predictions,".cat_edge=edges.cat)")
execGRASS("db.execute",
parameters = list(
sql = sql_str
))
sql_str<- paste0("UPDATE ", predictions, " SET netID=(SELECT netID FROM edges WHERE ",
predictions,".cat_edge=edges.cat)")
execGRASS("db.execute",
parameters = list(
sql = sql_str
))
sql_str<- paste0("UPDATE ", predictions, " SET str_edge=(SELECT stream FROM edges WHERE ",
predictions,".cat_edge=edges.cat)")
execGRASS("db.execute",
parameters = list(
sql = sql_str
))
# Calculate upDist ---------
message("Calculating upDist ...")
## MiKatt was not exact enough, results in identical upDist if two points lay
## in the same raster cell
# execGRASS("r.stream.distance",
# flags = c("overwrite", "quiet", "o"),
# parameters = list(
# stream_rast = "streams_r",
# direction = "dirs",
# method = "downstream",
# distance = "upDist"
# ))
#
# execGRASS("v.what.rast",
# flags = c("quiet"),
# parameters = list(
# map = predictions,
# raster = "upDist",
# column = "upDist"
# ))
# MiKatt: ! Round upDist to m
execGRASS("v.distance", flags = c("quiet"),
parameters =list(
from = predictions,
to = "edges",
to_type = "line",
upload = "to_along",
column = "distalong"
))
sql_str <- paste0('UPDATE ', predictions, ' SET upDist=',
'round(((SELECT upDist FROM edges WHERE edges.cat=',
predictions, '.cat_edge)-distalong),2)')
execGRASS("db.execute",
parameters = list(
sql=sql_str
))
# Calculate distRatio = distance from lower end of edge to site / length edge
message("Calculating distance ratio ...")
sql_str <- paste0('UPDATE ', predictions, ' SET ratio=1-',
'distalong/',
'(SELECT Length FROM edges WHERE edges.cat=', predictions, '.cat_edge)')
execGRASS("db.execute",
parameters = list(
sql=sql_str
))
execGRASS("v.db.dropcolumn",
map = predictions,
columns = "cat_edge")
}
#' Calculate offset
#' @description Recursive function to calculate the offset from the downstream
#' junction needed to place points with fixed distance along the streams.
#' It is called by \code{\link{calc_prediction_sites}} for each
#' outlet and should not be called by the user.
#'
#' @param dt data.table containing the attributes of the stream segments
#' @param id integer; 'stream' of the stream segment
#' @param offs number; offset from outlet of the stream segment (downstream);
#' equals the length of the segment if the point shall be placed directly at the
#' downstream junction.
#' @param dist number giving the distance between the points to create in map units.
#' @keywords internal
#'
#' @return Nothing; change 'offset' in dt.
#'
#' @author Mira Kattwinkel, \email{mira.kattwinkel@@gmx.net}
#'
#'@examples
#'\dontrun{
#' outlets <- dt.streams[next_str == -1, stream]
#' netID <- 1
#' for(i in outlets){
#' calc_offset(dt.streams, id = i, offs = 0, dist = 200)
#' }
#'}
calc_offset <- function(dt, id, offs, dist){
if(dt[stream == id, prev_str01,] == 0){ # check only one of prev01 and prev02 because they are always both 0
dt[stream == id, offset := floor(Length - offs)]
} else {
dt[stream == id, offset := floor(Length - offs)]
if(offs < 0){
offs <- offs + dt[stream == id, Length]
} else {
offs <- dist - (dt[stream == id, Length - offs] %% dist)
}
calc_offset(dt, dt[stream == id, prev_str01], offs, dist)
calc_offset(dt, dt[stream == id, prev_str02] ,offs, dist)
}
}
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.