Making Great Circle Arcs which look good on Web Mercator map?

Question

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:

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:

@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?

I have given some examples here: http://gis.stackexchange.com/questions/133026/great-circles-in-qgis-and-export-in-3857-webmap. Note that eqdc does not deliver true great-circles, but aeqd does. — AndreJ, Mar 14 '16 at 16:53

sema · Answer 1 · 2016-03-15T07:18:06.810

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):
                linestrings.append(coordinates)
                coordinates = []
        coordinates.append((direct['lon2'], direct['lat2']))

    linestrings.append(coordinates)
    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']:
    linestrings.append(linestring)

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

print(MultiLineString(linestrings))

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:

See http://geographiclib.sourceforge.net/html/python/examples.html#computing-waypoints for an alternative (and faster!) way of doing this with GeographicLib. Note also the used of the LONG_UNROLL flag to ensure continuity of longitudes. — cffk, Mar 17 '16 at 17:14

score 6 · Accepted Answer · answered Mar 14 '16 at 18:54

6

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

http://blog.cartodb.com/jets-and-datelines/

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

answered Mar 14 '16 at 18:54

Paul Ramsey

19,865
1
47
57

Making Great Circle Arcs which look good on Web Mercator map?

2 Answers2

Linked