How to plot arcs between two cities on a globe map using ggplot and sf-CodePudding

I have seen R codes generating a global map with ggplot and sf library, but I couldn't understand many of the details as I have limited knowledge in geography. Here is a visualization task I need to do with R.

If I have the coordinates (e.g. latitude and longitude) of two cities in the world, and I want to draw a red arc between them, which is the shorter arc on the great circle passing the two cities, how could I do that? The result should look like what we will get when we look up airline routes on google map.

Here is where I stopped at a fixed global map plot:

library("ggplot2")
library("sf")
world <- ne_countries(scale = "medium", returnclass = "sf")
ggplot(data = world)  
    geom_sf()  
    coord_sf(crs = "some crs code")

The crs part requires specific knowledge in geography, I don't know what to use and the meaning of different codes.

CodePudding user response：

You may want to have a look to https://geocompr.robinlovelace.net/ for getting an insight of what is spatial data.
The code is the reference to the map projection. The world is a sphere(-ish) and the representation on a two-dimension surface (screen, paper, etc.) is a challenge. That's why there are so many different projections (see more on https://geocompr.robinlovelace.net/spatial-class.html#crs-intro and https://geocompr.robinlovelace.net/reproj-geo-data.html#crs-in-r). Check also https://epsg.io/ to find the codes of each projection.

Now, getting back to your question, you may want to plot "Great Circles" between cities, that is the curve that represent the shortest path between points on a sphere (flights uses these lines for optimizing the flight time and fuel):

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(rnaturalearth)
library(tidyverse)
world <- ne_countries(scale = "medium", returnclass = "sf")

# Get three cities
cities <- ne_download(type = "populated_places", returnclass = "sf") %>%
  filter(NAME %in% c("San Francisco", "Rio de Janeiro", "Oslo"))

# On Mercator
ggplot(world)  
  geom_sf()  
  geom_sf(data=cities, color="red", size=3)  
  coord_sf(crs="EPSG:3857")

# On Robinson projection
ggplot(world)  
  geom_sf()  
  geom_sf(data=cities, color="red", size=3)  
  coord_sf(crs="ESRI:54030")


# Create great circles
library(geosphere)

init <- cities %>% as("Spatial")
dest <- bind_rows(cities[-1,], cities[1,]) %>% as("Spatial")

lines <- gcIntermediate(
  init,
  dest,
  sp = TRUE,
  addStartEnd = TRUE,
  n = 100,
  breakAtDateLine = FALSE
)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in Proj4
#> definition

lines_sf <- st_union(st_as_sf(lines))

# Now we have the lines
ggplot(lines_sf)  
  geom_sf()  
  geom_sf(data=cities, color="red")

# Create a nice plot

ggplot(world)  
  geom_sf()  
  geom_sf(data=lines_sf)  
  geom_sf(data=cities, color="red", size=3)  
  coord_sf(crs=" proj=laea  x_0=0  y_0=0  lon_0=-74  lat_0=40")

^{Created on 2022-06-10 by the reprex package (v2.0.1)}