0

I would like to calculate the slope for stream segments in R. The stream network has been derived from the DEM. I have split the stream segments into equal parts of 100m (or less) and I have the DEM (1m resolution). I have attempted to use the openSTARS package that runs through GRASS without success (I am finding the initial steps to start using GRASS very difficult).

I am aware of tools in QGIS and ArcGIS but I would like to make as many of my steps within R as possible (Calculate slope of line segments with QGIS ; Workflow for determining stream gradient?)

user3386170
  • 1,957
  • 1
  • 28
  • 53
  • 1
    Hi! I don't know the details but I think you may want to check slopes package. cc @robinlovelace – agila Jan 20 '21 at 20:56
  • If you sample your stream line vertex points over the DEM raster (which I assume you have?) then you have your stream points in 3d and you can compute slope by taking pairs of vertices and then slope is something like sqrt(dx^2+dy^2)/dz, and aspect by atan2(dx, dy) where the ds are differences in coordinates. – Spacedman Jan 20 '21 at 21:06
  • @agila slopes was so straightforward! Hours of searching culminating in a two-minute function. @Spacedman Could you provide more detail on determining the aspect? – user3386170 Jan 20 '21 at 21:24

1 Answers1

0

I managed to implement the approach suggested by @Spacedman to calculate slope. I am leaving the code here as a possible starter to calculate aspect. (the data called up here from the sample is a dummy and doesn't fit the names used)

I borrowed extensively from: https://www.r-spatial.org/r/2019/09/26/spatial-networks.html

    if(!"remotes" %in% installed.packages()) {
  install.packages("remotes")
}

cran_pkgs = c( #"sf", "tidygraph", "igraph", "osmdata", #"dplyr", #"tibble", #"ggplot2", "units",

"tmap",

#"rgrass7", "link2GI", "nabor" )

remotes::install_cran(cran_pkgs)

library(sf) library(tidygraph) library(igraph) library(dplyr) library(tibble) library(ggplot2) library(units) library(tmap) library(osmdata) library(rgrass7) library(link2GI) library(nabor) #for sample dataset taken from slopes package if(!require(slopes)){install.packages("slopes")} library(slopes) data(dem_lisbon_raster) data(lisbon_road_segments) dem<-dem_lisbon_raster stream_split<-lisbon_road_segments%>% select(OBJECTID) %>% rename(edgeID=OBJECTID)

#assign identifier to each stream segment stream_split <- stream_split %>% mutate(edgeID = c(1:n()))

#extract start and end nodes nodes<-stream_split %>%st_coordinates() %>% as_tibble %>% rename(edgeID=L1) %>% group_by(edgeID)%>%slice(c(1, n())) %>% ungroup() %>% mutate(start_end = rep(c('start', 'end'), times = n()/2))

#assign a unique index to each node (duplicate coordinates get the same index) nodes <- nodes %>% mutate(xy = paste(.$X, .$Y)) %>% mutate(nodeID = group_indices(., factor(xy, levels = unique(xy)))) %>% select(-xy)

#convert nodes to sf object nodes_sf <- nodes %>% distinct(nodeID, .keep_all = TRUE) %>% select(-c(edgeID, start_end)) %>% st_as_sf(coords = c('X', 'Y')) %>% st_set_crs(st_crs(stream_split))

#extract elevation at each node nodes_sf$elev<-raster::extract(dem, nodes_sf)

##add notes to each edge

source_nodes <- nodes %>% filter(start_end == 'start') %>% pull(nodeID)

target_nodes <- nodes %>% filter(start_end == 'end') %>% pull(nodeID)

stream_split = stream_split %>% mutate(from = source_nodes, to = target_nodes)

#convert elevations to tibble nodes.tib<-as_tibble(nodes_sf)%>%select(-geometry)

#add start and end elevations to segments stream_split<-inner_join(stream_split, nodes.tib, by=c("from"="nodeID")) %>% rename(elev_start=elev) stream_split<-inner_join(stream_split, nodes.tib, by=c("to"="nodeID")) %>% rename(elev_end=elev) stream_split$length<-st_length(stream_split)#caculate segment length stream_split$slope<-(stream_split$elev_end-stream_split$elev_start)/stream_split$length

However, I do not get the same result as when I use the slopes package.

stream_split$slope_package = slope_raster(stream_split, e = dem)
user3386170
  • 1,957
  • 1
  • 28
  • 53
  • @Spacedman This is my attempt to implement your suggestion. Since it is too long to include as a comment, I've posted as an answer but I'll pull it down if you post an answer. Perhaps, this code stub will be helpful. – user3386170 Jan 20 '21 at 22:31