Drawing smooth great circle?

Question

I want to draw smooth great circle in any scale and any distance on a Mercator map. The key problem here is how many intermedia points do I need for different location and distance to draw it smoothly and not too slow.

Someone mentioned in an 2008 blog post:

For those who are really paying attention, you may guess this isn't exactly the code used in Geoquery 2008... for a start 10 divisions was too coarse so Geoquery uses a variable number of points based on the length of the line and its proximity to the poles (where the 'curve' is more pronounced). Geoquery ALSO needs to stop/start drawing lines that cross Longitude 180° (so they wrap nicely on Mercator, and join when projected on a globe)... which is done by calculating the intersection of two great circles because SQL Server 2008 couldn't do it (but that's another story...)

But I'm still not sure how to find proper parameters.

Answer 1

FWIW here is some R code to explore the minimum distance to use between generated points.

This uses great circle distance and intermediate points functions in geosphere, and reprojects to Mercator using rgdal. There's a hack to cut Antarctica a bit to stop it spreading to Infinity.

## background data
library(maptools)
data(wrld_simpl)
dodge <- function() {
    ## fix Antarctica for Mercator
    objid <- which(wrld_simpl$NAME == "Antarctica")
    ## first poly is the offender
    ## don't do this at home
    bad <- which(!wrld_simpl@polygons[[objid]]@Polygons[[1]]@coords[,2] > -90)
    w <- wrld_simpl
    w@polygons[[objid]]@Polygons[[1]]@coords[bad ,2] <- -85.46362
    w
}
library(rgdal)
prj <- "+proj=merc +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
wrld <- spTransform(dodge(), CRS(prj))
## tools for great circle distance and points
library(geosphere)

set.seed(91)

x <- randomCoordinates(2)
## distance in metres (WGS84)
dst <- distMeeus(x[1,,drop = FALSE], x[2,,drop = FALSE])

## choose your minimum (metres)
mindist <- 15e4
nn <- dst/mindist

## great circle points
gcpts <- gcIntermediate(x[1,,drop = FALSE], x[2,,drop = FALSE], n = nn, addStartEnd = TRUE)
## gc points on Mercator
library(rgdal)
gcmerc <- spTransform(SpatialPoints(gcpts, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")), CRS(prj))


plot(wrld)
plot(gcmerc, add = TRUE)

See vignette("geosphere") for more.