In MiKatt/openSTARS: An Open Source Implementation of the 'ArcGIS' Toolbox 'STARS'
openSTARS
knitr::opts_chunk$set(
eval = TRUE,
warning = FALSE,
message = FALSE,
collapse = TRUE,
comment = "#>"
)
Introduction
openSTARS is an open source implementation of the STARS toolbox for ArcGIS (Peterson & Ver Hoef, 2014) using R and GRASS GIS.
It prepares the .ssn object needed for the SSN package.
A digital elevation model (DEM) is used to derive stream networks (in contrast to STARS that can clean an existing stream network). The reason for this is that existing stream networks (e.g. obtained as shape files) very often contain loops and dead ends that hinder building a valid topology for them.
For more information on STARS and SSN, see their web page.
Peterson, E. E., & Ver Hoef, J. M. (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), 1–17.
openSTARS can also be used for other purposes whenever observation sites at streams and their catchments need to be intersected with potential covariates. A detailed description of this package and a step-by-step manual can be found here: \url{https://doi.org/10.1371/journal.pone.0239237}.
Installation and loading
A functional installation of GRASS GIS (>=7.0) with installed add-ons r.stream.basins, r.stream.distance, r.stream.order and
r.hydrodem is needed.
These add-ons can be installed from within GRASS using the console and g.extension or in the GUI under 'Settings'/'Addons extensions'/'Install extensions from add-ons' under 'raster'.
Installation from CRAN repository:
install.packages("openSTARS")
library("openSTARS")
For the lastest development version of openSTARS on GitHub (carefull, might be experimental):
# install.packages("devtools")
devtools::install_github("MiKatt/openSTARS", ref = "dev")
library("openSTARS")
Step by step usage
Initiate and setup GRASS
First, a GRASS session must be initiated and setup. Adjust the paths to the GRASS installation (gisBase) and, if needed, to the GRASS data base where all GRASS files will be stored (gisDbase) to those on your system. The name of the GRASS location within can be given (location). The projection and extent of the GRASS location is based on that one of the digital elevation model (dem) and is used for all input and output files. On Windows systems, you might get a warinig 'WARNING: Concurrent mapset locking is not supported on Windows' that can be ignored.
Please use data in a metric coordinate reference system (CRS) appropriate for the study region (i.e. no long/lat CRS). Otherwise particularly the network correction on 'correct_colplex_confluences' might not work proberly. Suggestions for CRS can e.g. be found here http://epsg.io/.
# give paths to GRASS and where to store the GRASS data base
# Linux e.g.
grass_program_path <- "/usr/lib/grass78/"
# Windows e.g.
# grass_program_path <- "c:/Program Files/GRASS GIS 7.6"
working_dir <- file.path(tempdir(), "grass_workflow")
grass_db_path <- file.path(working_dir, "grassDB")
dir.create(working_dir)
setwd(tempdir())
# specify the path to the digital elevation model
dem_path <- system.file("extdata", "nc", "elev_ned_30m.tif", package = "openSTARS")
setup_grass_environment(dem = dem_path,
gisBase = grass_program_path,
gisDbase = grass_db_path,
location = "nc_openSTARS",
remove_GISRC = TRUE,
override = TRUE
)
gmeta()
Load data into GRASS
Use import_data to import all data into GRASS (DEM, observation sites and other optional data). Optional data includes a stream network to burn into the DEM (see derive_streams), prediction sites if they have been already created with a different program (alternatively, prediction sites can be created using calc_prediction_sites), and raster and vector maps of potential predictor variables for the SSN model that can later be intersected with the catchments of the sites (calc_attributes_edges and calc_attributes_sites_approx, or calc_attributes_sites_exact).
Note that vector data is reprojected on the fly in this import, while raster data is not as it should be done manually (due to the changes this might cause to the resolution etc.). It can be checked before the import if the current region (as defined by the dem) and other raster maps have the same projection. If this is not the case, they should be reprojected before, so that the dem and all raster files have identical projections. Please note that providing all data including vector data in the same coordinate reference system (CRS) is best to avoid errors during on the fly import.
preds_r_path <- system.file("extdata", "nc", "landuse_r.tif", package = "openSTARS")
check_projection(preds_r_path)
# the observation sites
sites_path <- system.file("extdata", "nc", "sites_nc.shp", package = "openSTARS")
# potential predictor in raster format
preds_r_path <- system.file("extdata", "nc", "landuse_r.tif", package = "openSTARS")
# potential predictors in vector format
preds_v_path <- c(system.file("extdata", "nc", "geology.shp", package = "openSTARS"),
system.file("extdata", "nc", "pointsources.shp", package = "openSTARS"))
# existing stream network
streams_path <- system.file("extdata", "nc", "streams.shp", package = "openSTARS")
import_data(dem = dem_path, sites = sites_path, streams = streams_path,
predictor_vector = preds_v_path, predictor_v_names = c("geology", "psources"),
predictor_raster = preds_r_path)
The DEM is loaded into the GRASS database as a raster map named dem, the sites as a vector map named sites_o and the (optional) stream network as a vector map named streams_o. Predictor sites are stored under their base file name, potential predictors either using their base file names or the ones provided in 'predictor_r_names' and 'predictor_v_names', respectively.
The data looks like this:
library(sp)
dem <- readRAST("dem", ignore.stderr = TRUE)
sites <- readVECT("sites_o", ignore.stderr = TRUE)
psources <- readVECT("psources", ignore.stderr = TRUE)
lu <- readRAST("landuse_r", ignore.stderr = TRUE)
plot(dem, col = gray(seq(0,1,length.out=20)))
cols <- colorRampPalette(c("blue", "red"))(length(sites$value))[rank(sites$value)]
points(sites, pch = 16, col = cols)
points(psources, pch = 19, col = 1, cex = 1.7)
legend(x = par("usr")[2]*0.991, y = par("usr")[4]*0.999, pch = c(16, 16, 19), ncol = 1, col = c(range(cols),1), legend = c(paste("value at sites:", c(range(sites$value))), "point sources"))
col <- adjustcolor(c("orange", "yellow", "green", "greenyellow", "green3","blue", "brown"), alpha.f = 1)
image(lu["landuse_r"], col = col)
legend("right", col = col, pch = 15, legend = as.factor(sort(unique(lu$landuse_r))),
title = "landuse")
Derive streams from DEM
Next, the streams must be derived from the DEM.
derive_streams()
An existing stream network (if provided to import_data before) can be burnt into the DEM to force the streams derived from the DEM to mapped ones. It is not possible to use a given stream network directly but it has to be derived from the DEM because otherwise it lacks topological information needed in the consecutive steps. Additional specifications on how the streams shall be created can be provided (see ?derive_streams and the GRASS function r.stream.extract for details).
derive_streams(accum_threshold = 100, condition = T, clean = TRUE, burn = 10)
dem <- readRAST("dem", ignore.stderr = TRUE)
streams <- readVECT("streams_v", ignore.stderr = TRUE)
plot(dem, col = terrain.colors(20))
lines(streams, col = "blue")
cols <- colorRampPalette(c("blue", "red"))(length(sites$value))[rank(sites$value)]
points(sites, pch = 16, col = cols)
Check and correct the network
Next, the stream network should be checked if there are stream segments with more than two inflows. These must be corrected because the .ssn object must not have such complex confluences. In the nc data set provided, there will be complex confluences only if accum_threshold is small (e.g. 150; note that this will take long to run and use a lot of memory).
cp <- check_compl_confluences()
if (cp)
correct_compl_confluences()
dem <- readRAST('dem', ignore.stderr = TRUE)
streams <- readVECT('streams_v', ignore.stderr = TRUE)
streams_orig <- readVECT('streams_v_o3', ignore.stderr = TRUE)
# zoom to a relevant part of the dem
plot(dem, col = terrain.colors(20), axes = TRUE,
xlim = c(640080,640170), ylim = c(219730,219820))
lines(streams_orig, col = 'red', lwd = 4)
lines(streams, col = 'blue', lty = 2, lwd = 2)
legend("topright", col = c("red", "blue"), lty = c(1,2), lwd = c(4,2),
legend = c("original", "corrected"))
The end node of one segment is moved a small fraction of the cell size upstream on a neighbouring segment. That segement is split into two parts at the new confluence. All features are corrected accordingly (cat, stream, length, prev_str01, prev_str02, next_str etc.).
Other topological errors as mentioned for the ArcGIS toolbox STARS do not occur
if the stream network is derived from a DEM.
Prepare edges
Now, information needed for the .ssn object can be derived for the streams and stored in a new vector map edges.
calc_edges()
edges <- readVECT("edges", ignore.stderr = TRUE)
head(edges@data, n = 4)
edges now holds the derived network plus attributes needed for the .ssn object
- network identifier (netID)
- reach identifier (rid)
- stream segment length (length)
- distance from the source (sourceDist)
- upstream distance, i.e. distance from the outlet of the network to the start (upstream node) of the stream segment (upDist)
- total catchment area (H2OArea)
- reach contributing area (rcaArea)
The additional fields hold information about the network: 'next_str' is the 'stream' this segment flows into, 'prev_str01' and 'prev_str02' are the two segments that flow into this segment.
Prepare sites
Often, survey sites do not lay exactly on the stream network (due to GPS imprecision, stream representation as lines, derivation of streams from a DEM, etc.). To assign an exact position of the sites on the network they are moved to the closest stream segment (snapped) using the GRASS function
v.distance. Additionally, attributes needed for .ssn object are assigned:
calc_sites()
sites <- readVECT("sites", ignore.stderr = TRUE)
head(sites@data, n = 4)
sites now holds the corrected sites plus attributes needed for the .ssn object
- point identifier (pid)
- location identifier (locID)
- network identifier (netID)
- reach identifier of the edge segment the point lies on (rid)
- upstream distance (upDist), i.e. the distance to the network outlet calculated using r.stream.distance.
- distance ratio, i.e. the ratio of the distance from the outflow of the edge to the point along the edge and the total length of the edge segment (distRatio).
Additional fields hold information on the snapping: distance of the original site to the closest edge ('dist'), i.e. how far the point was moved, and the new x and y coordinates ('xm', 'ym'). The filed 'cat_edge' gives the 'cat' of the stream segment the point lies on. It is used to identify the edge the point lies on to extract the 'rid'.
It is possible to give a maximum distance a point is allowed to be moved ('maxdist'). Points farther away from any stream are deleted.
calc_sites(maxdist = 100)
dem <- readRAST("dem", ignore.stderr = TRUE)
sites <- readVECT("sites", ignore.stderr = TRUE)
sites_orig <- readVECT("sites_o", ignore.stderr = TRUE)
edges <- readVECT("edges", ignore.stderr = TRUE)
plot(dem, col = terrain.colors(20))
lines(edges, col = "blue")
points(sites_orig, pch = 20, col = "black")
points(sites, pch = 21, cex=0.75, bg = "grey")
legend("topright", y = par("usr")[3]*1.01, col = 1, pt.bg = "grey", pch = c(21, 19), legend = c("snapped sites", "original sites"), ncol = 2)
Prepare prediction sites
Prediction sites can be created along the streams. Either the distance between the sites must be provided (dist) or the approximate number of sites that shall be created (nsites). Additionally, the creation can be restricted to certain networks (netIDs). The sites will be assigned regularly on the stream network. If prediction sites with specifec coordinates are needed, they should be created manually.
Similar as for the observation sites, attributes needed for .ssn object are assigned:
- point identifier (pid)
- location identifier (locID)
- network identifier (netID)
- reach identifier of the edge segment the point lies on (rid)
- upstream distance (upDist), i.e. the distance to the network outlet calculated using r.stream.distance.
- distance ratio, i.e. the ratio of the distance from the outflow of the edge to the point along the edge and the total length of the edge segment (distRatio).
The filed 'cat_edge' gives the 'cat' of the stream segment the point lies on (equivalent to 'rid').
calc_prediction_sites(predictions = "preds", nsites = 100, netIDs = 50 )
dem <- readRAST("dem", ignore.stderr = TRUE)
sites <- readVECT("sites", ignore.stderr = TRUE)
pred_sites <- readVECT("preds", ignore.stderr = TRUE)
edges <- readVECT("edges", ignore.stderr = TRUE)
plot(dem, col = terrain.colors(20))
lines(edges, col = "blue")
points(sites, pch = 21, cex=0.75, bg = "grey")
points(pred_sites, pch = 21, cex=0.75, bg = "royalblue")
legend("topright", ncol = 2, pt.bg = c("grey","royalblue"), pch = 21, legend = c("(snapped) observation sites","prediction sites"))
head(pred_sites@data, n = 5)
If prediction sites were prepared before in a different software and were loaded with 'import_data' they should be treated in the same way as the observation sites:
calc_sites(pred_site = name_of_prediction_sites)
Calculate attributes from raster and vector maps
Attributes (i.e. predictor variables for the .ssn object) can be calculated for observation and prediction sites. There are two ways to calculates attributes:
- approximately 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).
- exactly by intersecting the catchment of each point with raster maps;
For the approximate calculation, first attributes must be intersected with the sub-catchments of the stream segments and then they are assigned to each site based on the distance ratio of the point. Note that the sub-catchment area 'H2OArea' for each stream segment is calculated automatically in calc_edges.
# calculate slope from DEM as an example attribute
execGRASS("r.slope.aspect", flags = c("overwrite","quiet"),
parameters = list(
elevation = "dem",
slope = "slope"
))
# calculate average slope per sub-catchment of each stream segment using raster data
calc_attributes_edges(input_raster = c("slope", "landuse_r"), stat_rast = c("mean", "percent"),
attr_name_rast = c("avSlo","lu"),
round_dig = 4)
#calculate approx. catchment area, average slope per catchment and percentage of land use class of each site
calc_attributes_sites_approx(sites_map = "sites",
input_attr_name = c("avSlo",paste0("lup_",1:7)),
output_attr_name = c("avSloA",paste0("luA", 1:7)),
stat = c("mean", rep("percent", 7)))
sites <- readVECT("sites", ignore.stderr = TRUE)
head(sites@data, n = 5)
The exact calculation of attribute values for the total catchment of each point can take quite long (depending on the number of points): For each point, first the total catchment is delineated based on the DEM and then intersected with the map(s) provided. It must be decided on a case by case basis if the approximate calculation is good enough.
# calculate exact catchment area, average slope per catchment and number of point sources of each site
calc_attributes_sites_exact(sites_map = "sites",
input_raster = "slope",
stat_rast = "mean",
attr_name_rast = "avSloE",
input_vector = "psources",
stat_vect = "count",
attr_name_vect = "nps",
round_dig = 4)
sites <- readVECT("sites", ignore.stderr = TRUE)
head(sites@data, n = 5)
In both alternatives, the catchment area for each site is calculated automatically ('H2OAreaA' for calc_attributes_sites_appox and 'H2OArea' for calc_attributes_sites_exact).
If predictions are to be made in a later step, the same attributes must be calculated for the prediction sites.
Write all files to an ssn folder
All files needed (edges, sites and optionally prediction sites) are written to the file path provided and can then be read in by the SSN package.
ssn_dir <- file.path(tempdir(), 'nc.ssn')
export_ssn(ssn_dir)
list.files(ssn_dir)
Try with SSN package
library(SSN)
# import
ssn_obj <- importSSN(ssn_dir, o.write = TRUE)
plot(ssn_obj, 'value')
# Create Distance Matrix
createDistMat(ssn_obj, o.write = TRUE)
dmats <- getStreamDistMat(ssn_obj)
ssn_obj.Torg <- Torgegram(ssn_obj, "value", nlag = 20, maxlag = 15000)
plot(ssn_obj.Torg)
names(ssn_obj@data)
names(ssn_obj)
ssn_obj <- additive.function(ssn_obj, "H2OArea", "computed.afv")
# non-spatial model
ssn_obj.glmssn0 <- glmssn(value ~ upDist, ssn.object = ssn_obj,
CorModels = NULL)
summary(ssn_obj.glmssn0)
# same as
summary(lm(value ~ upDist, getSSNdata.frame(ssn_obj)))
# # # spatial model
ssn_obj.glmssn1 <- glmssn(value ~ upDist , ssn.object = ssn_obj,
CorModels = c("Exponential.taildown", "Exponential.tailup"),
addfunccol = "computed.afv")
summary(ssn_obj.glmssn1)
varcomp(ssn_obj.glmssn1)
Contributors
Want to contribute?
Checkout our contribution guide here.
Meta
- Please report any issues, bugs or feature requests.
- License: MIT
MiKatt/openSTARS documentation built on June 17, 2022, 5:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
openSTARS
knitr::opts_chunk$set( eval = TRUE, warning = FALSE, message = FALSE, collapse = TRUE, comment = "#>" )
Introduction
openSTARS is an open source implementation of the STARS toolbox for ArcGIS (Peterson & Ver Hoef, 2014) using R and GRASS GIS.
It prepares the .ssn object needed for the SSN package.
A digital elevation model (DEM) is used to derive stream networks (in contrast to STARS that can clean an existing stream network). The reason for this is that existing stream networks (e.g. obtained as shape files) very often contain loops and dead ends that hinder building a valid topology for them.
For more information on STARS and SSN, see their web page.
Peterson, E. E., & Ver Hoef, J. M. (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), 1–17.
openSTARS can also be used for other purposes whenever observation sites at streams and their catchments need to be intersected with potential covariates. A detailed description of this package and a step-by-step manual can be found here: \url{https://doi.org/10.1371/journal.pone.0239237}.
Installation and loading
A functional installation of GRASS GIS (>=7.0) with installed add-ons r.stream.basins, r.stream.distance, r.stream.order and r.hydrodem is needed. These add-ons can be installed from within GRASS using the console and g.extension or in the GUI under 'Settings'/'Addons extensions'/'Install extensions from add-ons' under 'raster'.
Installation from CRAN repository:
install.packages("openSTARS") library("openSTARS")
For the lastest development version of openSTARS on GitHub (carefull, might be experimental):
# install.packages("devtools") devtools::install_github("MiKatt/openSTARS", ref = "dev") library("openSTARS")
Step by step usage
Initiate and setup GRASS
First, a GRASS session must be initiated and setup. Adjust the paths to the GRASS installation (gisBase) and, if needed, to the GRASS data base where all GRASS files will be stored (gisDbase) to those on your system. The name of the GRASS location within can be given (location). The projection and extent of the GRASS location is based on that one of the digital elevation model (dem) and is used for all input and output files. On Windows systems, you might get a warinig 'WARNING: Concurrent mapset locking is not supported on Windows' that can be ignored.
Please use data in a metric coordinate reference system (CRS) appropriate for the study region (i.e. no long/lat CRS). Otherwise particularly the network correction on 'correct_colplex_confluences' might not work proberly. Suggestions for CRS can e.g. be found here http://epsg.io/.
# give paths to GRASS and where to store the GRASS data base # Linux e.g. grass_program_path <- "/usr/lib/grass78/" # Windows e.g. # grass_program_path <- "c:/Program Files/GRASS GIS 7.6" working_dir <- file.path(tempdir(), "grass_workflow") grass_db_path <- file.path(working_dir, "grassDB") dir.create(working_dir) setwd(tempdir()) # specify the path to the digital elevation model dem_path <- system.file("extdata", "nc", "elev_ned_30m.tif", package = "openSTARS") setup_grass_environment(dem = dem_path, gisBase = grass_program_path, gisDbase = grass_db_path, location = "nc_openSTARS", remove_GISRC = TRUE, override = TRUE )
gmeta()
Load data into GRASS
Use import_data to import all data into GRASS (DEM, observation sites and other optional data). Optional data includes a stream network to burn into the DEM (see derive_streams), prediction sites if they have been already created with a different program (alternatively, prediction sites can be created using calc_prediction_sites), and raster and vector maps of potential predictor variables for the SSN model that can later be intersected with the catchments of the sites (calc_attributes_edges and calc_attributes_sites_approx, or calc_attributes_sites_exact).
Note that vector data is reprojected on the fly in this import, while raster data is not as it should be done manually (due to the changes this might cause to the resolution etc.). It can be checked before the import if the current region (as defined by the dem) and other raster maps have the same projection. If this is not the case, they should be reprojected before, so that the dem and all raster files have identical projections. Please note that providing all data including vector data in the same coordinate reference system (CRS) is best to avoid errors during on the fly import.
preds_r_path <- system.file("extdata", "nc", "landuse_r.tif", package = "openSTARS") check_projection(preds_r_path)
# the observation sites sites_path <- system.file("extdata", "nc", "sites_nc.shp", package = "openSTARS") # potential predictor in raster format preds_r_path <- system.file("extdata", "nc", "landuse_r.tif", package = "openSTARS") # potential predictors in vector format preds_v_path <- c(system.file("extdata", "nc", "geology.shp", package = "openSTARS"), system.file("extdata", "nc", "pointsources.shp", package = "openSTARS")) # existing stream network streams_path <- system.file("extdata", "nc", "streams.shp", package = "openSTARS") import_data(dem = dem_path, sites = sites_path, streams = streams_path, predictor_vector = preds_v_path, predictor_v_names = c("geology", "psources"), predictor_raster = preds_r_path)
The DEM is loaded into the GRASS database as a raster map named dem, the sites as a vector map named sites_o and the (optional) stream network as a vector map named streams_o. Predictor sites are stored under their base file name, potential predictors either using their base file names or the ones provided in 'predictor_r_names' and 'predictor_v_names', respectively.
The data looks like this:
library(sp) dem <- readRAST("dem", ignore.stderr = TRUE) sites <- readVECT("sites_o", ignore.stderr = TRUE) psources <- readVECT("psources", ignore.stderr = TRUE) lu <- readRAST("landuse_r", ignore.stderr = TRUE) plot(dem, col = gray(seq(0,1,length.out=20))) cols <- colorRampPalette(c("blue", "red"))(length(sites$value))[rank(sites$value)] points(sites, pch = 16, col = cols) points(psources, pch = 19, col = 1, cex = 1.7) legend(x = par("usr")[2]*0.991, y = par("usr")[4]*0.999, pch = c(16, 16, 19), ncol = 1, col = c(range(cols),1), legend = c(paste("value at sites:", c(range(sites$value))), "point sources"))
col <- adjustcolor(c("orange", "yellow", "green", "greenyellow", "green3","blue", "brown"), alpha.f = 1) image(lu["landuse_r"], col = col) legend("right", col = col, pch = 15, legend = as.factor(sort(unique(lu$landuse_r))), title = "landuse")
Derive streams from DEM
Next, the streams must be derived from the DEM.
derive_streams()
An existing stream network (if provided to import_data before) can be burnt into the DEM to force the streams derived from the DEM to mapped ones. It is not possible to use a given stream network directly but it has to be derived from the DEM because otherwise it lacks topological information needed in the consecutive steps. Additional specifications on how the streams shall be created can be provided (see ?derive_streams and the GRASS function r.stream.extract for details).
derive_streams(accum_threshold = 100, condition = T, clean = TRUE, burn = 10)
dem <- readRAST("dem", ignore.stderr = TRUE) streams <- readVECT("streams_v", ignore.stderr = TRUE) plot(dem, col = terrain.colors(20)) lines(streams, col = "blue") cols <- colorRampPalette(c("blue", "red"))(length(sites$value))[rank(sites$value)] points(sites, pch = 16, col = cols)
Check and correct the network
Next, the stream network should be checked if there are stream segments with more than two inflows. These must be corrected because the .ssn object must not have such complex confluences. In the nc data set provided, there will be complex confluences only if accum_threshold is small (e.g. 150; note that this will take long to run and use a lot of memory).
cp <- check_compl_confluences() if (cp) correct_compl_confluences()
dem <- readRAST('dem', ignore.stderr = TRUE) streams <- readVECT('streams_v', ignore.stderr = TRUE) streams_orig <- readVECT('streams_v_o3', ignore.stderr = TRUE) # zoom to a relevant part of the dem plot(dem, col = terrain.colors(20), axes = TRUE, xlim = c(640080,640170), ylim = c(219730,219820)) lines(streams_orig, col = 'red', lwd = 4) lines(streams, col = 'blue', lty = 2, lwd = 2) legend("topright", col = c("red", "blue"), lty = c(1,2), lwd = c(4,2), legend = c("original", "corrected"))
The end node of one segment is moved a small fraction of the cell size upstream on a neighbouring segment. That segement is split into two parts at the new confluence. All features are corrected accordingly (cat, stream, length, prev_str01, prev_str02, next_str etc.).
Other topological errors as mentioned for the ArcGIS toolbox STARS do not occur if the stream network is derived from a DEM.
Prepare edges
Now, information needed for the .ssn object can be derived for the streams and stored in a new vector map edges.
calc_edges()
edges <- readVECT("edges", ignore.stderr = TRUE) head(edges@data, n = 4)
edges now holds the derived network plus attributes needed for the .ssn object
- network identifier (netID)
- reach identifier (rid)
- stream segment length (length)
- distance from the source (sourceDist)
- upstream distance, i.e. distance from the outlet of the network to the start (upstream node) of the stream segment (upDist)
- total catchment area (H2OArea)
- reach contributing area (rcaArea)
The additional fields hold information about the network: 'next_str' is the 'stream' this segment flows into, 'prev_str01' and 'prev_str02' are the two segments that flow into this segment.
Prepare sites
Often, survey sites do not lay exactly on the stream network (due to GPS imprecision, stream representation as lines, derivation of streams from a DEM, etc.). To assign an exact position of the sites on the network they are moved to the closest stream segment (snapped) using the GRASS function v.distance. Additionally, attributes needed for .ssn object are assigned:
calc_sites()
sites <- readVECT("sites", ignore.stderr = TRUE) head(sites@data, n = 4)
sites now holds the corrected sites plus attributes needed for the .ssn object
- point identifier (pid)
- location identifier (locID)
- network identifier (netID)
- reach identifier of the edge segment the point lies on (rid)
- upstream distance (upDist), i.e. the distance to the network outlet calculated using r.stream.distance.
- distance ratio, i.e. the ratio of the distance from the outflow of the edge to the point along the edge and the total length of the edge segment (distRatio).
Additional fields hold information on the snapping: distance of the original site to the closest edge ('dist'), i.e. how far the point was moved, and the new x and y coordinates ('xm', 'ym'). The filed 'cat_edge' gives the 'cat' of the stream segment the point lies on. It is used to identify the edge the point lies on to extract the 'rid'.
It is possible to give a maximum distance a point is allowed to be moved ('maxdist'). Points farther away from any stream are deleted.
calc_sites(maxdist = 100)
dem <- readRAST("dem", ignore.stderr = TRUE) sites <- readVECT("sites", ignore.stderr = TRUE) sites_orig <- readVECT("sites_o", ignore.stderr = TRUE) edges <- readVECT("edges", ignore.stderr = TRUE) plot(dem, col = terrain.colors(20)) lines(edges, col = "blue") points(sites_orig, pch = 20, col = "black") points(sites, pch = 21, cex=0.75, bg = "grey") legend("topright", y = par("usr")[3]*1.01, col = 1, pt.bg = "grey", pch = c(21, 19), legend = c("snapped sites", "original sites"), ncol = 2)
Prepare prediction sites
Prediction sites can be created along the streams. Either the distance between the sites must be provided (dist) or the approximate number of sites that shall be created (nsites). Additionally, the creation can be restricted to certain networks (netIDs). The sites will be assigned regularly on the stream network. If prediction sites with specifec coordinates are needed, they should be created manually.
Similar as for the observation sites, attributes needed for .ssn object are assigned:
- point identifier (pid)
- location identifier (locID)
- network identifier (netID)
- reach identifier of the edge segment the point lies on (rid)
- upstream distance (upDist), i.e. the distance to the network outlet calculated using r.stream.distance.
- distance ratio, i.e. the ratio of the distance from the outflow of the edge to the point along the edge and the total length of the edge segment (distRatio).
The filed 'cat_edge' gives the 'cat' of the stream segment the point lies on (equivalent to 'rid').
calc_prediction_sites(predictions = "preds", nsites = 100, netIDs = 50 )
dem <- readRAST("dem", ignore.stderr = TRUE) sites <- readVECT("sites", ignore.stderr = TRUE) pred_sites <- readVECT("preds", ignore.stderr = TRUE) edges <- readVECT("edges", ignore.stderr = TRUE) plot(dem, col = terrain.colors(20)) lines(edges, col = "blue") points(sites, pch = 21, cex=0.75, bg = "grey") points(pred_sites, pch = 21, cex=0.75, bg = "royalblue") legend("topright", ncol = 2, pt.bg = c("grey","royalblue"), pch = 21, legend = c("(snapped) observation sites","prediction sites")) head(pred_sites@data, n = 5)
If prediction sites were prepared before in a different software and were loaded with 'import_data' they should be treated in the same way as the observation sites:
calc_sites(pred_site = name_of_prediction_sites)
Calculate attributes from raster and vector maps
Attributes (i.e. predictor variables for the .ssn object) can be calculated for observation and prediction sites. There are two ways to calculates attributes:
- approximately 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).
- exactly by intersecting the catchment of each point with raster maps;
For the approximate calculation, first attributes must be intersected with the sub-catchments of the stream segments and then they are assigned to each site based on the distance ratio of the point. Note that the sub-catchment area 'H2OArea' for each stream segment is calculated automatically in calc_edges.
# calculate slope from DEM as an example attribute execGRASS("r.slope.aspect", flags = c("overwrite","quiet"), parameters = list( elevation = "dem", slope = "slope" )) # calculate average slope per sub-catchment of each stream segment using raster data calc_attributes_edges(input_raster = c("slope", "landuse_r"), stat_rast = c("mean", "percent"), attr_name_rast = c("avSlo","lu"), round_dig = 4) #calculate approx. catchment area, average slope per catchment and percentage of land use class of each site calc_attributes_sites_approx(sites_map = "sites", input_attr_name = c("avSlo",paste0("lup_",1:7)), output_attr_name = c("avSloA",paste0("luA", 1:7)), stat = c("mean", rep("percent", 7))) sites <- readVECT("sites", ignore.stderr = TRUE) head(sites@data, n = 5)
The exact calculation of attribute values for the total catchment of each point can take quite long (depending on the number of points): For each point, first the total catchment is delineated based on the DEM and then intersected with the map(s) provided. It must be decided on a case by case basis if the approximate calculation is good enough.
# calculate exact catchment area, average slope per catchment and number of point sources of each site calc_attributes_sites_exact(sites_map = "sites", input_raster = "slope", stat_rast = "mean", attr_name_rast = "avSloE", input_vector = "psources", stat_vect = "count", attr_name_vect = "nps", round_dig = 4) sites <- readVECT("sites", ignore.stderr = TRUE) head(sites@data, n = 5)
In both alternatives, the catchment area for each site is calculated automatically ('H2OAreaA' for calc_attributes_sites_appox and 'H2OArea' for calc_attributes_sites_exact).
If predictions are to be made in a later step, the same attributes must be calculated for the prediction sites.
Write all files to an ssn folder
All files needed (edges, sites and optionally prediction sites) are written to the file path provided and can then be read in by the SSN package.
ssn_dir <- file.path(tempdir(), 'nc.ssn') export_ssn(ssn_dir) list.files(ssn_dir)
Try with SSN package
library(SSN) # import ssn_obj <- importSSN(ssn_dir, o.write = TRUE) plot(ssn_obj, 'value') # Create Distance Matrix createDistMat(ssn_obj, o.write = TRUE) dmats <- getStreamDistMat(ssn_obj) ssn_obj.Torg <- Torgegram(ssn_obj, "value", nlag = 20, maxlag = 15000) plot(ssn_obj.Torg) names(ssn_obj@data) names(ssn_obj) ssn_obj <- additive.function(ssn_obj, "H2OArea", "computed.afv") # non-spatial model ssn_obj.glmssn0 <- glmssn(value ~ upDist, ssn.object = ssn_obj, CorModels = NULL) summary(ssn_obj.glmssn0) # same as summary(lm(value ~ upDist, getSSNdata.frame(ssn_obj))) # # # spatial model ssn_obj.glmssn1 <- glmssn(value ~ upDist , ssn.object = ssn_obj, CorModels = c("Exponential.taildown", "Exponential.tailup"), addfunccol = "computed.afv") summary(ssn_obj.glmssn1) varcomp(ssn_obj.glmssn1)
Contributors
Want to contribute?
Checkout our contribution guide here.
Meta
- Please report any issues, bugs or feature requests.
- License: MIT
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.