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.
%matplotlib inline
import math
import folium
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
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)
folium.PolyLine([la, nyc, berlin, singapore, sydney],
weight=2,
color="red"
).add_to(map)
<folium.features.PolyLine at 0x10fe8ab70>
map
from geographiclib.geodesic import Geodesic
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"]
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
for lat, lon in intermediatePoints(berlin, potsdam):
print(lat, lon)
52.516071 13.37698 52.39962 13.04784
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)
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
# 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}
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)
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
Test using longer segment lengths:
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
Test using longer minimum length:
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
Let's check the distances (in km):
Geodesic.WGS84.Inverse(*(la + berlin))["s12"] / 1000
9331.178262315017
Geodesic.WGS84.Inverse(*(berlin + sydney))["s12"] / 1000
16090.2938772803