Using truly geodesic polylines with Folium

This notebook shows how long straight lines are unusable on map projections like the very popular Mercator one, also used in Folium. And it implements a more detailed version with intermediate points in order to get a more realistic visualization of long lines on such map projections. These might look unusual (except for pilots), but they are much closer to the truth.

You will need to pip-install folium and geographiclib for this notebook to work as expected. If you look at this inside a GitHub repo you won't see the maps. Then try on the notebook viewer, there you will.



In [1]:

    
%matplotlib inline

import math

import folium



In [2]:

    
la  = 34.05351, -118.24531
nyc = 40.71453, -74.00713
berlin = 52.516071, 13.37698
potsdam = 52.39962, 13.04784
singapore = 1.29017, 103.852
sydney = -33.86971, 151.20711

Straight lines



In [3]:

    
map = folium.Map()



In [4]:

    
folium.Marker(location=la, popup="Los Angeles").add_to(map)
folium.Marker(location=nyc, popup="New York").add_to(map)
folium.Marker(location=berlin, popup="Berlin").add_to(map)
folium.Marker(location=singapore, popup="Singapore").add_to(map)
folium.Marker(location=sydney, popup="Sydney").add_to(map)

folium.PolyLine([la, nyc, berlin, singapore, sydney],
    weight=2,
    color="red"
).add_to(map)









    Out[4]:





<folium.features.PolyLine at 0x10fe8ab70>



In [5]:

    
map









    Out[5]:

Geodesic lines



In [6]:

    
from geographiclib.geodesic import Geodesic



In [7]:

    
def intermediatePoints(start, end, min_length_km=2000, segment_length_km=1000):
    geod = Geodesic.WGS84
    if geod.Inverse(*(start + end))["s12"] / 1000.0 < min_length_km:
        yield start
        yield end
    else:
        inv_line = geod.InverseLine(*(start + end))
        segment_length_m = 1000 * segment_length_km
        n = int(math.ceil(inv_line.s13 / segment_length_m))
        for i in range(n + 1):
            s = min(segment_length_m * i, inv_line.s13)
            g = inv_line.Position(s, Geodesic.STANDARD | Geodesic.LONG_UNROLL)
            yield g["lat2"], g["lon2"]



In [8]:

    
for lat, lon in intermediatePoints(berlin, sydney):
    print(lat, lon)









    



52.516071000000004 13.37698
53.98060697563462 28.186710757145015
53.54352761856573 43.36110532382015
51.269764128107006 57.605966972807366
47.45591983586865 70.09894541405228
42.47711273544207 80.67820310705217
36.66746005578033 89.58694414588376
30.28155147707397 97.19250070851366
23.50236612011881 103.84466410408318
16.461311189816005 109.8356359866446
9.256358051557065 115.40317997219232
1.9655955488759838 120.74640030861372
-5.342834815344675 126.04322573894493
-12.602793633961847 131.4671688473819
-19.742610059885617 137.20381613244894
-26.67641934166494 143.46796206741342
-33.292748740783914 150.52093336626385
-33.869710000000005 151.20710999999997



In [9]:

    
for lat, lon in intermediatePoints(berlin, potsdam):
    print(lat, lon)









    



52.516071 13.37698
52.39962 13.04784

Mapping geodesic lines



In [10]:

    
from itertools import tee

def pairwise(iterable):
    "s -> (s0,s1), (s1,s2), (s2, s3), ..."
    a, b = tee(iterable)
    next(b, None)
    return zip(a, b)



In [11]:

    
map = folium.Map()

folium.Marker(location=la, popup="Los Angeles").add_to(map)
folium.Marker(location=nyc, popup="New York").add_to(map)
folium.Marker(location=berlin, popup="Berlin").add_to(map)
folium.Marker(location=singapore, popup="Singapore").add_to(map)
folium.Marker(location=sydney, popup="Sydney").add_to(map)

for start, end in pairwise([la, nyc, berlin, singapore, sydney]):
    folium.PolyLine(list(intermediatePoints(start, end)), weight=2).add_to(map)

for start, end in pairwise([la, nyc, berlin, singapore, sydney]):
    folium.PolyLine([start, end], weight=2, color="red").add_to(map)

map









    Out[11]:

Building a geodesic PolyLine subclass



In [12]:

    
# not used anymore, but still nice ;-)
def extract_dict(aDict, keys):
    "Return a sub-dict of ``aDict`` by keys and remove those from ``aDict``."
    return {k: aDict.pop(k) for k in keys}



In [13]:

    
class GeodesicPolyLine(folium.PolyLine):
    """
    A geodesic version of a PolyLine inserting intermediate points when needed. 
    
    This will calculate intermediate points with some segment length whenever
    the geodesic length between two adjacent locations is above some maximal
    value.
    """
    def __init__(self, locations, min_length_km=2000, segment_length_km=1000, **kwargs):
        kwargs1 = dict(min_length_km=min_length_km, segment_length_km=segment_length_km)
        geodesic_locs = [intermediatePoints(start, end, **kwargs1) for start, end in pairwise(locations)]
        super().__init__(geodesic_locs, **kwargs)



In [14]:

    
map = folium.Map()
folium.Marker(location=la, popup="Los Angeles").add_to(map)
folium.Marker(location=berlin, popup="Berlin").add_to(map)
folium.Marker(location=sydney, popup="Sydney").add_to(map)
GeodesicPolyLine([la, berlin, sydney]).add_to(map)
map









    Out[14]:

Test using longer segment lengths:



In [15]:

    
map = folium.Map()
folium.Marker(location=la, popup="Los Angeles").add_to(map)
folium.Marker(location=berlin, popup="Berlin").add_to(map)
folium.Marker(location=sydney, popup="Sydney").add_to(map)
GeodesicPolyLine([la, berlin, sydney], segment_length_km=2000).add_to(map)
map









    Out[15]:

Test using longer minimum length:



In [16]:

    
map = folium.Map()
folium.Marker(location=la, popup="Los Angeles").add_to(map)
folium.Marker(location=berlin, popup="Berlin").add_to(map)
folium.Marker(location=sydney, popup="Sydney").add_to(map)
GeodesicPolyLine([la, berlin, sydney], min_length_km=10000).add_to(map)
map









    Out[16]:

Let's check the distances (in km):



In [17]:

    
Geodesic.WGS84.Inverse(*(la + berlin))["s12"] / 1000









    Out[17]:





9331.178262315017



In [18]:

    
Geodesic.WGS84.Inverse(*(berlin + sydney))["s12"] / 1000









    Out[18]:





16090.2938772803