# 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.PolyLine([la, nyc, berlin, singapore, sydney],
weight=2,
color="red"

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

for start, end in pairwise([la, nyc, berlin, singapore, sydney]):

for start, end in pairwise([la, nyc, berlin, singapore, sydney]):

map

Out[11]:

## Building a geodesic PolyLine subclass¶

In [12]:
# not used anymore, but still nice ;-)
"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()
map

Out[14]:

Test using longer segment lengths:

In [15]:
map = folium.Map()
map

Out[15]:

Test using longer minimum length:

In [16]:
map = folium.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