Compute the degrees between pairs of coordinates in R

Question

I am trying to calculate the angle (with 0 positioned west) between each pair of points.

Here below is a representative example of the dataset I have.

data<-data.frame(from_id=c("id1","id2","id3","id4"),to_id=c("id2","id5","id4","id1"), from_latitude=c(2.8033,3.1565,4.2032,3.5352), from_longitude=c(24.7305,26.8886,22.8083,25.8845), to_latitude=c(3.1565,0.0364,3.5352,2.8033), to_longitude=c(26.8886,18.2507,25.8845,24.7305))

I have tried to use stplanr::line_bearing(), however, I only get missing values and it is very demanding since I have many observations. Do you have any idea of how I could proceed with the data I have? Thank you!

Here below is the output I would like (the degrees are wrong but it is just to illustrate)

 from_id to_id from_latitude from_longitude to_latitude to_longitude degrees
1     id1   id2        2.8033        24.7305      3.1565      26.8886      93
2     id2   id5        3.1565        26.8886      0.0364      18.2507     352
3     id3   id4        4.2032        22.8083      3.5352      25.8845      53
4     id4   id1        3.5352        25.8845      2.8033      24.7305     182

Here below my code in case:

# Load Libraries
library(sf)
library(sp)

# I transform the points into Lines
line_data<-data
b = line_data[, c("from_longitude", "from_latitude")]
names(b) = c("long", "lat")
e = line_data[, c("to_longitude", "to_latitude")]
names(e) = c("long", "lat")

line_data $geometry = do.call(
  "c", 
  lapply(seq(nrow(b)), function(i) {
    st_sfc(
      st_linestring(
        as.matrix(
          rbind(b[i, ], e[i, ])
        )
      ),
      crs = 4326
    )
  }))

data_sf = st_as_sf(line_data)

# I create I shapefile of starting points with the "from" coordinates
coordinates(data)=~from_longitude+from_latitude
data <-st_as_sf(data)
data <-data%>%sf::st_set_crs(4326)

crop_buffer <- st_buffer(data, 1000)
crop_buffer <- st_crop(data_sf, crop_buffer)
crop_buffer <- st_intersection(data_sf, crop_buffer)

bearing<-stplanr::line_bearing(crop_buffer)

Interesting that you calculate zero from west... Usual definition of azimuth is from north. But I digress... Apart from the solution by @tospig consider two approaches from answers to this earlier question: https://stackoverflow.com/questions/72276233/calculate-the-angle-between-two-sf-points-in-r-to-calculate-the-direction-of-roa — Jindra Lacko, May 23 '23 at 11:28

tospig · Accepted Answer · 2023-05-23T10:29:58.887

You could write your own function to calculate the bearing

data <- data.frame(
  from_id = c("id1","id2","id3","id4")
  , to_id = c("id2","id5","id4","id1")
  , from_latitude = c(2.8033,3.1565,4.2032,3.5352)
  , from_longitude = c(24.7305,26.8886,22.8083,25.8845)
  , to_latitude = c(3.1565,0.0364,3.5352,2.8033)
  , to_longitude = c(26.8886,18.2507,25.8845,24.7305)
  )

bearing <- function(aLon, aLat, bLon, bLat) {
  
  y <- sin(bLon - aLon) * cos(bLat);
  x <- (cos(aLat) * sin(bLat)) - (sin(aLat) * cos(bLat) * cos(bLon - aLon));
  
  ## constrain to a compass bearing, 0 == north
  b <- (atan2(y, x) * 180 / pi) + 360;
  b <- b %% 360
  
  return(b) 
}

data$bearing <- bearing(
  aLon = data$from_longitude
  , aLat = data$from_latitude
  , bLon = data$to_longitude
  , bLat = data$to_latitude
  )

# > data
#   from_id to_id from_latitude from_longitude to_latitude to_longitude  bearing
# 1     id1   id2        2.8033        24.7305      3.1565      26.8886 258.4687
# 2     id2   id5        3.1565        26.8886      0.0364      18.2507 266.2079
# 3     id3   id4        4.2032        22.8083      3.5352      25.8845 356.5173
# 4     id4   id1        3.5352        25.8845      2.8033      24.7305 117.7054

reference - https://www.movable-type.co.uk/scripts/latlong.html

Thank you very much! it works and it is way more efficient :) — Coralie, May 23 '23 at 10:21

Compute the degrees between pairs of coordinates in R

1 Answers1