{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Using truly geodesic polylines with Folium\n", "\n", "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.\n", "\n", "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). " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import math\n", "\n", "import folium" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "la = 34.05351, -118.24531\n", "nyc = 40.71453, -74.00713\n", "berlin = 52.516071, 13.37698\n", "potsdam = 52.39962, 13.04784\n", "singapore = 1.29017, 103.852\n", "sydney = -33.86971, 151.20711" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Straight lines" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "map = folium.Map()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "folium.Marker(location=la, popup=\"Los Angeles\").add_to(map)\n", "folium.Marker(location=nyc, popup=\"New York\").add_to(map)\n", "folium.Marker(location=berlin, popup=\"Berlin\").add_to(map)\n", "folium.Marker(location=singapore, popup=\"Singapore\").add_to(map)\n", "folium.Marker(location=sydney, popup=\"Sydney\").add_to(map)\n", "\n", "folium.PolyLine([la, nyc, berlin, singapore, sydney],\n", " weight=2,\n", " color=\"red\"\n", ").add_to(map)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Geodesic lines" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from geographiclib.geodesic import Geodesic" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def intermediatePoints(start, end, min_length_km=2000, segment_length_km=1000):\n", " geod = Geodesic.WGS84\n", " if geod.Inverse(*(start + end))[\"s12\"] / 1000.0 < min_length_km:\n", " yield start\n", " yield end\n", " else:\n", " inv_line = geod.InverseLine(*(start + end))\n", " segment_length_m = 1000 * segment_length_km\n", " n = int(math.ceil(inv_line.s13 / segment_length_m))\n", " for i in range(n + 1):\n", " s = min(segment_length_m * i, inv_line.s13)\n", " g = inv_line.Position(s, Geodesic.STANDARD | Geodesic.LONG_UNROLL)\n", " yield g[\"lat2\"], g[\"lon2\"]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "52.516071000000004 13.37698\n", "53.98060697563462 28.186710757145015\n", "53.54352761856573 43.36110532382015\n", "51.269764128107006 57.605966972807366\n", "47.45591983586865 70.09894541405228\n", "42.47711273544207 80.67820310705217\n", "36.66746005578033 89.58694414588376\n", "30.28155147707397 97.19250070851366\n", "23.50236612011881 103.84466410408318\n", "16.461311189816005 109.8356359866446\n", "9.256358051557065 115.40317997219232\n", "1.9655955488759838 120.74640030861372\n", "-5.342834815344675 126.04322573894493\n", "-12.602793633961847 131.4671688473819\n", "-19.742610059885617 137.20381613244894\n", "-26.67641934166494 143.46796206741342\n", "-33.292748740783914 150.52093336626385\n", "-33.869710000000005 151.20710999999997\n" ] } ], "source": [ "for lat, lon in intermediatePoints(berlin, sydney):\n", " print(lat, lon)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "52.516071 13.37698\n", "52.39962 13.04784\n" ] } ], "source": [ "for lat, lon in intermediatePoints(berlin, potsdam):\n", " print(lat, lon)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mapping geodesic lines" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "from itertools import tee\n", "\n", "def pairwise(iterable):\n", " \"s -> (s0,s1), (s1,s2), (s2, s3), ...\"\n", " a, b = tee(iterable)\n", " next(b, None)\n", " return zip(a, b)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "map = folium.Map()\n", "\n", "folium.Marker(location=la, popup=\"Los Angeles\").add_to(map)\n", "folium.Marker(location=nyc, popup=\"New York\").add_to(map)\n", "folium.Marker(location=berlin, popup=\"Berlin\").add_to(map)\n", "folium.Marker(location=singapore, popup=\"Singapore\").add_to(map)\n", "folium.Marker(location=sydney, popup=\"Sydney\").add_to(map)\n", "\n", "for start, end in pairwise([la, nyc, berlin, singapore, sydney]):\n", " folium.PolyLine(list(intermediatePoints(start, end)), weight=2).add_to(map)\n", "\n", "for start, end in pairwise([la, nyc, berlin, singapore, sydney]):\n", " folium.PolyLine([start, end], weight=2, color=\"red\").add_to(map)\n", "\n", "map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Building a geodesic PolyLine subclass" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# not used anymore, but still nice ;-)\n", "def extract_dict(aDict, keys):\n", " \"Return a sub-dict of ``aDict`` by keys and remove those from ``aDict``.\"\n", " return {k: aDict.pop(k) for k in keys}" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "class GeodesicPolyLine(folium.PolyLine):\n", " \"\"\"\n", " A geodesic version of a PolyLine inserting intermediate points when needed. \n", " \n", " This will calculate intermediate points with some segment length whenever\n", " the geodesic length between two adjacent locations is above some maximal\n", " value.\n", " \"\"\"\n", " def __init__(self, locations, min_length_km=2000, segment_length_km=1000, **kwargs):\n", " kwargs1 = dict(min_length_km=min_length_km, segment_length_km=segment_length_km)\n", " geodesic_locs = [intermediatePoints(start, end, **kwargs1) for start, end in pairwise(locations)]\n", " super().__init__(geodesic_locs, **kwargs)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "map = folium.Map()\n", "folium.Marker(location=la, popup=\"Los Angeles\").add_to(map)\n", "folium.Marker(location=berlin, popup=\"Berlin\").add_to(map)\n", "folium.Marker(location=sydney, popup=\"Sydney\").add_to(map)\n", "GeodesicPolyLine([la, berlin, sydney]).add_to(map)\n", "map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test using longer segment lengths:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "map = folium.Map()\n", "folium.Marker(location=la, popup=\"Los Angeles\").add_to(map)\n", "folium.Marker(location=berlin, popup=\"Berlin\").add_to(map)\n", "folium.Marker(location=sydney, popup=\"Sydney\").add_to(map)\n", "GeodesicPolyLine([la, berlin, sydney], segment_length_km=2000).add_to(map)\n", "map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test using longer minimum length:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "map = folium.Map()\n", "folium.Marker(location=la, popup=\"Los Angeles\").add_to(map)\n", "folium.Marker(location=berlin, popup=\"Berlin\").add_to(map)\n", "folium.Marker(location=sydney, popup=\"Sydney\").add_to(map)\n", "GeodesicPolyLine([la, berlin, sydney], min_length_km=10000).add_to(map)\n", "map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check the distances (in km):" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "9331.178262315017" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Geodesic.WGS84.Inverse(*(la + berlin))[\"s12\"] / 1000" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "16090.2938772803" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Geodesic.WGS84.Inverse(*(berlin + sydney))[\"s12\"] / 1000" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [default]", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.2" } }, "nbformat": 4, "nbformat_minor": 2 }