README.Rmd

---
output:
  github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

openSTARS
=============

```{r echo=FALSE}
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](http://www.fs.fed.us/rm/boise/AWAE/projects/SpatialStreamNetworks.shtml).

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)](https://grass.osgeo.org/#) with installed add-ons [r.stream.basins](https://grass.osgeo.org/grass78/manuals/addons/r.stream.basins.html), [r.stream.distance](https://grass.osgeo.org/grass78/manuals/addons/r.stream.distance.html), [r.stream.order](https://grass.osgeo.org/grass78/manuals/addons/r.stream.order.html) and 
[r.hydrodem](https://grass.osgeo.org/grass78/manuals/addons/r.hydrodem.html) 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:
```{r install, eval=FALSE}
install.packages("openSTARS")
library("openSTARS")
```

For the lastest development version of openSTARS on GitHub (carefull, might be experimental):
```{r install_gh, eval = TRUE}
# 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/.

```{r setup_grass}
# 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
                        )
```

```{r gmeta}
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.

```{r check_prij, message = TRUE}
preds_r_path <- system.file("extdata", "nc", "landuse_r.tif", package = "openSTARS")
check_projection(preds_r_path)
```

```{r import_data, message = TRUE}
# 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:

```{r plot_data1.1, message=FALSE, results='hide'}
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"))

```
```{r plot_data1.2, message=FALSE, results='hide'}
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.
```{r derive_streams, results='hide', message=TRUE, eval = FALSE}
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](https://grass.osgeo.org/grass78/manuals/r.stream.extract.html) for details).

```{r derive_streams2, results='hide', message=TRUE, eval=TRUE}
derive_streams(accum_threshold = 100, condition = T, clean = TRUE, burn = 10)
```

```{r plot_data2, message=FALSE, warning=FALSE, results="hide"}
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).

```{r compl_confl, eval = TRUE}
cp <- check_compl_confluences()

if (cp)
  correct_compl_confluences()
```

```{r plot_compl_confl, eval = TRUE}
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"))
```


<!-- An example of a complex junction and the correction would look like this: -->

<!-- ![Original network with complex junction. Arrows indicate flow direction.](README_files/compl_junction1.PNG)  ![Corrected network. Arrows indicate flow direction.](README_files/compl_junction2.PNG) -->

<!-- *Left* Original network with complex junction (i.e. three inflows to one outflow).  -->
<!-- *Right* Corrected network. Arrows indicate flow direction. -->


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

```{r prep_edges, results='hide', message=FALSE, warning=FALSE}
calc_edges()
```

```{r 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](https://grass.osgeo.org/grass78/manuals/v.distance.html). Additionally, attributes needed for .ssn object are assigned: 

```{r prep_sites, results='hide', message=FALSE, warning=FALSE}
calc_sites()
```

```{r show_sites, message=FALSE, warning=FALSE}
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](https://grass.osgeo.org/grass78/manuals/addons/r.stream.distance.html).
* 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.

```{r prep_sites_max, results='hide', message=FALSE, warning=FALSE, eval = FALSE}
calc_sites(maxdist = 100)
```

```{r plot_data3}
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](https://grass.osgeo.org/grass78/manuals/addons/r.stream.distance.html).
* 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').

```{r prep_pred_sites, results='hide', message=FALSE, warning=FALSE}
calc_prediction_sites(predictions = "preds", nsites = 100, netIDs = 50 )
```

```{r plot_data4}
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:

```{r prep_predsites, results='hide', message=FALSE, warning=FALSE, eval = FALSE}
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: 

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

```{r attributes_approx}
# 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.


```{r attributes_exact}
# 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.

```{r export}
ssn_dir <- file.path(tempdir(), 'nc.ssn')
export_ssn(ssn_dir)
list.files(ssn_dir)
```


#### Try with SSN package
```{r ssn_test}
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

+ [Mira Kattwinkel](https://github.com/MiKatt)
+ [Eduard Szöcs](https://github.com/EDiLD)

### Want to contribute?

Checkout our [contribution guide here](https://github.com/edild/openSTARS/blob/master/CONTRIBUTING.md).

### Meta

* Please [report any issues, bugs or feature requests](https://github.com/MiKatt/openSTARS/issues).
* License: MIT