Erstellen Sie Großkreisbögen, die auf der Web Mercator-Karte gut aussehen?

Ich versuche, eine Karte mit Daten zu einigen Flügen zu erstellen und möchte Great Circle Arcs verwenden, um die Quelle und die Ziele zu verbinden.

Grundsätzlich möchte ich etwas Ähnliches wie die berühmte Facebook-Karte machen:

Ich habe die in diesem Beitrag angegebenen Funktionen verwendet: /gis//a/5205/442 (dh dieser Blog-Artikel: http://anitagraser.com/2011/08/20/visualizing-global-connections / ) und ich könnten Linien bekommen, aber sie überqueren die internationale Datumsgrenze und bündeln sich an den Polen:

@underdark hat im verlinkten Blogpost erwähnt, dass diese Zeilen geteilt werden müssen, aber ich weiß nicht, wie ich sie in PostGIS automatisch teilen soll.

Zusätzlich muss auch das Bündeln der Linien in der Nähe der Pole gelöst werden.

Wie mache ich beides, wenn ich die Punktpositionen für Start und Ende dieser Flüge habe?

Die Prinzipien in diesem Blog-Beitrag werden auf Allzweck-PostGIS übertragen.

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

Verwenden Sie ST_Segmentizeim Grunde genommen die Geografie und ein bisschen Magie, um Datumsgrenzen zu schneiden.

Sie könnten die Geodäten berechnen. Wenn Sie die Geodät von A nach B anzeigen möchten, können Sie zuerst die Entfernung und den Azimut von A nach B berechnen (inverses geodätisches Problem) und dann Punkte von A bis zu mehreren Punkten zwischen A und B berechnen (direktes geodätisches Problem). Ich habe ein einfaches Skript in Python hinzugefügt, das GeographicLib verwendet und nur das Zeug in GeoJSON ausgibt:

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))

Das Ergebnis ist die wahre Geodät zwischen den Punkten in WGS-84. Natürlich können Sie dann die Koordinaten in jede gewünschte Projektion umwandeln. Das auf geojson.io visualisierte Ergebnis sieht folgendermaßen aus: