Output has been removed from this notebook to reduce file sizes in the repo



In [ ]:

    
%load_ext autoreload
%autoreload 2
import pandas as pd
import numpy as np
import pandana as pdna
import geopandas.io.osm as osm
%matplotlib inline

Download OpenStreetMap restaurants for a good part of the Bay Area

Note: used http://boundingbox.klokantech.com/ to get the bounding box



In [ ]:

    
gdf = osm.query_osm('node', 
                    bbox=[-122.8662,37.1373,-121.4798,38.2158],
                    tags='amenity=restaurant')



In [ ]:

    
gdf = gdf[gdf.type == 'Point'].to_crs(epsg=3740)
print gdf.geometry.head(3)
print len(gdf)



In [ ]:

    
x, y = zip(*[(p.x, p.y) for (i, p) 
             in gdf.geometry.iteritems()])
x = pd.Series(x)
y = pd.Series(y)

Get OpenStreetMap networks for Bay Area that I had previously - someday soon we'll have direct OSM import



In [ ]:

    
store = pd.HDFStore('data/osm_bayarea.h5', "r")
nodes = store.nodes
edges = store.edges
print nodes.head(3)
print edges.head(3)

Initialize and preprocess the network



In [ ]:

    
net=pdna.Network(nodes.x, 
                       nodes.y, 
                       edges["from"], 
                       edges.to, 
                       edges[["weight"]])
net.precompute(3000)

Nearest point-of-interest queries



In [ ]:

    
net.init_pois(num_categories=1, max_dist=2000, max_pois=10)



In [ ]:

    
net.set_pois("restaurants", x, y)



In [ ]:

    
a = net.nearest_pois(2000, "restaurants", num_pois=10)
print a.head(1)



In [ ]:

    
from shapely.geometry import Point
from fiona.crs import from_epsg
import geopandas as gpd
bbox=[-122.539365,37.693047,-122.347698,37.816069]
bbox = gpd.GeoSeries([Point(bbox[0], bbox[1]),
                      Point(bbox[2], bbox[3])], 
                     crs=from_epsg(4326))
bbox = bbox.to_crs(epsg=3740)
bbox = [bbox[0].x, bbox[0].y, bbox[1].x, bbox[1].y]

Here's a map of the distance to the nearest restaurant



In [ ]:

    
net.plot(a[1], bbox=bbox, scheme="diverging", 
         color="BrBG")

Here's a map of the distance to the 5th nearest restaurant



In [ ]:

    
net.plot(a[5], bbox=bbox, scheme="diverging", 
         color="BrBG")

Here's a map of the distance to the 10th nearest restaurant



In [ ]:

    
net.plot(a[10], bbox=bbox, scheme="diverging", 
         color="BrBG")

A similar workflow is used to do general network aggregations

Relate the x-ys to nodes



In [ ]:

    
node_ids = net.get_node_ids(x, y)

Assign the variable (in this case just location) to the network



In [ ]:

    
net.set(node_ids)

This is it - run the queries!



In [ ]:

    
%time s = net.aggregate(500, type="sum", decay="linear")
%time t = net.aggregate(1000, type="sum", decay="linear")
%time u = net.aggregate(2000, type="sum", decay="linear")
%time v = net.aggregate(3000, type="sum", decay="linear")
%time w = net.aggregate(3000, type="count", decay="flat")

Here's a map of access to restaurants with a 500m radius



In [ ]:

    
net.plot(s, bbox=bbox, scheme="diverging", 
         color="BrBG", log_scale=True)

Or 1000 meters



In [ ]:

    
net.plot(t, bbox=bbox, scheme="diverging", 
         color="BrBG", log_scale=True)

Or 2000 meters radius



In [ ]:

    
net.plot(u, bbox=bbox, scheme="diverging", 
         color="BrBG", log_scale=True)

Or 3000m radius



In [ ]:

    
net.plot(v, bbox=bbox, scheme="diverging", 
         color="BrBG", log_scale=True)

Or the whole Bay Area region - someone please help me with this visualization!



In [ ]:

    
net.plot(v, scheme="diverging", 
         color="BrBG", log_scale=True)



In [ ]: