# Drawing smooth great circle?

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.

• Would you be able to edit your question to verify that Geoquery is the GIS software that you are using, please?
– PolyGeo
Sep 16, 2014 at 3:54
• Sep 16, 2014 at 4:07
• @mdsumner : thanks for the amazing and awesome page, I will take a close look at this adaptive algorithm (if I can figure out where its related code is). Sep 16, 2014 at 5:13

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
w <- wrld_simpl
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])

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)