# 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
dem_path <- system.file("extdata", "nc", "elev_ned_30m.tif", package = "openSTARS")
sites_path <- system.file("extdata", "nc", "sites_nc.shp", package = "openSTARS")
streams_path <- system.file("extdata", "nc", "streams.shp", package = "openSTARS")
preds_v_path <- system.file("extdata", "nc", "pointsources.shp", package = "openSTARS")
preds_r_path <- system.file("extdata", "nc", "landuse_r.tif", package = "openSTARS")
                 
import_data(dem = dem_path, sites = sites_path, streams = streams_path,
            predictor_vector = preds_v_path, predictor_raster = preds_r_path)
# Derive streams from DEM
# burn in 'streams' 10 meters
derive_streams(burn = 10, 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()
}
# calculate slope as potential predictor
execGRASS("r.slope.aspect", flags = c("overwrite","quiet"),
parameters = list(
  elevation = "dem",
    slope = "slope"
    )) 
# Prepare edges
calc_edges()
calc_attributes_edges(input_raster = c("slope", "landuse_r"), 
                      stat_rast = c("max", "percent"), 
                      attr_name_rast = c("maxSlo", "luse"),
                      input_vector = "pointsources", stat_vect = "count",
                      attr_name_vect = "psource")
                      
# Plot eges with percentage of forest in the catchment (lusep_5) as line width
edges <- readVECT('edges', ignore.stderr = TRUE)
head(edges@data)
lu <- readRAST("landuse_r", ignore.stderr = TRUE, plugin = FALSE)
# plot landuse data
library(raster)
op <- par()
par(xpd = FALSE)
plot(raster(lu), legend = FALSE, xaxt = "n", yaxt = "n", bty = "n",
col = adjustcolor(c("red", "goldenrod", "green", "forestgreen",
"darkgreen", "blue", "lightblue"), alpha.f = 0.7))
par(xpd = TRUE)
legend("bottom", cex = 0.75,
  legend = c("developed", "agriculture", "herbaceous", "shrubland", "forest", "water", "sediment"),
  fill = c("red", "goldenrod", "green", "forestgreen","darkgreen", "blue", "lightblue"),
  horiz = TRUE, inset = -0.175)
  # line width is relative to the area of land use class 5 (forest) in the rca of the edge segment
plot(edges, lwd = edges$lusep_5_c * 10, add = TRUE)    
par <- op
# }
Run the code above in your browser using DataLab