[R-sig-Geo] How to plot great circles on a Robinson projection map in ggplot2?

[Tony Gray]

Hi! 

I want to plot great circles on a map in ggplot2 similar to the way it is done in this post: http://www.stanford.edu/~cengel/cgi-bin/anthrospace/great-circles-on-a-recentered-worldmap-in-ggplot -- though I am not so interested in re-centring. 

The wrinkle is that I want to do it on a Robinson projected map. I can use gcIntermediate on a Mercator map, but I can't figure it out when the coordinates have been translated using project().

The base code follows. Does anyone know how I can get great circles to connect the B locations to the A location? Any advice would be hugely appreciated.

Thank you! 

-- Tony

Here is the code:

library(rgdal)
library(ggplot2)
library(geosphere)

# B locations
lat <- c(36.0834857,40.96342,59.32893)
lng <- c(140.0766423,-76.6127329,18.06491)
lnglat <- cbind(lng,lat)
lnglat_df <- as.data.frame(lnglat)
names(lnglat_df) <- c('longitude','latitude')
lnglat_robin_df  <- project(cbind(lnglat_df$longitude, lnglat_df$latitude), proj="+init=ESRI:54030")
lnglat_robin_df <- as.data.frame(lnglat_robin_df)
names(lnglat_robin_df) <- c('longitude','latitude')

# A location
ALatitude <- 43.661102
ALongitude <- -79.395929
ACoords <- project(cbind(ALongitude, ALatitude), proj="+init=ESRI:54030")
ACoords_df <- as.data.frame(ACoords)
names(ACoords_df) <- c("longitude", "latitude")

# read shapefile
wmap <- readOGR(dsn="ne_110m_land", layer="ne_110m_land")
# project
wmap_robin <- spTransform(wmap, CRS("+init=ESRI:54030"))
# convert 
wmap_robin_df <- fortify(wmap_robin)

# create a blank ggplot theme -- thanks rpsychologist!

theme_opts <- list(theme(panel.grid.minor = element_blank(),
                    panel.grid.major = element_blank(),
                    panel.background = element_blank(),
                    plot.background = element_rect(fill="gray0"),
                    panel.border = element_blank(),
                    axis.line = element_blank(),
                    axis.text.x = element_blank(),
                    axis.text.y = element_blank(),
                    axis.ticks = element_blank(),
                    axis.title.x = element_blank(),
                    axis.title.y = element_blank(),
                    plot.title = element_text(size=22)))

# OK Plot dots
ggplot (wmap_robin_df, aes(long,lat, group=group)) +
geom_polygon(fill="gray10") +
geom_point(data=lnglat_robin_df, aes(longitude, latitude, group=NULL, fill=NULL), color="#32caf6", alpha=0.5, size=4) +
coord_equal() +
theme_opts +
scale_fill_manual(values=c("black", "black"), guide="none")

#save
ggsave
("map.png", width=12.5, height=8.25, dpi=72)