Skip to content

Instantly share code, notes, and snippets.

@johnbaums
Last active November 7, 2022 03:03
Show Gist options
  • Save johnbaums/24478656055c48e24bc4393ac782a4e3 to your computer and use it in GitHub Desktop.
Save johnbaums/24478656055c48e24bc4393ac782a4e3 to your computer and use it in GitHub Desktop.
Create arrows along great circles between pairs of points
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))
}
@johnbaums
Copy link
Author

johnbaums commented Nov 7, 2022

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)

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment