Recently, I enjoyed playing around with mapnificent to see their definition of areas one can reach with public transport in a given time. However, I was wanted to (1) expand this to include more cities and (2) create an illustration where one can compare the differences easily.

Public Transport / Drive (15 minutes) of World Cities

World Cities - 15min Public Transport

I will briefly go through the method used.

Step 1 (Generating the Data):

I used the code I described in an earlier post to generate a 15-minute drive-time and public-transport isochrone using 2PM (local-time) for around 30 cities (the centre point was chosen automatically the Google API as I simply passed the city name to the geocoder).

For example, this is how London looked like:

london

We can compare this to a screenshot of Mapnificient’s rendering:

london map.PNG

I thought they looked pretty similar – bear in mind that they use a different data-source and public transport varies by time of day / day of the week, etc.

And Barcelona:

barce

Ooops! One of the nice things about the Geocoding/Distance API is that if a point is not available it selects the closest available point. This means that isochrones are smoother and look less jerky (in general), however this doesn’t quite work for the ocean! I had two options:

A – Geocode each co-ordinate and then replace the initial co-ordinate with the geocoded one (this would move all the sea points to land), however this would also eliminate the useful ‘smoothing’ feature

B – Find some water shape-files and drop the points not in land using a point-in-polygon test

Step 2 (Cleaning the Data):

I found some water here and here (actually lakes). I then wrapped shapely’s “shape.contains(point)” function and used that to kill the points in the water. I also took advantage of the re-run (on the .csv of points output in the previous step) to tweak the map how I wanted (constant zoom to allow comparison, I wanted one version in black to emphasise the shape and another in colour, some useful statistics, etc).

Ideally one would use a spatial-index (like “r-tree”) for the point in polygon test but I was in a slight rush and it didn’t take that long anyway (when I get some more time I will revisit this).

So, here is Barcelona again:

barce.PNG

The stats created by the javascript:

statas

And then the version I used in my final illustration:

barce

Python Scripts:

import fiona
from shapely.geometry import Point, shape

# Shape files (Small Lakes Europe, Lakes, Coast)
le_loc = ".../ne_10m_lakes_europe.shp"
lw_loc = ".../ne_10m_lakes.shp"
c_loc = ".../water_polygons.shp"

print("Opening Shape-Files")
le_poly = [shape(pol['geometry']) for pol in fiona.open(le_loc)]
lw_poly = [shape(pol['geometry']) for pol in fiona.open(lw_loc)]
c_poly = [shape(pol['geometry']) for pol in fiona.open(c_loc)]

def gen_shapes():
    for one in [c_poly, lw_poly, le_poly]:
        for shape in one:
            yield shape
            
def point_in_water(lat, lon):
    xy_point = Point(lon, lat) 
    for shape in gen_shapes():
        if shape.contains(xy_point):
            return True
    return False

def html_isochrone(coord_set,
                   radius_km,
                   map_name,
                   origin=["51", "0"],
                   black = False):
    """
    Create a JS map of the isochrones
    """
    htmltext = """<!DOCTYPE html >
      <style type="text/css">
                html, body {
                    height: 100%;
                    width: 100%;
                    padding: 0px;
                    margin: 0px;
                }
    </style>
    <head>
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no" />
    <meta http-equiv="content-type" content="text/html; charset=UTF-8"/>
    <title>Isochrone Analysis</title>
    <xml id="myxml">
    <markers>
    """
    for coord in coord_set:
        # print(coord)
        rowcord = '<marker name = "' + coord[2] + '" lat = "' + coord[0] + '" lng = "' + coord[1] + '"/>\n'
        htmltext += rowcord
    htmltext += """
    </markers>
    </xml>
    <script type="text/javascript" src="https://maps.googleapis.com/maps/api/js?&libraries=geometry"></script>
    <script type="text/javascript">
    var XML = document.getElementById("myxml");
    if(XML.documentElement == null)
    XML.documentElement = XML.firstChild;
    var MARKERS = XML.getElementsByTagName("marker");
    """
    htmltext += "var RADIUS_KM = " + str(radius_km) + ";"
    htmltext += """
    var map;
    var geocoder = new google.maps.Geocoder();
    var BUFFER = 1.1;
    var circles_transit = [];
    var circles_driving = [];
    var circles_walking = [];
    function drawCircle(point, radius, dir){
        // Function from stackoverflow
        var d2r = Math.PI / 180;   // degrees to radians
        var r2d = 180 / Math.PI;   // radians to degrees
        var earthsradius = 3963; // 3963 is the radius of the earth in miles
        var points = 32;
        // find the radius in lat/lon
        var rlat = (radius / earthsradius) * r2d;
        var rlng = rlat / Math.cos(point.lat() * d2r);

        var extp = new Array();
        if (dir==1) {var start=0;var end=points+1} // one extra here makes sure we connect the
        else{var start=points+1;var end=0}
        for (var i=start; (dir==1 ? i < end : i > end); i=i+dir)
        {
            var theta = Math.PI * (i / (points/2));
            ey = point.lng() + (rlng * Math.cos(theta)); // center a + radius x * cos(theta)
            ex = point.lat() + (rlat * Math.sin(theta)); // center b + radius y * sin(theta)
            extp.push(new google.maps.LatLng(ex, ey));
        }
        return extp;
    }"""
    htmltext += """
    function load() {{
        // Initialize
        var my_lat = {0};
        var my_lng = {1};""".format(origin[0], origin[1])
    htmltext += """
        var orig =  new google.maps.LatLng(my_lat,my_lng);
        var max_transit = 0;
        var num_transit = 0;
        var max_driving = 0;
        var num_driving = 0;
        var max_walking = 0;
        var num_walking = 0;"""
    if not black:
        htmltext += """
            // Custom formatting
            var mapOptions = {
                    center: orig,
                    zoom: 12,
                    backgroundColor: 'black',
                    styles: [
                    //{"featureType":"all", "stylers":[{"visibility":"off"}]},
                    //{"featureType":"road","stylers":[{"saturation":0}]},
                    ]
            };"""
    else:
        htmltext += """
            // Custom formatting
            var mapOptions = {
                    center: orig,
                    zoom: 12,
                    backgroundColor: 'black',
                    styles: [
                    {"featureType":"all", "stylers":[{"visibility":"off"}]},
                    {"featureType":"road","stylers":[{"saturation":0}]}
                    ]
            };"""      
    htmltext += """
        map = new google.maps.Map(document.getElementById('map'),
            mapOptions);
        var bounds = new google.maps.LatLngBounds();
        for (var i = 0; i < MARKERS.length; i++) {
            var name = MARKERS[i].getAttribute("name");
            if (name == 'transit') {
                num_transit += 1;
                var point_i = new google.maps.LatLng(
                    parseFloat(MARKERS[i].getAttribute("lat")),
                    parseFloat(MARKERS[i].getAttribute("lng")));
                circles_transit.push(drawCircle(point_i,RADIUS_KM*(1000/1609.344)*BUFFER,1))
                bounds.extend(point_i);
                // Compute distance
                var distance = google.maps.geometry.spherical.computeDistanceBetween (orig, point_i);
                if (distance > max_transit) {max_transit = distance};
            }
            if (name == 'driving') {
                num_driving += 1;
                var point_i = new google.maps.LatLng(
                    parseFloat(MARKERS[i].getAttribute("lat")),
                    parseFloat(MARKERS[i].getAttribute("lng")));
                circles_driving.push(drawCircle(point_i,RADIUS_KM*(1000/1609.344)*BUFFER,1))
                bounds.extend(point_i);
                // Compute distance
                var distance = google.maps.geometry.spherical.computeDistanceBetween (orig, point_i);
                if (distance > max_driving) {max_driving = distance};
            }
            if (name == 'walking') {
                num_walking += 1;
                var point_i = new google.maps.LatLng(
                    parseFloat(MARKERS[i].getAttribute("lat")),
                    parseFloat(MARKERS[i].getAttribute("lng")));
                circles_walking.push(drawCircle(point_i,RADIUS_KM*(1000/1609.344)*BUFFER,1))
                bounds.extend(point_i);
                // Compute distance
                var distance = google.maps.geometry.spherical.computeDistanceBetween (orig, point_i);
                if (distance > max_walking) {max_walking = distance};
            }
        };"""
    htmltext += """
        // Colours
        var col_transit = "{0}";
        var col_driving = "{1}";
        var col_walking = "{2}";
    """.format("red", "blue", "green")
    if not black:
        htmltext += """
            // Maximum Bound Circles
            if (num_driving > 0) {
                console.log("Driving (metres)", max_driving)
                console.log("Driving (count)", num_driving)
                var cityCircle = new google.maps.Circle({
                    strokeColor: col_driving,
                    strokeOpacity: 0.8,
                    strokeWeight: 1,
                    fillOpacity: 0,
                    map: map,
                    center: orig,
                    radius: max_driving
                });	
            }		
            if (num_transit > 0) {
                console.log("Transit (metres)", max_transit)
                console.log("Transit (count)", num_transit)
                var cityCircle = new google.maps.Circle({
                    strokeColor: col_transit,
                    strokeOpacity: 0.8,
                    strokeWeight: 1,
                    fillOpacity: 0,
                    map: map,
                    center: orig,
                    radius: max_transit
                });		
            }
            if (num_walking > 0) {
                console.log("Walking (metres)", max_walking)
                console.log("Walking (count)", num_walking)
                var cityCircle = new google.maps.Circle({
                    strokeColor: col_walking,
                    strokeOpacity: 0.8,
                    strokeWeight: 1,
                    fillOpacity: 0,
                    map: map,
                    center: orig,
                    radius: max_walking
                });		
            }"""
    htmltext += """
        // Isochrone polygons
        if (circles_driving.length > 0) {
            var joined_driving = new google.maps.Polygon({
                    paths: circles_driving,
                    strokeColor: col_driving,
                    strokeWeight: 0,
                    fillColor: col_driving,
                    fillOpacity: 0.3
            });
            joined_driving.setMap(map);
        }
        if (circles_transit.length > 0) {
            var joined_transit = new google.maps.Polygon({
                    paths: circles_transit,
                    strokeColor: col_transit,
                    strokeWeight: 0,
                    fillColor: col_transit,
                    fillOpacity: 0.3
            });
            joined_transit.setMap(map);
        }
        if (circles_walking.length > 0) {
            var joined_walking = new google.maps.Polygon({
                    paths: circles_walking,
                    strokeColor: col_walking,
                    strokeWeight: 0,
                    fillColor: col_walking,
                    fillOpacity: 0.3
            });
            joined_walking.setMap(map);
        }
        //Scaling and zoom
        map.fitBounds(bounds);
        var listener = google.maps.event.addListener(map, "idle", function() { 
            map.setZoom(12); 
            google.maps.event.removeListener(listener); 
        });
    }
    </script>
    </head>
    <body onload="load()">
    <center>
    <div style="padding-top: 20px; padding-bottom: 20px;">
    <div id="map" style="width:90%; height:1024px;"></div>
    </center>
    </body>
    </html>
    """
    with open('%s_land_on_water.html' % map_name, 'w') as f:
        f.write(htmltext)
    f.close()
import csv
import os
import time
from multiprocessing import Pool
from geopy.geocoders import GoogleV3

geolocator = GoogleV3()
path_2_csvs = ".../Google_Isochrone/"
duration = 15
radius_km = 0.075
counter = 0

for i in os.listdir(path_2_csvs):
    if i.endswith("_15min.csv"):
        
        # A: GEOCODE CITY
        centre = geolocator.geocode(i.replace("_15min.csv",""))
        centre_coord = ["%.5f" % centre.latitude , "%.5f" % centre.longitude]
        print("Num: %d, Name: %s" % (counter, i))
        print("Geocoded:", centre, centre_coord)
        counter += 1
        time.sleep(1)
        
        # B: REMOVE WATER FROM CSV
        in_csv = path_2_csvs + i
        map_coords_set = []
        # Extract points within 15min
        # Remove points in water
        with open(in_csv) as f:
            for x in csv.reader(f):
                # Lat, Lng, Mins, Mode
                if float(x[2]) <= duration:
                    if not point_in_water(float(x[0]), float(x[1])):
                        map_coords_set.append([x[0], x[1], x[3]])
                    else:
                        print("Skipping point in water at: %.5f, %.5f" % (float(x[0]), float(x[1])))  
                        
        # C: CREATE FINAL MAPS (full and simple)
        html_isochrone(map_coords_set, radius_km, in_csv.replace(".csv", "_"), centre_coord)
        html_isochrone(map_coords_set, radius_km, in_csv.replace(".csv", "__black"), centre_coord, black=True)

Step 3 (Creating the Final Output:

Obviously when you have around 25 HTML files you don’t want to manually screen-shot and photo-shop each one so I put together this script, which would automatically open each HTML, take a screenshot, then piece the screenshots together.

Python Scripts:

import os
import ntpath
from selenium import webdriver
from selenium.webdriver.support.ui import WebDriverWait
from selenium.webdriver.support import expected_conditions as EC
import time
from PIL import Image

def crop_image(path_nm):
    """Crop the screenshot to 700 by 700"""
    img = Image.open(path_nm)
    img_crop = img.crop((600,200,1300,900))
    img_crop.save(path_nm.replace('.png','_crop.png'))    

def screenshot_page(imlist):
    """Take screenshot of html; ideally add js code to change title to 'Map Loaded'"""
    options = webdriver.ChromeOptions()
    options.add_argument("--start-maximized")
    browser = webdriver.Chrome('C:/chromedriver.exe',chrome_options=options)
    for html in imlist:
        output = path + os.path.splitext(ntpath.basename(html))[0] + '.png'
        if not os.path.isfile(output):
            print('Converting %s to %s' % (html, output))
            browser.get('file://' + html)
            wait = WebDriverWait(browser, 10)
            # Require this in the javascript (tells selenium when to take the screenshot)
            #element = wait.until(EC.title_is('Map Loaded'))
            # Can also work with if stupidly time.sleep
            time.sleep(10)
            browser.save_screenshot(output)
            # Crop to square
            crop_image(output)
        else:
            print('Already exists: %s' % output)
    browser.quit()
    
if __name__ == '__main__':
    time.sleep(15)
    path = '////map_col/'
    htmllist = [os.path.join(root, f) for root,
              subdirs,files in os.walk(path) for f in files if os.path.splitext(f)[1].lower() in ('*.html')]
    print('List contains %d files' % len(htmllist))
    screenshot_page(htmllist)
"""Combine images into Grid"""
import os
from PIL import Image
from PIL import ImageFont
from PIL import ImageDraw 
from PIL import ImageEnhance

impath = '.../map_col/'
imlist = [impath+i for i in os.listdir(impath) if i.endswith("_crop.png")]
imlabel = [
            "Amsterdam - 425 (3872m)",
            "Bangkok - 834 (4996m)",
            "Barcelona - 327 (4228m)",
            "Berlin - 971 (4754m)",
            "Boston - 562 (7640m)",
            "Brussels - 1037 (5943m)",
            "Copenhagen - 744 (5544m)",
            "Geneva - 928 (4893m)",
            "London - 854 (4190m)",
            "Madrid - 380 (7452m)",
            "Moscow - 500 (5396m)",
            "Munich - 882 (6265m)",
            "New Delhi - 347 (3545m)",
            "New York City - 570 (4481m)",
            "Paris - 565 (4191m)",
            "Prague - 988 (6817m)",
            "Rome - 247 (1667m)",
            "San Francisco - 1268 (6907m)",
            "Sao Paulo - 403 (4650m)",
            "Sydney - 387 (6397m)",
            "Taipei - 130 (2264m)",
            "Vienna - 197 (1787m)",
            "Washington - 274 (1895m)",
            "Zurich - 1049 (5552m)"
        ]
counter = 0
col = 4
row = 6
output_img = Image.new('RGB',(700*col,700*row))

for i in range(0,700*col,700):
    for j in range(0,700*row,700):
        # Open Image
        img = Image.open((imlist)[counter])
        # Enhance (for the B+W version)
        #contr = ImageEnhance.Contrast(img)
        #bright = ImageEnhance.Brightness(img)
        #img = contr.enhance(1)
        #img = bright.enhance(4)
        # Add label
        draw = ImageDraw.Draw(img)
        font = ImageFont.truetype("arial.ttf", 30)
        draw.text((200,650),imlabel[counter],(0,0,0),font=font)
        # Add to canvas
        output_img.paste(img, (i,j)) 
        counter += 1
        
output_img.save(impath + 'transport_isochrone_15min_col.png') 
print('Done')

Finally, I wanted to emphasize the shape and have the viewer focus on the isochrones and hence the simple style. If you prefer having everything included check out the below (I think my favourite version would involve adding a white city boundary to each one of the black maps)

World Cities - 15min Public Transport (Colour)