<jupyter_start><jupyter_text>Using truly geodesic polylines with FoliumThis notebook shows how long straight lines are unusable on map projections like the very popular Mercator one, also used in [Folium](https://github.com/python-visualization/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](http://nbviewer.jupyter.org/github/deeplook/notebooks/blob/master/mapping/geodesic_polylines.ipynb).<jupyter_code>%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<jupyter_output><empty_output><jupyter_text>Straight lines<jupyter_code>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)
map<jupyter_output><empty_output><jupyter_text>Geodesic lines<jupyter_code>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)
for lat, lon in intermediatePoints(berlin, potsdam):
    print(lat, lon)<jupyter_output>52.516071 13.37698
52.39962 13.04784<jupyter_text>Mapping geodesic lines<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Building a geodesic PolyLine subclass<jupyter_code># 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<jupyter_output><empty_output><jupyter_text>Test using longer segment lengths:<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Test using longer minimum length:<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Let's check the distances (in km):<jupyter_code>Geodesic.WGS84.Inverse(*(la + berlin))["s12"] / 1000
Geodesic.WGS84.Inverse(*(berlin + sydney))["s12"] / 1000<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12899034892353378, 'fraction_numerical': 0.034521158129175944, 'mean_word_length': 5.778616352201258, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '151478', 'license': 'MIT', 'n_tokens_mistral': 2058, 'n_tokens_neox': 1881, 'n_words': 419, 'path': 'mapping/geodesic_polylines.ipynb', 'repo_name': 'deeplook/notebooks', 'size': '5388'}