Learn R Programming

openSTARS (version 1.1.0)

calc_attributes_sites_approx: Calculate attributes of the sites.

Description

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).

Usage

calc_attributes_sites_approx(sites_map = "sites", input_attr_name,
  output_attr_name = NULL, stat, round_dig = 2, calc_basin_area = TRUE)

Arguments

sites_map

character; name of the sites the attributes shall be calculated for. "sites" refers to the observation sites.

input_attr_name

character vector; input column name in the edges attribute table.

output_attr_name

character vector (optional); output column name appended to the site attribute data table. If not provided it is set to input_attr_name. Attribute names must not be longer than 10 characters.

stat

name or character vector giving the statistics to be calculated. See details below.

round_dig

integer; number of digits to round results to.

calc_basin_area

boolean; shall the catchment area be calculated? (Useful to set to FALSE if the function has been called before.)

Value

Nothing. The function appends new columns to the sites_map attribute table

  • 'H2OAreaA': Total watershed area of the watershed upstream of each site.

  • attr_name: Additional optional attributes calculated based on input_attr_name.

Details

The approximate total catchment area (H2OAreaA) is always calculated if calc_basin_area is TRUE. If 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 junction 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 plant or total catchment area.

Examples

Run this code
# NOT RUN {
# Initiate GRASS session
if(.Platform$OS.type == "windows"){
  gisbase = "c:/Program Files/GRASS GIS 7.4.0"
  } else {
  gisbase = "/usr/lib/grass74/"
  }
initGRASS(gisBase = gisbase,
    home = tempdir(),
    override = TRUE)

# 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")
pred_path <- system.file("extdata", "nc", "landuse.shp", package = "openSTARS")
setup_grass_environment(dem = dem_path)
import_data(dem = dem_path, sites = sites_path,
 predictor_vector = pred_path, predictor_v_names = "landuse")
gmeta()

# Derive streams from DEM
derive_streams(burn = 0, accum_threshold = 700, condition = TRUE, clean = TRUE)

# Check and correct complex junctions (there are no complex juctions in this 
# example date set)
cj <- check_compl_junctions()
if(cj){
  correct_compl_junctions()
}

# 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 = "landuse", stat_vect = "percent", attr_name_vect = "landuse")

calc_sites() 
 
# approximate potential predictor variables for each site based on edge values
calc_attributes_sites_approx(input_attr_name = c("maxSlo", "agri", "forest", "grass", "urban"), 
  output_attr_name = c("maxSloA", "agriA", "forestA", "grassA", "urbanA"),
  stat = c("max", rep("percent", 4)))

# Plot data with maximum slope per edge as color ramp (steep slopes in red)
dem <- readRAST('dem', ignore.stderr = TRUE)
edges <- readVECT('edges', ignore.stderr = TRUE)
sites <- readVECT('sites', ignore.stderr = TRUE)
lu <- readVECT("landuse", ignore.stderr = TRUE)
plot(dem, col = gray(seq(0,1,length.out=20))) 
col <- adjustcolor(c("red", "green", "blue", "yellow"), alpha.f = 0.3)
plot(lu, add = TRUE, col = col[as.numeric(as.factor(lu$landuse))])
 mm <- range(c(edges$agri_c), na.rm = TRUE) 
b <- seq(from=mm[1],to=mm[2]+diff(mm)*0.01,length.out=10)
c_ramp <- colorRampPalette(c("blue", "red"))
cols <- c_ramp(length(b))[as.numeric(cut(edges$agri_c, breaks = b, right = FALSE))]
plot(edges,col=cols, add = TRUE , lwd=2)
mm <- range(c(sites$agriA), na.rm = TRUE) 
b <- seq(from=mm[1],to=mm[2]+diff(mm)*0.01,length.out=10)
c_ramp <- colorRampPalette(c("blue", "red"))
cols <- c_ramp(length(b))[as.numeric(cut(sites$agriA, breaks = b, right = FALSE))]
plot(sites ,col=cols, add = TRUE, pch = 19)
legend("topleft", col = col, pch = 15, legend = as.factor(sort(unique(lu$landuse))), 
  title = "landuse", ncol = 4)
legend("topright", col = cols[c(1, length(cols))], lwd = 2, 
  legend = paste("precent agri", c(min(sites$agriA), max(sites$agriA))), pch = 19)
# }
# NOT RUN {
                 
# }

Run the code above in your browser using DataLab