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

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