I'm trying to create a map, with data indicating some flights, and want to use Great Circle Arcs, to connect the source and destinations.

Basically, I want to do something similar to the famous Facebook map: enter image description here

I used the functions given in this post: https://gis.stackexchange.com/a/5205/442, (i.e. this blog article: http://anitagraser.com/2011/08/20/visualizing-global-connections/) and I could get lines, but they cross the International Date Line, as well as bunch up, at the poles:

enter image description here

@underdark, has mentioned in the linked blogpost that these lines need need to be split, but I don't know how to split them automatically in PostGIS.

Additionally the bunching of the lines near the poles need to be resolved as well.

How do I do both of these, when I have the point locations for start and end of these flights?


The principles in this blog post transfer over to general purpose PostGIS.


Basically, use ST_Segmentize on geography, and a bit of magic to slice date-line crossing lines.

| improve this answer | |

You could calculate the Geodesics. Saying you want to show the geodesic from A to B, you could first calculate the distance and azimuth from A to B (inverse Geodesic problem) and then calculate points from A to several points between A and B (direct Geodesic problem). I have added a simple script in Python using GeographicLib just outputting the stuff in GeoJSON:

from geographiclib.geodesic import Geodesic
from geojson import MultiLineString

def geodesic(lat1, lon1, lat2, lon2, steps):
    inverse = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2)
    linestrings = []
    coordinates = []

    for i in range(0, steps + 1):
        direct = Geodesic.WGS84.Direct(inverse['lat1'], inverse['lon1'], inverse['azi1'], (i / float(steps)) * inverse['s12'])
        if len(coordinates) > 0:
            if (coordinates[-1][0] < -90 and direct['lon2'] > 90) or (coordinates[-1][0] > 90 and direct['lon2'] < -90):
                coordinates = []
        coordinates.append((direct['lon2'], direct['lat2']))

    geojson = MultiLineString(linestrings)
    return geojson

linestrings = []

# San Francisco: 37.7793, -122.4192
# Bangalore: 12.9, 77.616667
for linestring in geodesic(37.7793, -122.4192, 12.95, 77.616667, 100)['coordinates']:

# Boston: 42.357778, -71.059444
# Bangalore: 12.9, 77.616667
for linestring in geodesic(42.357778, -71.059444, 12.95, 77.616667, 100)['coordinates']:


The result is the true geodesic between the points in WGS-84. Of course, you could then transform the coordinates to whatever projection you need. The result visualized on geojson.io looks like this:

enter image description here

| improve this answer | |

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.