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.

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

  • well, I'm not familiar with R, it may take a while for me to figure out what this code actually doing – feverzsj Sep 16 '14 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 Sep 16 '14 at 11:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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