2

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.

3
  • Would you be able to edit your question to verify that Geoquery is the GIS software that you are using, please?
    – PolyGeo
    Commented Sep 16, 2014 at 3:54
  • FYI: bost.ocks.org/mike/example
    – mdsumner
    Commented 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).
    – feverzsj
    Commented Sep 16, 2014 at 5:13

1 Answer 1

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.

2
  • well, I'm not familiar with R, it may take a while for me to figure out what this code actually doing
    – feverzsj
    Commented Sep 16, 2014 at 5:16
  • Yeah, sorry - just fwiw. If you want to use R I will expand this, just happened to be looking at something similar today.
    – mdsumner
    Commented Sep 16, 2014 at 11:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.