In [7]:
import sys
import json
from pathlib import Path
from math import pi

import requests
from shapely.geometry import shape, Point
import geopandas as gpd
import spectra

sys.path.append('.')

import gtfs_tools as gt


DATA_DIR = Path('../data')
OUT_DIR = Path('../output')

%matplotlib inline

Problem 2


In [8]:
# Played around with https://mapzen.com/mobility/explorer/ and copied API query

url = 'https://matrix.mapzen.com/isochrone'
query = {"locations":[{"lat":"-36.84466680273936","lon":"174.7675895690918"}],"costing":"pedestrian","denoise":0.3,"polygons":True,"generalize":50,"costing_options":{"pedestrian":{"use_ferry":0}},"contours":[{"time":60}]}
api_key = 'valhalla-Mc6zgDA'
params = {'json': json.dumps(query), 'api_key': api_key}
r = requests.get(url, params=params)
r.json()
json.dumps(r.json())


Out[8]:
'{"type": "FeatureCollection", "features": [{"type": "Feature", "properties": {"contour": 60, "fill-opacity": 0.33, "fill": "#bf4040"}, "geometry": {"type": "Polygon", "coordinates": [[[174.749405, -36.832985], [174.750839, -36.832985], [174.750854, -36.838047], [174.745544, -36.838226], [174.750916, -36.838337], [174.750992, -36.840069], [174.752762, -36.840122], [174.752838, -36.841888], [174.75293, -36.838314], [174.754684, -36.838261], [174.754745, -36.836514], [174.756516, -36.836418], [174.756561, -36.834717], [174.761795, -36.834713], [174.761856, -36.840057], [174.761963, -36.838295], [174.769058, -36.838295], [174.769089, -36.840057], [174.770874, -36.840099], [174.77095, -36.841873], [174.771027, -36.83831], [174.778091, -36.838303], [174.778091, -36.841766], [174.776321, -36.841866], [174.783508, -36.841934], [174.783554, -36.843678], [174.787125, -36.843746], [174.787201, -36.845509], [174.787323, -36.841949], [174.790726, -36.841949], [174.790787, -36.843685], [174.792542, -36.843754], [174.792603, -36.845497], [174.794342, -36.84557], [174.794403, -36.847305], [174.796143, -36.847385], [174.796204, -36.849121], [174.799728, -36.849216], [174.799866, -36.85096], [174.807205, -36.850967], [174.807373, -36.849342], [174.809082, -36.849079], [174.807724, -36.848457], [174.807724, -36.846085], [174.81073, -36.845699], [174.812057, -36.845943], [174.812531, -36.849258], [174.816452, -36.84938], [174.81752, -36.848644], [174.819763, -36.848568], [174.820663, -36.852695], [174.81781, -36.856174], [174.817688, -36.857841], [174.814285, -36.858078], [174.814835, -36.85944], [174.813995, -36.861393], [174.809921, -36.860741], [174.809021, -36.856209], [174.805542, -36.856075], [174.805313, -36.854408], [174.805115, -36.856125], [174.803406, -36.856319], [174.803314, -36.857937], [174.801682, -36.858028], [174.800064, -36.857937], [174.799911, -36.856258], [174.798233, -36.856155], [174.798111, -36.854458], [174.796219, -36.854458], [174.796097, -36.857964], [174.794449, -36.858044], [174.794281, -36.859779], [174.792587, -36.859882], [174.792572, -36.861794], [174.796356, -36.861847], [174.796707, -36.860382], [174.798065, -36.860157], [174.804642, -36.860592], [174.806274, -36.86174], [174.805283, -36.863564], [174.806915, -36.86536], [174.800735, -36.869858], [174.79985, -36.871441], [174.798065, -36.871971], [174.794449, -36.875916], [174.792892, -36.874161], [174.791595, -36.874405], [174.79129, -36.87849], [174.788284, -36.879093], [174.785355, -36.881603], [174.783417, -36.881813], [174.781784, -36.881058], [174.778198, -36.883465], [174.776352, -36.883736], [174.774551, -36.881401], [174.773087, -36.883446], [174.767319, -36.883499], [174.765427, -36.885258], [174.761887, -36.886536], [174.75827, -36.883648], [174.756454, -36.884098], [174.754654, -36.8829], [174.749222, -36.88588], [174.743256, -36.88036], [174.740173, -36.880112], [174.738388, -36.877811], [174.733887, -36.875278], [174.729324, -36.870155], [174.725571, -36.869129], [174.726624, -36.865364], [174.722824, -36.863548], [174.722153, -36.859863], [174.720825, -36.858124], [174.720947, -36.85564], [174.718491, -36.853844], [174.717682, -36.850887], [174.719193, -36.849808], [174.72403, -36.849201], [174.722488, -36.848675], [174.722488, -36.845871], [174.724045, -36.845592], [174.724304, -36.844051], [174.725815, -36.843773], [174.726059, -36.842197], [174.729416, -36.841923], [174.729553, -36.840263], [174.731216, -36.840111], [174.731384, -36.838478], [174.734497, -36.838478], [174.734756, -36.84013], [174.734924, -36.838398], [174.738419, -36.838276], [174.738541, -36.83477], [174.740219, -36.83466], [174.740326, -36.832947], [174.749405, -36.832985]]]}}]}'

Problem 3


In [9]:
path = DATA_DIR/'auckland_gtfs_20161017.zip'
feed = gt.read_gtfs(path)

In [10]:
s = feed['stops']
cond = s['stop_name'].str.contains('train', case=False)
cond &= s['parent_station'].isnull()
cond &= s['location_type'] == 0
train_stations = s[cond].copy()
train_stations


Out[10]:
stop_lat zone_id stop_lon stop_id parent_station stop_desc stop_name location_type stop_code
738 -37.01346 NaN 174.87467 0099 NaN NaN Homai Train Station 0 99
739 -37.02327 NaN 174.89617 0098 NaN NaN Manurewa Train Station 0 98
740 -37.06429 NaN 174.94611 0097 NaN NaN Papakura Train Station 0 97
903 -37.20331 NaN 174.91016 0134 NaN NaN Pukekohe Train Station 0 134
905 -36.89778 NaN 174.84967 0130 NaN NaN Panmure Train Station 0 130
906 -36.84429 NaN 174.76848 0133 NaN NaN Britomart Train Station 0 133
1263 -36.89848 NaN 174.80811 0112 NaN NaN Ellerslie Train Station 0 112
1264 -36.88966 NaN 174.79742 0113 NaN NaN Greenlane Train Station 0 113
1265 -36.93454 NaN 174.83130 0111 NaN NaN Westfield Train Station 0 111
1266 -36.86243 NaN 174.80950 0116 NaN NaN Orakei Train Station 0 116
1267 -36.86632 NaN 174.82075 0117 NaN NaN Meadowbank Train Station 0 117
1268 -36.88138 NaN 174.78542 0114 NaN NaN Remuera Train Station 0 114
1269 -36.86972 NaN 174.77883 0115 NaN NaN Newmarket Train Station 0 115
1270 -36.86796 NaN 174.75899 0118 NaN NaN Mt Eden Train Station 0 118
1271 -36.87768 NaN 174.72045 0119 NaN NaN Baldwin Ave Train Station 0 119
2370 -36.92013 NaN 174.80137 0606 NaN NaN Te Papapa Train Station 0 606
2372 -36.92539 NaN 174.78680 0605 NaN NaN Onehunga Train Station 0 605
2912 -36.86561 NaN 174.76984 0277 NaN NaN Grafton Train Station 0 277
6070 -36.86618 NaN 174.57630 0127 NaN NaN Swanson Train Station 0 127
6071 -36.87347 NaN 174.62083 0126 NaN NaN Sturges Rd Train Station 0 126
6072 -36.88101 NaN 174.63098 0125 NaN NaN Henderson Train Station 0 125
6073 -36.91028 NaN 174.65317 0124 NaN NaN Glen Eden Train Station 0 124
6074 -36.91067 NaN 174.66707 0123 NaN NaN Fruitvale Rd Train Station 0 123
6075 -36.87251 NaN 174.74453 0122 NaN NaN Kingsland Train Station 0 122
6076 -36.89687 NaN 174.63201 0121 NaN NaN Sunnyvale Train Station 0 121
6077 -36.88478 NaN 174.71409 0120 NaN NaN Mt Albert Train Station 0 120
6078 -36.90939 NaN 174.68405 0129 NaN NaN New Lynn Train Station 0 129
6079 -36.86794 NaN 174.60334 0128 NaN NaN Ranui Train Station 0 128
6451 -36.96304 NaN 174.83897 0109 NaN NaN Middlemore Train Station 0 109
6452 -36.98976 NaN 174.85608 0108 NaN NaN Puhinui Train Station 0 108
6453 -36.89714 NaN 174.69911 0105 NaN NaN Avondale Train Station 0 105
6454 -36.87499 NaN 174.73521 0104 NaN NaN Morningside Train Station 0 104
6455 -37.03118 NaN 174.90609 0107 NaN NaN Te Mahia Train Station 0 107
6457 -36.94669 NaN 174.83321 0101 NaN NaN Otahuhu Train Station 0 101
6459 -36.87880 NaN 174.85412 0103 NaN NaN Glen Innes Train Station 0 103
6460 -36.91009 NaN 174.81570 0102 NaN NaN Penrose Train Station 0 102
8172 -36.91464 NaN 174.84260 0244 NaN NaN Sylvia Park Train Station 0 244
8910 -37.04112 NaN 174.91945 0106 NaN NaN Takanini Train Station 0 106
8918 -36.97767 NaN 174.84925 0100 NaN NaN Papatoetoe Train Station 0 100
9450 -36.99388 NaN 174.87740 9218 NaN NaN Manukau Train Station 0 9218

Problem 4


In [11]:
def get_walking_catchment(lon, lat, distance=1, polygons=True):
    """
    Issue a GET request to the Mapzen Isochrone API to get an approximation
    of the distance ``distance`` kilometer walking catchment centered at the 
    given WGS84 longitude and latitude.
    
    Return the result as a (decoded) GeoJSON Feature, which will contain 
    the catchment as a polygon if ``polygons``; otherwise the catchment boundary
    will be represented by a linestring.
    """
    url = 'https://matrix.mapzen.com/isochrone'
    api_key = 'valhalla-Mc6zgDA'
    query = {
        'locations':[{'lon': lon, 'lat':lat}],
        'polygons': polygons,
        'costing': 'pedestrian',
        'costing_options': {'pedestrian': {'walking_speed': 5.0, 'use_ferry': 0}},
        'contours': [{'time': 12*distance}],
        'denoise': 0.3,  
        'generalize': 50,
    }
    params = {'api_key': api_key, 'json': json.dumps(query) }
    r = requests.get(url, params=params)
    collection = r.json()
    return collection['features'][0]

lon = 174.76848
lat = -36.84429
json.dumps(get_walking_catchment(lon, lat, distance=2))


Out[11]:
'{"type": "Feature", "properties": {"contour": 24, "fill-opacity": 0.33, "fill": "#bf4040"}, "geometry": {"type": "Polygon", "coordinates": [[[174.776031, -36.837906], [174.777237, -36.837906], [174.777405, -36.841217], [174.779037, -36.841526], [174.779358, -36.843197], [174.779984, -36.841843], [174.781158, -36.84153], [174.782349, -36.841843], [174.782974, -36.84333], [174.785706, -36.84391], [174.786591, -36.845188], [174.788452, -36.844772], [174.789459, -36.845566], [174.790451, -36.84845], [174.788193, -36.850262], [174.784775, -36.852074], [174.781174, -36.856636], [174.775726, -36.856514], [174.773926, -36.859104], [174.77211, -36.858376], [174.763077, -36.860146], [174.760757, -36.859802], [174.75795, -36.857494], [174.755753, -36.857494], [174.754517, -36.85207], [174.749756, -36.850262], [174.748398, -36.842823], [174.752426, -36.841442], [174.754044, -36.838615], [174.75946, -36.839844], [174.761337, -36.839481], [174.76149, -36.837841], [174.76828, -36.837814], [174.768402, -36.8395], [174.770142, -36.839581], [174.770309, -36.841297], [174.770752, -36.839855], [174.772217, -36.839516], [174.772446, -36.837925], [174.776031, -36.837906]]]}}'

In [12]:
w = get_walking_catchment(lon, lat)
shape(w['geometry'])


Out[12]:

In [13]:
# Collect all walking catchments in a geodataframe.

WGS84_CRS = {'init': 'epsg:4326'}
NZTM_CRS = {'no_defs': True, 'init': 'epsg:2193'}

def my_apply(row):
    lon, lat = row
    w = get_walking_catchment(lon, lat)
    return shape(w['geometry'])
    
f = train_stations.copy()
f['geometry'] = f[['stop_lon', 'stop_lat']].apply(my_apply, axis=1)
f = gpd.GeoDataFrame(f[['stop_name', 'stop_lon', 'stop_lat', 'geometry']], crs=WGS84_CRS)
walking_catchments = f.copy()
f.plot()


Out[13]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f2fd32a6470>

In [14]:
# Collect all flying catchments in a geodataframe.
    
f = train_stations.copy()
f['geometry'] = f[['stop_lon', 'stop_lat']].apply(lambda row: Point(*row), axis=1)
f = gpd.GeoDataFrame(f[['stop_name', 'stop_lon', 'stop_lat', 'geometry']], crs=WGS84_CRS)
f = f.to_crs(NZTM_CRS)
f['geometry'] = f['geometry'].buffer(1000)
f = f.to_crs(WGS84_CRS)
flying_catchments = f.copy()
f.plot()


Out[14]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f2fd657aa90>

In [16]:
# Collect train stations into geodataframe

f = train_stations.copy()
f['geometry'] = f[['stop_lon', 'stop_lat']].apply(lambda row: Point(*row), axis=1)
f = gpd.GeoDataFrame(f[['stop_name', 'stop_lon', 'stop_lat', 'geometry']], crs=WGS84_CRS)
# Buffer some to form polygons
f['geometry'] = f['geometry'].to_crs(NZTM_CRS).buffer(20).to_crs(WGS84_CRS)
train_stations = f.copy()
f.plot()


Out[16]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f2fceb2d2e8>

Problem 5


In [17]:
walking_catchments['area_ratio'] = walking_catchments.to_crs(NZTM_CRS)['geometry'].area/(pi*1000**2) 
walking_catchments


Out[17]:
stop_name stop_lon stop_lat geometry area_ratio
738 Homai Train Station 174.87467 -37.01346 POLYGON ((174.875229 -37.007957, 174.877594 -3... 0.522511
739 Manurewa Train Station 174.89617 -37.02327 POLYGON ((174.898285 -37.016212, 174.900833 -3... 0.682678
740 Papakura Train Station 174.94611 -37.06429 POLYGON ((174.946564 -37.057079, 174.947922 -3... 0.622172
903 Pukekohe Train Station 174.91016 -37.20331 POLYGON ((174.911987 -37.194645, 174.916809 -3... 0.649253
905 Panmure Train Station 174.84967 -36.89778 POLYGON ((174.855591 -36.890594, 174.855865 -3... 0.611970
906 Britomart Train Station 174.76848 -36.84429 POLYGON ((174.767593 -36.838326, 174.768326 -3... 0.592793
1263 Ellerslie Train Station 174.80811 -36.89848 POLYGON ((174.811844 -36.889858, 174.813614 -3... 0.619149
1264 Greenlane Train Station 174.79742 -36.88966 POLYGON ((174.802963 -36.882847, 174.804749 -3... 0.602659
1265 Westfield Train Station 174.83130 -36.93454 POLYGON ((174.835449 -36.927303, 174.839493 -3... 0.368494
1266 Orakei Train Station 174.80950 -36.86243 POLYGON ((174.811737 -36.857128, 174.815018 -3... 0.305153
1267 Meadowbank Train Station 174.82075 -36.86632 POLYGON ((174.815979 -36.862606, 174.817047 -3... 0.314657
1268 Remuera Train Station 174.78542 -36.88138 POLYGON ((174.789185 -36.87455, 174.793823 -36... 0.663345
1269 Newmarket Train Station 174.77883 -36.86972 POLYGON ((174.782593 -36.861095, 174.782227 -3... 0.615930
1270 Mt Eden Train Station 174.75899 -36.86796 POLYGON ((174.761566 -36.860523, 174.769104 -3... 0.678439
1271 Baldwin Ave Train Station 174.72045 -36.87768 POLYGON ((174.726013 -36.870857, 174.728943 -3... 0.614967
2370 Te Papapa Train Station 174.80137 -36.92013 POLYGON ((174.806976 -36.913258, 174.806595 -3... 0.565115
2372 Onehunga Train Station 174.78680 -36.92539 POLYGON ((174.785553 -36.916321, 174.788635 -3... 0.621842
2912 Grafton Train Station 174.76984 -36.86561 POLYGON ((174.771759 -36.857006, 174.775833 -3... 0.617033
6070 Swanson Train Station 174.57630 -36.86618 POLYGON ((174.584518 -36.862106, 174.585449 -3... 0.320911
6071 Sturges Rd Train Station 174.62083 -36.87347 POLYGON ((174.621033 -36.868412, 174.626358 -3... 0.458427
6072 Henderson Train Station 174.63098 -36.88101 POLYGON ((174.627762 -36.872128, 174.631088 -3... 0.645281
6073 Glen Eden Train Station 174.65317 -36.91028 POLYGON ((174.65538 -36.901363, 174.656143 -36... 0.700176
6074 Fruitvale Rd Train Station 174.66707 -36.91067 POLYGON ((174.671326 -36.903339, 174.674377 -3... 0.602174
6075 Kingsland Train Station 174.74453 -36.87251 POLYGON ((174.750519 -36.865269, 174.754456 -3... 0.525485
6076 Sunnyvale Train Station 174.63201 -36.89687 POLYGON ((174.633896 -36.887314, 174.632889 -3... 0.492508
6077 Mt Albert Train Station 174.71409 -36.88478 POLYGON ((174.710632 -36.876144, 174.713211 -3... 0.628699
6078 New Lynn Train Station 174.68405 -36.90939 POLYGON ((174.690063 -36.902092, 174.692444 -3... 0.721825
6079 Ranui Train Station 174.60334 -36.86794 POLYGON ((174.601898 -36.859074, 174.606567 -3... 0.674085
6451 Middlemore Train Station 174.83897 -36.96304 POLYGON ((174.835587 -36.954288, 174.838882 -3... 0.442938
6452 Puhinui Train Station 174.85608 -36.98976 POLYGON ((174.85466 -36.980831, 174.866196 -36... 0.550285
6453 Avondale Train Station 174.69911 -36.89714 POLYGON ((174.701065 -36.888493, 174.703934 -3... 0.642214
6454 Morningside Train Station 174.73521 -36.87499 POLYGON ((174.739304 -36.867844, 174.743271 -3... 0.637240
6455 Te Mahia Train Station 174.90609 -37.03118 POLYGON ((174.90979 -37.024353, 174.909515 -37... 0.450125
6457 Otahuhu Train Station 174.83321 -36.94669 POLYGON ((174.835114 -36.938076, 174.835815 -3... 0.439327
6459 Glen Innes Train Station 174.85412 -36.87880 POLYGON ((174.852554 -36.870079, 174.854233 -3... 0.650843
6460 Penrose Train Station 174.81570 -36.91009 POLYGON ((174.810791 -36.902885, 174.811996 -3... 0.542044
8172 Sylvia Park Train Station 174.84260 -36.91464 POLYGON ((174.843033 -36.907524, 174.844498 -3... 0.426502
8910 Takanini Train Station 174.91945 -37.04112 POLYGON ((174.92366 -37.033768, 174.925552 -37... 0.621110
8918 Papatoetoe Train Station 174.84925 -36.97767 POLYGON ((174.853806 -36.969967, 174.854141 -3... 0.616410
9450 Manukau Train Station 174.87740 -36.99388 POLYGON ((174.883148 -36.98682, 174.887894 -36... 0.627952

In [18]:
walking_catchments.plot(column='area_ratio', cmap='OrRd')


Out[18]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f2fd0087a90>

Problem 6


In [19]:
# Color-code by area ratio using Spectra

w = walking_catchments.copy()
cuts = [0, 0.25, 0.5, 0.75, 1]
colors = reversed(['#d7191c','#fdae61','#ffffbf','#abdda4','#2b83ba'])
scale = spectra.scale(colors).colorspace('lch').domain(cuts)
w['fill'] = w['area_ratio'].map(lambda x: scale(x).hexcode)
w['fill-opacity'] = 0.8

f = flying_catchments.copy()
f['fill'] = 'gray'
f['fill-opacity'] = 0.3

t = train_stations.copy()

# Export to GeoJSON and plot in geojsonio for a closer look  
path = OUT_DIR/'walking_catchments.geojson'
geo = w.to_json()
with path.open('w') as tgt:
    tgt.write(geo)
    
path = OUT_DIR/'flying_catchments.geojson'
geo = f.to_json()
with path.open('w') as tgt:
    tgt.write(geo)

path = OUT_DIR/'train_stations.geojson'
geo = t.to_json()
with path.open('w') as tgt:
    tgt.write(geo)

In [ ]: