Using loop to calculate bearings over a list

Question

My goal is to apply the geosphere::bearing function to a very large data frame, yet because the data frame concerns multiple individuals, I split it using the purrr package and split function.

I have seen the use of 'lists' and 'forloops' in the past but I have no experience with these.

Below is a fraction of my dataset, I have split the dataframe by ID, into a list with 43 elements. I have attached long and lat in wgs84 to the initial data frame.

ID         Date        Time         Datetime      Long        Lat   x            y
10_17   4/18/2017   15:02:00    4/18/2017 15:02 379800.5    5181001 -91.72272   46.35156
10_17   4/20/2017   6:00:00     4/20/2017 6:00  383409      5179885 -91.7044    46.34891
10_17   4/21/2017   21:02:00    4/21/2017 21:02 383191.2    5177960 -91.72297   46.35134
10_24   4/22/2017   10:03:00    4/22/2017 10:03 383448.6    5179918 -91.72298   46.35134
10_17   4/23/2017   12:01:00    4/23/2017 12:01 378582.5    5182110 -91.7242    46.34506
10_24   4/24/2017   1:00:00     4/24/2017 1:00  383647.4    5180009 -91.72515   46.34738
10_24   4/25/2017   16:01:00    4/25/2017 16:01 383407.9    5179872 -91.7184    46.32236
10_17   4/26/2017   18:02:00    4/26/2017 18:02 380691.9    5179353 -91.65361   46.34712
10_36   4/27/2017   20:00:00    4/27/2017 20:00 382521.9    5175266 -91.66127   46.3485
10_36   4/29/2017   11:01:00    4/29/2017 11:01 383443.8    5179909 -91.70303   46.35451
10_36   4/30/2017   0:00:00     4/30/2017 0:00  383060.8    5178361 -91.6685    46.32941
10_40   4/30/2017   13:02:00    4/30/2017 13:02 383426.3    5179873 -91.70263   46.35481
10_40   5/2/2017    17:02:00    5/2/2017 17:02  383393.7    5179883 -91.67099   46.34138
10_40   5/3/2017    6:01:00     5/3/2017 6:01   382875.8    5179376 -91.66324   46.34763
10_88   5/3/2017    19:02:00    5/3/2017 19:02  383264.3    5179948 -91.73075   46.3684
10_88   5/4/2017    8:01:00     5/4/2017 8:01   378554.4    5181966 -91.70413   46.35429
10_88   5/4/2017    21:03:00    5/4/2017 21:03  379830.5    5177232 -91.66452   46.37274

I then try this function

library(geosphere)
library(sf)
library(magrittr)

dis_list <- split(data, data$ID) 
 answer <- lapply(dis_list, function(df) {
  start <- df[-1 , c("x", "y")] %>% 
    st_as_sf(coords = c('x', 'y')) 
  
  end <- df[-nrow(df), c("x", "y")] %>% 
    st_as_sf(coords = c('x', 'y')) 
  
  
  angles <-geosphere::bearing(start, end) 
  
  df$angles <- c(NA, angles)
  df
})

answer

which gives the error

Error in .pointsToMatrix(p1) : 
'list' object cannot be coerced to type 'double'

Answer 1

A google search on "pass sf points to geosphere bearings" brings up this SE::GIS answer that seems to address the issue which I would characterize as "how to extract numeric vectors from items that are sf-classed POINTS": https://gis.stackexchange.com/questions/416316/compute-east-west-or-north-south-orientation-of-polylines-sf-linestring-in-r

I needed to work with a single section first and then apply the lessons from @Spacedman to this task:

> st_coordinates( st_as_sf(dis_list[[1]], coords = c('x', 'y')) )
          X        Y
1 -91.72272 46.35156
2 -91.70440 46.34891
3 -91.72297 46.35134
4 -91.72420 46.34506
5 -91.65361 46.34712

So st_coordinates wilL extract the POINTS classed values into a two column matrix that can THEN get passed to geosphere::bearings

 dis_list <- split(dat, dat$ID) 
 answer <- lapply(dis_list, function(df) {
     start <- df[-1 , c("x", "y")] %>% 
         st_as_sf(coords = c('x', 'y')) %>% st_coordinates
     
     end1 <- df[-nrow(df), c("x", "y")] %>% 
         st_as_sf(coords = c('x', 'y')) %>% st_coordinates
     
     
     angles <-geosphere::bearing(start, end1) 
     
     df$angles <- c(NA, angles)
     df
 })
 answer
#------------------------
$`10_17`
     ID      Date     Time      date  time     Long     Lat         x        y
1 10_17 4/18/2017 15:02:00 4/18/2017 15:02 379800.5 5181001 -91.72272 46.35156
2 10_17 4/20/2017  6:00:00 4/20/2017  6:00 383409.0 5179885 -91.70440 46.34891
3 10_17 4/21/2017 21:02:00 4/21/2017 21:02 383191.2 5177960 -91.72297 46.35134
5 10_17 4/23/2017 12:01:00 4/23/2017 12:01 378582.5 5182110 -91.72420 46.34506
8 10_17 4/26/2017 18:02:00 4/26/2017 18:02 380691.9 5179353 -91.65361 46.34712
         Datetime     angles
1 4/18/2017 15:02         NA
2  4/20/2017 6:00 -78.194383
3 4/21/2017 21:02 100.694352
5 4/23/2017 12:01   7.723513
8 4/26/2017 18:02 -92.387473

$`10_24`
     ID      Date     Time      date  time     Long     Lat         x        y
4 10_24 4/22/2017 10:03:00 4/22/2017 10:03 383448.6 5179918 -91.72298 46.35134
6 10_24 4/24/2017  1:00:00 4/24/2017  1:00 383647.4 5180009 -91.72515 46.34738
7 10_24 4/25/2017 16:01:00 4/25/2017 16:01 383407.9 5179872 -91.71840 46.32236
         Datetime    angles
4 4/22/2017 10:03        NA
6  4/24/2017 1:00  20.77910
7 4/25/2017 16:01 -10.58228

$`10_36`
      ID      Date     Time      date  time     Long     Lat         x        y
9  10_36 4/27/2017 20:00:00 4/27/2017 20:00 382521.9 5175266 -91.66127 46.34850
10 10_36 4/29/2017 11:01:00 4/29/2017 11:01 383443.8 5179909 -91.70303 46.35451
11 10_36 4/30/2017  0:00:00 4/30/2017  0:00 383060.8 5178361 -91.66850 46.32941
          Datetime    angles
9  4/27/2017 20:00        NA
10 4/29/2017 11:01 101.72602
11  4/30/2017 0:00 -43.60192

$`10_40`
      ID      Date     Time      date  time     Long     Lat         x        y
12 10_40 4/30/2017 13:02:00 4/30/2017 13:02 383426.3 5179873 -91.70263 46.35481
13 10_40  5/2/2017 17:02:00  5/2/2017 17:02 383393.7 5179883 -91.67099 46.34138
14 10_40  5/3/2017  6:01:00  5/3/2017  6:01 382875.8 5179376 -91.66324 46.34763
          Datetime     angles
12 4/30/2017 13:02         NA
13  5/2/2017 17:02  -58.48235
14   5/3/2017 6:01 -139.34297

$`10_88`
      ID     Date     Time     date  time     Long     Lat         x        y
15 10_88 5/3/2017 19:02:00 5/3/2017 19:02 383264.3 5179948 -91.73075 46.36840
16 10_88 5/4/2017  8:01:00 5/4/2017  8:01 378554.4 5181966 -91.70413 46.35429
17 10_88 5/4/2017 21:03:00 5/4/2017 21:03 379830.5 5177232 -91.66452 46.37274
         Datetime     angles
15 5/3/2017 19:02         NA
16  5/4/2017 8:01  -52.55217
17 5/4/2017 21:03 -123.91920

The help page for st_coordinates characterizes its function as "retrieve coordinates in matrix form".

Answer 2

Given the data is all ready in longitude and latitude form.
Then just using bearing(data[, c("Long", "Lat")]) and distGeo(data[, c("Long", "Lat")]) from geosphere on the split data frames will work. No need to create a start and end points.

library(geosphere)

dfs <- split(data, data$ID)
library(geosphere)
answer <- lapply(dfs, function(df) {
   df$distances <-c(distGeo(df[,c("Long", "Lat")]))
   df$bearings <- c(bearing(df[,c("Long", "Lat")]))
   
   df
})
answer

The sf package is useful for converting between coordinate systems, but with the data set above, that step can be skipped. I find the geosphere package more straight forward and simpler to use.