starbucks.py

#!/usr/bin/env python
import sys
import os
import json
import random

MIN_DISTANCE = 0.0
MAX_DISTANCE = 1000000.0
MAX_SAMPLE = 1000000

# approximate
# https://stackoverflow.com/questions/4913349/
from math import radians, cos, sin, asin, sqrt
def haversine(lon1, lat1, lon2, lat2):
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles
    return c * r

def haversine_feature(o, d):
    return haversine(
        o['geometry']['coordinates'][0],
        o['geometry']['coordinates'][1],
        d['geometry']['coordinates'][0],
        d['geometry']['coordinates'][1]
    )

def generate_od(bbox, features, sample1=None, sample2=None, min_distance=MIN_DISTANCE, max_distance=MAX_DISTANCE):
    # simple bbox test
    left, bottom, right, top = bbox
    features = [
        i for i in features 
        if left <= i['geometry']['coordinates'][0] <= right and bottom <= i['geometry']['coordinates'][1] <= top
    ]

    # find and sample od pairs
    pairs = []
    for o in features:
        d = sorted(features, key=lambda d:haversine_feature(o,d) )[1]
        print "%s -> %s"%(o, d)
        pairs.append([o,d])
    return pairs

def load_starbucks(filename):
    features = []
    with open(filename, 'rb') as f:
        data = json.load(f)
        for i in data:
            feature = {
                "type":"Feature", 
                "properties":{
                    "city":i["city"], 
                    "name":i["name"], 
                    "country":i["country"],
                    "store_id":i["store_id"]
                },
                "geometry": {
                    "type":"LineString",
                    "coordinates": [i["longitude"], i["latitude"]]
                }
            }
            features.append(feature)
    return {"type": "FeatureCollection", "features": features}

if __name__ == "__main__":
    # load data
    # https://github.com/mmcloughlin/starbucks/blob/master/locations.json
    data = load_starbucks('starbucks.json')

    bbox = '-122.532063,37.716418,-122.355080,37.811140'
    # bbox = '-180,-90,180,90'
    bbox = map(float, bbox.split(','))

    # generate pairs
    od = generate_od(bbox, data['features'], sample1=0, sample2=0)

    # write output
    features = []
    for o,d in od:
        features.append({
            "type": "Feature",
            "properties": {
                "name": "%s (%s) -> %s (%s)"%(
                    o['properties']['name'], 
                    o['properties']['store_id'], 
                    d['properties']['name'], 
                    d['properties']['store_id']
                ),
                "from_store_id": o['properties']['store_id'],
                "to_store_id": d['properties']['store_id'],
                "distance": haversine_feature(o, d),
                "stroke": "#000ad3",
                "stroke-width": "2.0",
                "stroke-opacity": 1,
            },
            "geometry": {
                "type": "LineString",
                "coordinates": [
                    o['geometry']['coordinates'],
                    d['geometry']['coordinates'],
                ]
            }
        })
    features = sorted(features, key=lambda x:x['properties']['distance'])
    fc = {
        "type": "FeatureCollection",
        "features": features
    }
    with open('test.geojson', 'wb') as f:
        json.dump(fc, f)