Previously I wrote about generating isochrones by methodically extracting data from the Google Distance Matrix API. As an alternative to the concave-hull approach I wanted to try and make something that would look similar to isochrones produced using ArcGIS. For those not very familiar with ESRI’s software, they look something like this (10 minutes drive-time):

arcgis_10min_ST47QD

My end-result using Python looks like this:

ilia_10min_ST47QD_v2

Here is another example of the Python-method (30 minutes isochrone around Brixton, South London):

se59ef

I incorporated this method of generating an isochrone in the data-visualisation below:

GLA Pollution (NO2) + Brixton 30min Drive

The data visualisation shows the road-network in London (with road colours going from green to red as the level NO2 pollution increases in that area) with a white, 30 minute drive-time isochrone around South London. An attempt to see if one can escape the smog.

As an aside:

You will notice the Google Isochrones are smaller than their ArcGIS equivalents – here is another example (mine above and ArcGIS below)

result

The reason for this is two-fold:

First, Google really does give a slower travel-time for short distances (those less than 100 minutes), however after 100 minutes ArcGIS seems to get faster. Below is an in-house analysis of some sample routes with a line of best-fit.

gog

Second, imperfections in the method (filling a space with many points and routing to those) mean that sometimes a point lands on the other side of the road. This results in Google routing us to the nearest roundabout and asking the car to reverse.

Step-by-Step Method (with Python Script at the end)

Here is my attempt at producing isochrones that look similar to the ones produced by ArcGIS. The general idea is to take all the points that are reached in X minutes; fatten them up and merge them into a big polygon to show the rough area. Then, extract all the roads within that area to create the isochrone.

Step A: Generate the data by following the exact same routine described here to generate the data for the isochrone

Step B: Find yourself a road-network shapefile (I used the OS Open Roads network) and verify the coordinate system matches your data (the OS Open Roads data was actually in OS British National Grid 1936 and hence you can see in the code remnants of a conversion using pyproj, however I ended up converting the whole shape-file into latitude and longitude and thus commented those bits out later).

01_iso_roads_and_points

Now that we have the isochrone data and the road-network we can begin!

Step 1: We create a template polygon shape-file by buffering our points so that no gaps are left in the hexagonal-grid spacing (i.e. we want to buffer by the radius size we used to generate the isochrone and add on 10%). We then form a union of all the polygons (buffered points) and dissolve them to remove any holes in-between. We end up with a multi-polygon shape. (Technically we would want to buffer in a coordinate system where the spacing between the x and y axis is constant e.g. OSGB36, otherwise we will not get circles but elongated ellipses – which you can actually see below. However since we are creating just a template to cut with, for simplicity I ignored this)

02_points_expanded

Step 2: We clip our road-network (lines) with the template polygon and produce a shape-file (lines) which contains only the roads that fell within our template. I tried many different commands and packages to do this actually (gIntersection, gClip, gDifference, rTree, pyclipper, SQL commands in ogr, clipper and ogr2ogr) and none were as good as the Intersection command in ArcMap. I would love to hear of a better way to do this! I also couldn’t get rTree to work with my version of Python (3.4) which is a same as the examples I found here and here looked interesting. Hence, I settled on using ogrinfo to create a index and then ogr2ogr with the -spat parameter to clip only within a bounding box.

03_intersect_with_road

Step 3: Now we buffer our multi-line shape to convert it into a multi-polygon shape. You will notice that we get gaps – I’m not sure it is entirely obvious that they should be filled (imagine a park or a lake surrounded by roads), however I wanted to replicate ArcGIS isochrones and from what I have seen, they fill them in. The only issue is that gaps can be (1) true gaps or (2) fake gaps erroneously caused by imperfections in the algorithm, and only the latter should be closed.

04_road_line_expanded

Finally we create a union of all the polygons and dissolve them:

ilia_10min_ST47QD_v2

import fiona
import subprocess
import time
from shapely.geometry import Point, mapping, shape, MultiPolygon, Polygon
from shapely.ops import unary_union
from fiona.crs import from_epsg

# Preliminary Globals:
os.environ["GDAL_DATA"] = 'C:/Anaconda3/Lib/site-packages/osgeo/data/gdal'
my_roads_shape = 'K:/compdata/OS_Shape_Files/os_open_roads/trim/ROADLINK_Project.shp'
produce_poly_shape = NAME + '_rough.shp'
interim_poly_shape = NAME + "_interim.shp"
final_poly_shape = NAME + "_final.shp"

# NAME
# RADIUS_KM
# duration
# produce_csv_file
# <- were defined in the previous script used to generate the data

# Projections:
# from pyproj import Proj, transform
#wgs84=Proj("+init=EPSG:4326")  # LatLon with WGS84 datum used by GPS units and Google Earth
#osgb36=Proj("+init=EPSG:27700")  # UK Ordnance Survey, 1936 datum

#########
# STEP 1
#########
# Create shape-file of buffered points
# Extract only points from CSV that are within X minutes
pts = []
with open(produce_csv_file) as f:
    for x in csv.reader(f):
        # Lat, Lng, Minutes and keep those points within
        if float(x[2]) <= float(duration):
            # Shapely:
            pts.append(Point(
                # Transform from WGS84 to OSGB36 (if the roads shape-file is in OSGB36)
                # Here both are in WGS84 so we don't transform
                # transform(wgs84, osgb36, float(x[1]), float(x[0]))
                float(x[1]), float(x[0])
            ))

# Buffer points (divide spacing by 100) + 10% -> polygons
buffer = [point.buffer(1.1*float(RADIUS_KM)/100) for point in pts]
# Union -> multi-polygons
merged = unary_union(buffer)
# Remove holes (dissolve)
exterior_polys = []
for poly_l in merged:
    # Take only exterior coordinates
    exterior_polys.append(Polygon(poly_l.exterior))
merged = MultiPolygon(exterior_polys)

# Save the interim shape-file (a rough template to cut with)
schema = {'geometry': 'Polygon', 'properties': {'id': 'int'}}
with fiona.open(produce_poly_shape, "w", driver='ESRI Shapefile', crs=from_epsg(4326), schema=schema) as output:
    output.write({'geometry': mapping(merged), 'properties': {'id': 123}})

#########
# STEP 2
#########
# Find the bounds of our isochrone
# So that we first clip a rectangle and then we clip by the polygon (faster)
bounds = merged.bounds
x_min = '%.5f' % bounds[0]
y_min = '%.5f' % bounds[1]
x_max = '%.5f' % bounds[2]
y_max = '%.5f' % bounds[3]
print(x_min, x_max, y_min, y_max)

# Clip using GDAL
# Clip of lines using polygon -> lines
# http://www.gisinternals.com/release.php
print("Generation Time to Rough Shapefile: %.0f seconds" % (time.clock() - start))
print("Beginning Clip")
start = time.clock()

# Make sure shape-file has an index
# The R-Tree spatial-index looks great however haven't found a workin version for
# Python 3.4 (https://pypi.python.org/pypi/Rtree/)
# So I have used OGRINFO to create an index
#ogrinfo "../ROADLINK_Project.shp" -sql "CREATE SPATIAL INDEX ON ROADLINK_Project"
subprocess.call(["C:\Program Files\GDAL\ogr2ogr",
                 "-f",
                 "ESRI Shapefile",
                 "-spat",
                 x_min, y_min, x_max, y_max,
                 "-clipsrc",
                 produce_poly_shape,
                 interim_poly_shape,
                 my_roads_shape])

print('Clipped in %.0f seconds' % (time.clock()-start))

#########
# STEP 3
#########
road_lines = ([shape(lin['geometry']) for lin in fiona.open(interim_poly_shape)])
buffer = []
# Buffer by constant amount (0.001 looks most like arcGIS)
for road in road_lines:
    buffer.append(road.buffer(float(0.1/100)))
# Union
merged = unary_union(buffer)
# Remove holes
exterior_polys = []
for poly_l in merged:
    exterior_polys.append(Polygon(poly_l.exterior))
merged = MultiPolygon(exterior_polys)

# Save the final shape-file
schema = {'geometry': 'Polygon','properties': {'id': 'int'}}
with fiona.open(final_poly_shape, "w", driver='ESRI Shapefile', crs=from_epsg(4326), schema=schema) as output:
    output.write({'geometry': mapping(merged), 'properties': {'id': 123}})

print('Complete! Total time-taken: %.0f' % (time.clock() - overall_start))