Calculate the angle between two sf points in R to calculate the direction of road line strings

Question

Trying to calculate the angle of an sf vector of linestrings using the start and end points of each sf linestring. I don't want to split the line up.

I used the formula in this similar example for the degrees calculation Calculate angle between two Latitude/Longitude points

Example below downloads a small road network, creates a function 'angles' which adds a column to the data frame with the angle and then plots it in mapview so each line is coloured by angle. However, the values are not correct. I think simple error with the formula?

library(osmdata)
library(sf)
library(dplyr)


bb <- st_bbox(c(xmin = 5.22, xmax = 5.246, ymax = 52.237, ymin = 52.227), crs = st_crs(4326))

x <- opq(bbox = bb) %>%
  add_osm_feature(key = c('highway')) %>%
  osmdata_sf()

##extract building polygons
roads <- x$osm_lines %>% 
  filter(highway == "motorway") %>% 
  select(osm_id, geometry)

angles <- function(x){
  
  ## define line finish coordinates
  f <- x %>%
    st_transform(28992) %>% 
    st_line_sample(sample = 1) %>%
    st_cast("POINT") %>% 
    st_transform(4326) %>% 
    st_coordinates()
  
  ## define line start coordinates
  s <- x %>%
    st_transform(28992) %>% 
    st_line_sample(sample = 0) %>% 
    st_cast("POINT") %>% 
    st_transform(4326) %>% 
    st_coordinates()
  
  ## get latitude of finish points
  lat2 <- f[,2]
  ## get latitudes of start points
  lat1 <- s[,2]
  

  
  ## get longitudes of start points
  lon2 <- f[,1]
  ## get longitudes of start points
  lon1 <- s[,1]
  
  ## delta longitudes
  #dlon <- f[,1]-s[,1]
  #theta <- (atan2(sin(dlon)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon)))*180/pi
  
  theta <- atan2(lat2-lat1, lon2-lon1)*180/pi
  
  x$angle_deg <- theta
  
  return(x)
  
}

a <- angles(roads) %>% 
  select(angle_deg)

library(mapview)

mapview(a)

Answer 1

I suggest you consider lwgeom::st_geod_azimuth() - it will give bearing (in radians).

Two things to consider:

• it requires geometry type POINT - because bearing from one point to another is crystal clear, but slope of a LINESTRING may / will vary one section to the next
• it will give a vector of one element less than the original points - the last point has no defined bearing, as there is no "next" point; to resolve this I am adding artificial NA to the last element.

Should you desire output in degrees just pipe the calculation to a units::set_units("degrees") call.

library(osmdata)
library(sf)
library(dplyr)


bb <- st_bbox(c(xmin = 5.22, xmax = 5.246, ymax = 52.237, ymin = 52.227), crs = st_crs(4326))

x <- opq(bbox = bb) %>%
  add_osm_feature(key = c('highway')) %>%
  osmdata_sf()

##extract building polygons
roads <- x$osm_lines %>% 
  filter(highway == "motorway") %>% 
  select(osm_id, geometry)  %>% 
  st_cast("POINT") %>% 
  mutate(bearing = c(lwgeom::st_geod_azimuth(.),
                     units::set_units(NA, "radians"))) # 1 extra point align vector length

mapview::mapview(roads["bearing"])

Answer 2

Many thanks for this. I realise maybe I wasn't completely clear in the question. I don't want to split the linestrings up as I am using the angle to populate these lines with data from other geometry. However, your code helped a lot. Here is the final code I have which adapts yours to result in angle of each unique linestring.

library(osmdata)
library(sf)
library(dplyr)

    bb <- st_bbox(c(xmin = 5.22, xmax = 5.246, ymax = 52.237, ymin = 52.227), crs = st_crs(4326))
    
    x <- opq(bbox = bb) %>%
      add_osm_feature(key = c('highway')) %>%
      osmdata_sf()
    
pts <- x$osm_lines %>% ## extract polylines from osm data
  filter(highway == "motorway") %>% ## only motorway network
  select(osm_id, geometry)  %>% ## remove excess columns
  st_cast("POINT") %>% ## split into constituent points
  group_by(osm_id) %>% 
  filter(row_number()==1 | row_number()==n()) %>% ## take start and end points of each osm id
  mutate(bearing = c(lwgeom::st_geod_azimuth(geometry))) %>% ## calculate bearing of start and end points
  mutate(bearing = units::set_units(bearing, "degrees")) %>% ## convert to degrees
  distinct(osm_id, .keep_all = TRUE) ## just take one point for each osm id
    
    roads <- x$osm_lines %>% 
      filter(highway == "motorway") %>% 
      select(osm_id, geometry)  %>% 
      mutate(bearing = pts$bearing)
    
    mapview::mapview(roads["bearing"])

Answer 3

You might be interested in using the sfnetworks package for spatial networks. It includes the function edge_azimuth() to calculate the azimuth of the edges in a network.

Do note that the OSM data might need some pre-processing before they form a clean, routable network. See e.g. the sfnetworks vignette on pre-processing and cleaning: https://luukvdmeer.github.io/sfnetworks/articles/sfn02_preprocess_clean.html

Reprex:

library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(sf)
#> Linking to GEOS 3.10.1, GDAL 3.4.0, PROJ 8.2.0; sf_use_s2() is TRUE
library(dplyr)
library(sfnetworks)

bb <- st_bbox(c(xmin = 5.22, xmax = 5.246, ymax = 52.237, ymin = 52.227), crs = st_crs(4326))

x <- opq(bbox = bb) %>%
  add_osm_feature(key = c('highway')) %>%
  osmdata_sf()

net <- x$osm_lines %>% 
  filter(highway == "motorway") %>% 
  as_sfnetwork() %>%
  activate("edges") %>%
  mutate(bearing = edge_azimuth()) %>%
  mutate(bearing = units::set_units(bearing, "degrees"))

plot(st_as_sf(net)[, "bearing"])

^{Created on 2022-08-10 by the reprex package (v2.0.1)}