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

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

• Could you elaborate on the final step? I don't understand how slope or aspect is assigned to each vertex. Shouldn't it be along a segment? (even if the segment is each line between vertices) Feb 16, 2021 at 0:04
• Please, see edits in the original response! Feb 16, 2021 at 3:01
• With your method, it doesn't distinguish between the start and end correctly so I get negative slopes (which I can fix) but also aspect which is oriented in the opposite direction of stream flow. I'm also not sure whether aspect defined as direction of line for each line segments obtains the general aspect of the stream segment (based on hill-slope) that I am looking for. I'm sorry if my original question was not sufficiently clear. Feb 16, 2021 at 23:00
• Oh, I thought it does. Could you provide a data example? If you see my example, the order of the vertices on the line represent the from/to of a segment. And of course I've used an example of up-slopes and down-slopes on `L1` and only up-slopes on L2. Both calculation, slope and aspect, are dependent of the order of vertices in a line. BTW, Negative slopes can be easily to positive to reverse the order, and as well aspect can be reversed with asp + 180 (if asp < 180) or asp - 180 (if asp > 180). Reversing means changing the order of start and end vertices Feb 17, 2021 at 2:05