Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
Finding the Furthest Place in Boston from a Dunks
#!/usr/bin/env python3
import sys
import json
import requests
import geopandas as gpd
import shapely
from shapely import wkt
from shapely.geometry import MultiPoint, Point, box
from shapely.ops import voronoi_diagram
import matplotlib.pyplot as plt
from functools import reduce
import contextily as ctx
Other sources
* [City of Boston Boundary](
def get_dunks_locations(location, radius:float, maxMatches:int):
"""Gets all Dunkin Donut store locations and information
From JP Bulman:
with open('cache.json', 'r') as f:
newJson =
print('Got cache!')
return newJson
print('No cache')
url = ''
response =, allow_redirects=False, data={
'service': 'DSL',
'origin': location,
'radius': str(int(round(radius))),
'maxMatches': str(maxMatches)
content = response.json()
newJson = {"data": []}
for entry in content["data"]["storeAttributes"]:
entry["geoJson"] = json.loads(entry["geoJson"])
with open('cache.json', 'w') as f:
json.dump(newJson, f)
return newJson
def get_encompassing_circle(shape):
center = shape.centroid
center.reset_index(drop=True, inplace=True)
corners = gpd.GeoSeries([Point(shape.total_bounds[0], shape.total_bounds[1]),
Point(shape.total_bounds[2], shape.total_bounds[1]),
Point(shape.total_bounds[0], shape.total_bounds[3]),
Point(shape.total_bounds[2], shape.total_bounds[3])]) =
furthest_dist = center.distance(corners).max() / 5280
return center, furthest_dist
def find_furthest(dunks, area):
voron = gpd.GeoSeries(voronoi_diagram(MultiPoint([d for d in dunks]))) =
# TODO(brycew): this isn't _really_ correct, it fails on islands.
# If an island is entirely within a voronoi polygon, then there's no intersection with the
# exterior boundaries of either, so it's not considered a candidate, even though it is.
# The only other solution is to just consider the entire exterior of each city polygon
# as candidates, but then we'd have infinite (not really, but still a crapton)
# points on the resulting lines to check, which breaks the brute force distance stuff.
# I'm not writing another fucking distance function for this.
myis = [area_part.exterior.intersection(poly.exterior) for area_part in area.geometry[0] for poly in voron.item()]
intersections = [list(inters) for inters in myis if not inters.is_empty]
list_inters = reduce(lambda l, y: y + l, intersections, [])
#showable_inters = gpd.GeoSeries([i[0] for i in y])
# Get individual points from all the voronoi polygons
corner_candidates = set(reduce(lambda l, y: y.exterior.coords[:] + l, voron.item(), []))
border_candidates = set(reduce(lambda l, y: y.coords[:] + l, list_inters, []))
candidates = corner_candidates.union(border_candidates)
candidate_dunk = None
best_candidate = None
best_dist = -1 * float('Infinity')
# Find the voronoi point that's furthest away from all dunks
for p in candidates:
closest_dunk = None
min_dist = float('Infinity')
pp = Point(p[0], p[1])
if not area.contains(pp).any():
for d in dunks:
dist = d.distance(pp)
if dist < min_dist:
min_dist = dist
closest_dunk = d
if min_dist > best_dist:
candidate_corner = pp
best_dist = min_dist
candidate_dunk = closest_dunk
close = gpd.GeoSeries(candidate_corner)
close_dunk = gpd.GeoSeries(candidate_dunk) = =
return close, close_dunk
def draw_map(area, dunks, close, close_dunk, area_name:str):
# Change to Web Mercator, which is what most Web map tiles use
area = area.to_crs(epsg=3857)
dunks = dunks.to_crs(epsg=3857)
close = close.to_crs(epsg=3857)
close_dunk = close_dunk.to_crs(epsg=3857)
# voron = voron.to_crs(epsg=3857)
base = None
# base = voron.plot(ax=base if base else None, color='black')
base = area.boundary.plot(ax=base if base else None, color='orange', label=area_name)
base = dunks.plot(ax=base, color='red', marker='o', markersize=20, label='Dunks Stores')
base = close.plot(ax=base, color='blue', marker='o', markersize=60, label='Furthest Point from Dunks')
base = close_dunk.plot(ax=base, color='black', marker='o', label="It's closest dunks")
plt.legend(loc='upper left')
ctx.add_basemap(base) #, url=ctx.providers.Stamen.TonerLite)
def main():
# Some ideas of other files to try it on:,
# And choosing the name of a neighborhood
if len(sys.argv) > 1:
shapefile = sys.argv[1]
if len(sys.argv) > 3:
name = sys.argv[2]
display_name = sys.argv[3]
name = None
display_name = ''
# Downloaded from
#shapefile = 'zip://'
# Downloaded from the same location above, but manually edited to remove the Harbor islands
shapefile = 'boston_no_islands.geojson'
name = None
display_name = 'Boston City Limits'
area = gpd.read_file(shapefile)
if name is not None:
area = area[area['Name'] == name]
area.reset_index(inplace=True, drop=True)
area = area.to_crs(epsg=2249)
center, furthest_dist = get_encompassing_circle(area)
print(f'Furthest distance from corner to Boston: {furthest_dist}')
#import code
lat_long = center.to_crs('EPSG:4326')
all_locations = get_dunks_locations(f'{lat_long[0].y},{lat_long[0].x}', furthest_dist, 100000)
dunks = gpd.GeoSeries([Point(loc['geoJson']['coordinates'][1],
loc['geoJson']['coordinates'][0]) for loc in all_locations["data"]]) = 'EPSG:4326'
dunks = dunks.to_crs(
close, close_dunk = find_furthest(dunks, area)
draw_map(area, dunks, close, close_dunk, display_name)
if __name__ == '__main__':
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment