Finding the Furthest Place in Boston from a Dunks

donuts.py

#!/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](https://data.boston.gov/dataset/city-of-boston-boundary)
* https://github.com/Toblerity/Shapely/blob/master/shapely/ops.py
* https://www.dunkindonuts.com/en/locations?location=02125
* https://geopandas.org/docs.html
"""


def get_dunks_locations(location, radius:float, maxMatches:int):
  """Gets all Dunkin Donut store locations and information
  From JP Bulman: https://www.kaggle.com/jpbulman/usa-dunkin-donuts-stores?select=dunkin.py
  """
  try:
    with open('cache.json', 'r') as f:
      newJson = json.read(f) 
      print('Got cache!')
      return newJson
  except:
    print('No cache')
    pass

  url = 'https://www.dunkindonuts.com/bin/servlet/dsl'
  response = requests.post(url, 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"])
    newJson["data"].append(entry)

  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])])
  corners.crs = shape.crs
  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])))
  voron.crs = dunks.crs

  # 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 lot) 
  # points on the resulting lines to check, which breaks the brute force distance stuff.
  # I'm not writing another 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():
      continue
    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)
  close.crs = dunks.crs
  close_dunk.crs = dunks.crs
  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)

  plt.show()


def main():
  # Some ideas of other files to try it on: https://data.boston.gov/dataset/boston-neighborhoods, 
  # 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]
    else:
      name = None
      display_name = ''
  else:
    # Downloaded from https://data.boston.gov/dataset/city-of-boston-boundary
    #shapefile = 'zip://City_of_Boston_Boundary.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)
  print(area.crs)
  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
  #code.interact(local=locals())
  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"]])
  dunks.crs = 'EPSG:4326'
  dunks = dunks.to_crs(area.crs)
  close, close_dunk = find_furthest(dunks, area)
  draw_map(area, dunks, close, close_dunk, display_name) 

if __name__ == '__main__':
  main()