Calculate aspect for stream segments in R from 2D linestring and DEM

Question

I would like to calculate the general aspect for stream segments in R. The stream network has been derived from the DEM and is in 2D. I have split the stream segments into equal parts of 100m (or less) and I have the DEM (1m resolution). While the aspect may change (slightly) along each segment, I am interested in the overall or dominant aspect for that segment based on the hill slope from start to finish. [I can see it might possible to use the DEM to create a 3D linestring but I don't know how]

When I use a basic tool to calculate aspect for each cell/pixel of a DEM raster, it shows opposite facing-slopes on either side of the stream (logically). Here, I want the aspect of the stream, i.e. of the general slope/aspect of the hill along which the stream flows. If the stream is flowing along a hill that is north-facing, the aspect that I am hoping to obtain is North (or 0). I don't want the direction of flow of the stream (from x1,x2 to y1,y2) although it may overlap given that the stream will probably flow in the most direct route downhill.

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

I had initially asked a question about slope and aspect but only got an answer for slope so I'm asking it specifically about aspect here.

Answer 1

I am not sure what you are looking for. The calculation of slope or aspect on a multi-segment linestring may be misunderstood due to the directional flexibility of each segment. However, the calculation is not so difficult when you get the 3D coordinates of the vertices. Here is an example:

Notes:

• I will assume you have a 3D LINESTRING and your data belongs to a projected CRS
• If you have 2D points you have to retrieve the third dimension from another resource, e.j. DEM

Edit after request of @user3386170

• For convenience I created a function to build a index for re-arrange the vertices of a line.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

# function to calculate slope  of several pair of vertices
slp <- function(x, y, z){ # x,y,z: numeric vectors with at least 2 values 
  dx <- diff(x)
  dy <- diff(y)
  dz <- diff(z)
  sqrt(rowSums(cbind(dx, dy) ^ 2)) / dz
}

# function to calculate aspect of several pair of vertices
asp <- function(x, y){ # x,y: numeric vectors with at least 2 values
  dx <- diff(x)
  dy <- diff(y)
  a <- atan2(dx, dy) * 180/pi
  if_else(a < 0, a + 360, a)
  }

# Smart build of index for pair vertex of a segment
build_index <- function(x) {
  len <- seq_along(x)
  sort(c(len[-length(len)], len[-1]))
}

# example of how the index of vertices in a line will be ordered
build_index(c(1:5)) 
#> [1] 1 2 2 3 3 4 4 5

# Example data: With  3D points (X, Y, Z)
l <- st_as_sfc(c("LINESTRING(760661.28 9814877.20 1000, 760978.04 9815259.80 1500, 761233.8 9815038.0 1200)",
                 "LINESTRING(760970.4 9814709.5 900, 760776.3 9814878.6 1000, 761391.8 9815341.6 1200)"), 
          crs = 32717) %>%
  st_sf()

# extract vertices coordinates
vertex <- st_coordinates(l) %>% 
  as_tibble() %>% 
  group_by(L1) %>%
  slice(build_index(L1)) %>%
  mutate(segment = ceiling(row_number() / 2)) # Create a segment index

# Calculate slope and aspect by segment
slp_asp <- vertex %>%
  group_by(L1, segment) %>% # important use Lines and segments indexes
  mutate(slope = c(NA, slp(X, Y, Z)),
         aspect = c(NA, asp(X, Y)))

# convert to segmented linestring
segment_lines <- slp_asp %>%
  st_as_sf(coords = c("X", "Y", "Z"), crs = st_crs(l)) %>%
  #let's convert to multipoint and aggregate the slope and aspect.
  summarise(slope = mean(slope, na.rm = TRUE), 
            aspect = mean(aspect, na.rm = TRUE), .groups = "drop") %>%
  st_cast("LINESTRING") %>% # convert to LINESTRING
  mutate(label = glue::glue("L{L1} - S{segment}"))

# Let's see what happened
segment_lines
#> Simple feature collection with 4 features and 5 fields
#> geometry type:  LINESTRING
#> dimension:      XYZ
#> bbox:           xmin: 760661.3 ymin: 9814710 xmax: 761391.8 ymax: 9815342
#> z_range:        zmin: 900 zmax: 1500
#> projected CRS:  WGS 84 / UTM zone 17S
#> # A tibble: 4 x 6
#>      L1 segment  slope aspect                                    geometry label 
#> * <dbl>   <dbl>  <dbl>  <dbl>                            <LINESTRING [m]> <glue>
#> 1     1       1  0.993   39.6 Z (760661.3 9814877 1000, 760978 9815260 1… L1 - …
#> 2     1       2 -1.13   131.  Z (760978 9815260 1500, 761233.8 9815038 1… L1 - …
#> 3     2       1  2.57   311.  Z (760776.3 9814879 1000, 760970.4 9814710… L2 - …
#> 4     2       2  3.85    53.0 Z (760776.3 9814879 1000, 761391.8 9815342… L2 - …

library(tmap)
s <- tm_shape(segment_lines) +
  tm_lines(col = "slope", style = "cat") +
  tm_text("label")
a <- tm_shape(segment_lines) +
  tm_lines(col = "aspect", style = "cat") +
  tm_text("label")
tmap_arrange(s, a)

^{Created on 2021-02-15 by the reprex package (v1.0.0)}

The output produced is a Linestring object, where each Line is a segment of a parent line (from original lines).

After calculate slope and aspect for your vertices, you could summarise / aggregate in a more representative value for your linestring. e.j mean slope and prevailing aspect. And finally join the result to your original spatial data frame

Best regards