#!/usr/bin/env python
# coding: utf-8

# # 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](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). 

# In[1]:


get_ipython().run_line_magic('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)


# In[5]:


map


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


# In[9]:


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


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


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


# 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


# 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


# Let's check the distances (in km):

# In[17]:


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


# In[18]:


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