Search code examples
rspatialintersectiongeosphere

geosphere/R: Calculate intersection points between two great circles defined with 4 points


I'm trying to use geosphere to calculate intersection points between two great circles that are given in a dataframe format like this:

library(dplyr)
library(geosphere)

df <- data.frame(
  # first line
  ln1_lonb = 1:4,
  ln1_lone = 2:5,
  ln1_latb = 10:13,
  ln1_late = 11:14,
  # second line
  ln2_lonb = 5:8,
  ln2_lone = 6:9,
  ln2_latb = 14:17,
  ln2_late = 15:18
)

I tried using the gcintersect function from geosphere that takes matrices as inputs. In order to use it in a dataframe I used cbind, but seems mutate doesn't work well with this:

df %>% 
  mutate(
    int_points = gcIntersect(
      cbind(ln1_lonb, ln1_latb),
      cbind(ln1_lone, ln1_late),
      cbind(ln2_lonb, ln2_latb),
      cbind(ln2_lone, ln2_late)
    )
  )
>Error: Column `int_points` must be length 4 (the number of rows) or one, not 16

Probably the error is due to getting longer output than expected (the number of rows of the dataframe). So I tried putting it in a list:

df %>% 
  mutate(
    int_points = list(gcIntersect(
      cbind(ln1_lonb, ln1_latb),
      cbind(ln1_lone, ln1_late),
      cbind(ln2_lonb, ln2_latb),
      cbind(ln2_lone, ln2_late)
    ))
  )

Again here I see the output is all combinations instead of getting the 4 coordinates of the 2 points of intersections per row.

Expected output would be a list in a cell in a new column that will contain the coordinates of both points.

Is there a solution without using loops (or purrr) since this will be significantly slower than mutate.

The value for int_points in the first row should be the same as the output of this:

gcIntersect(cbind(1,2), cbind(10,11), cbind(5,6), cbind(14,15))

Solution

  • We can do a rowwise

    library(tidyverse)
    library(geosphere)
    df %>% 
          rowwise %>% 
          mutate( int_points = list(gcIntersect(
            cbind(ln1_lonb, ln1_latb),
            cbind(ln1_lone, ln1_late),
            cbind(ln2_lonb, ln2_latb),
            cbind(ln2_lone, ln2_late)
         ))) %>% 
         ungroup %>%
        mutate(int_points = map(int_points, as_tibble)) %>% 
        unnest(int_points)
    # A tibble: 4 x 12
    #  ln1_lonb ln1_lone ln1_latb ln1_late ln2_lonb ln2_lone ln2_latb ln2_late  lon1  lat1  lon2  lat2
    #     <int>    <int>    <int>    <int>    <int>    <int>    <int>    <int> <dbl> <dbl> <dbl> <dbl>
    #1        1        2       10       11        5        6       14       15 -176. -12.6  3.61  12.6
    #2        2        3       11       12        6        7       15       16 -175. -13.6  4.60  13.6
    #3        3        4       12       13        7        8       16       17 -174. -14.6  5.60  14.6
    #4        4        5       13       14        8        9       17       18 -173. -15.6  6.59  15.6