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