Create arrows along great circles between pairs of points

sf_arrow.R

sf_arrow <- function(from, to, len, deg, prefixFrom='from.', prefixTo='to.') {
  # from, to: sf objects. `from` should be either a single point, or the same
  #   number of points as in `to`.
  # len: the length of the arrow head, in map units.
  # deg: the angle (in degrees) of the arrow head.
  # prefixFrom, prefixTo: the attributes of `from` and `to` are added to the 
  #  returned sf object. These prefixes will be prepended to the existing names
  #  of `from` and `to`.
  crs <- st_crs(to)
  n <- nrow(to)
  names(from)[1:(length(from)-1)] <- paste0(prefixFrom, names(from)[1:(length(from)-1)])
  names(to)[1:(length(to)-1)] <- paste0(prefixTo, names(to)[1:(length(to)-1)])
  if(nrow(from) != n) {
    if(nrow(from) == 1) {
      from <- from[rep(1, n), ]
    } else {
      stop('Number of features in `from` must be either 1 or equal to number ',
           'of features in `to`.')
    }
  }
  
  lapply(seq_len(n), function(i) {
    from_ <- sf::st_transform(from[i, ], 4326)
    to_ <- sf::st_transform(to[i, ], 4326)
    stem <- geosphere::gcIntermediate(
      sf::as_Spatial(from_), sf::as_Spatial(to_), n=500, sp=TRUE, 
      breakAtDateLine=TRUE, addStartEnd=TRUE
    ) %>% sf::st_as_sfc() %>% 
      sf::st_transform(crs)
    
    xy <- sf::st_coordinates(tail(sf::st_cast(stem, 'POINT'), 2))
    x1 <- xy[1, 1]; x2 <- xy[2, 1]; y1 <- xy[1, 2]; y2 <- xy[2, 2]
    dirn <- x2 > x1
    theta <- atan((y2 - y1)/(x2 - x1))
    if(dirn) {
      arr1 <- c(x2 + len*cos(theta - (180-deg)*pi/180),
                y2 + len*sin(theta - (180-deg)*pi/180))
      arr2 <- c(x2 + len*cos(theta + (180-deg)*pi/180), 
                y2 + len*sin(theta + (180-deg)*pi/180))
    } else {
      arr1 <- c(x2 + len*cos(theta - deg*pi/180), 
                y2 + len*sin(theta - deg*pi/180))
      arr2 <- c(x2 + len*cos(theta + deg*pi/180), 
                y2 + len*sin(theta + deg*pi/180))
    }
    
    l1 <- sf::st_union(
      to[i, ]$geometry, 
      sf::st_as_sf(data.frame(x=arr1[1], y=arr1[2]), coords=c('x', 'y'), crs=sf::st_crs(from))) %>% 
      sf::st_cast('LINESTRING')
    l2 <- sf::st_union(
      to[i, ]$geometry, 
      sf::st_as_sf(data.frame(x=arr2[1], y=arr2[2]), coords=c('x', 'y'), crs=sf::st_crs(from))) %>% 
      sf::st_cast('LINESTRING')
    
    sf::st_as_sf(sf::st_geometry(sf::st_union(c(stem, l1, l2))))
  }) %>% 
    do.call(rbind, .) %>% 
    cbind(sf::st_drop_geometry(from), sf::st_drop_geometry(to))
}

Example

world <- rnaturalearth::ne_countries(returnclass='sf', scale=50) %>% 
  sf::st_transform('+proj=robin') %>% 
  sf::st_make_valid()

centroids <- sf::st_centroid(world)

from <- dplyr::filter(centroids, admin=='United Kingdom')

to <- dplyr::filter(centroids, admin != 'United Kingdom') %>% 
  .[sample(nrow(.), 6), ]

arrows <- sf_arrow(from, to, len=3e5, deg=35) %>% 
  sf::st_transform('+proj=robin') %>% 
  sf::st_make_valid() %>% 
  dplyr::mutate(volume=runif(nrow(.)))

o.mode <- tmap::tmap_mode('plot')
o.s2 <- sf::sf_use_s2(TRUE)

world %>% 
  tmap::tm_shape() +
  tmap::tm_polygons() +
  tmap::tm_shape(arrows) +
  tmap::tm_lines(lwd='volume', scale=4, legend.lwd.show=FALSE) +
  tmap::tm_style('bw')

tmap::tmap_mode(o.mode)
sf::sf_use_s2(o.s2)