# NOT RUN {
# 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)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab