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

3
  • 1
    Hi! I don't know the details but I think you may want to check slopes package. cc @robinlovelace – agila Jan 20 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 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 at 21:24
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)
1
  • @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 at 22:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.