{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Analysis of Gas Prices in the GTA\n", "We obtain gas station price data for the Greater Toronto Area with data from yellowpages.ca and gasbuddy.com and cluster them based on location and price per litre of regular gas. We also test the price effects of having another gas station or Starbucks location within a 100m radius of a gas station." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from bs4 import BeautifulSoup\n", "import urllib\n", "import requests\n", "import math\n", "import re\n", "import pandas as pd\n", "import numpy as np\n", "import datetime\n", "import sys\n", "import time\n", "import json\n", "import ast\n", "from geopy.geocoders import GoogleV3\n", "import seaborn as sns\n", "from time import sleep\n", "geolocator = GoogleV3()\n", "%matplotlib inline\n", "today = datetime.datetime.now()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Get the names and addressess of gas stations listed in Yellow Pages in the Toronto region\n", "\n", "def get_addresses (page_num):\n", " with open('gas_dict.json', 'r') as infile:\n", " gas_dict = json.load(infile)\n", " \n", " url = 'https://www.yellowpages.ca/search/si/'+str(page_num)+'/Gas%20Stations/Toronto+ON'\n", " try:\n", " html = urllib.request.urlopen(url).read()\n", " except:\n", " return page_num\n", " soup = BeautifulSoup(html, 'html.parser')\n", " \n", " content = soup.findAll('div',{'class':'listing__content__wrap'})\n", "\n", " for n in range(len(content)):\n", " try:\n", " name = content[n].findAll('a')[0].get('title').split('-')[0]\n", " except:\n", " name = '#NO NAME'\n", " \n", " try:\n", " addr_ele = content[n].findAll('span',{'class':'listing__address--full'})[0].findAll('span')\n", " addr = ','.join([addr_ele[i].string for i,e in enumerate(addr_ele)])\n", " except:\n", " addr = '#NO ADDRESS#'\n", " \n", " gas_dict[addr] = {'name':name}\n", " \n", " with open('gas_dict.json', 'w') as outfile:\n", " json.dump(gas_dict,outfile)\n", " \n", " return page_num + 1\n", "\n", "def address_main():\n", " # If file doesn't exist, create it\n", " try:\n", " with open('gas_dict.json', 'r') as infile:\n", " gas_dict = json.load(infile)\n", " except:\n", " with open('gas_dict.json', 'w') as outfile:\n", " gas_dict={}\n", " json.dump(gas_dict,outfile)\n", " \n", " page_num = 1\n", " while page_num < 22:\n", " print('Fetching page ',page_num)\n", " page_num = get_addresses(page_num)\n", " sleep(np.random.random(1)*2+2)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Get rid of incomplete or empty addresses\n", "def clean_addresses():\n", " with open('gas_dict.json', 'r') as infile:\n", " gas_dict = json.load(infile)\n", " \n", " gas_dict.pop('#NO ADDRESS#',None)\n", " gas_dict.pop('ON',None)\n", " gas_dict.pop('12001 Hwy 400,Maple,ON,L7B 1A8',None) # This is actually both redundant and the wrong city ...\n", " \n", " with open('gas_dict.json', 'w') as outfile:\n", " json.dump(gas_dict,outfile)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# for ontariogasprices.com scraping, we need a list of postal codes to search with\n", "# the site focuses on finding the lowest price in an area and does not prioritize towards specific addresses\n", "# but this works to our advantage since we can potentially find stations not included in the yellowpages data\n", "# gas types are\n", "# A: regular, B: mid, C: premium, D: diesel\n", "\n", "def get_area_prices(gas_type, postal_code):\n", " with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'r') as infile:\n", " gas_stations = json.load(infile)\n", " \n", " url = 'http://www.ontariogasprices.com/GasPriceSearch.aspx?fuel='+gas_type+'&qsrch='+postal_code[0]+'%20'+postal_code[1]\n", " print ('Fetching from ', url) \n", " html = requests.get(url).text\n", " soup = BeautifulSoup(html, 'html.parser') \n", "\n", " if len(soup.findAll('tbody')) == 0:\n", " pass\n", " else:\n", " data_rows = soup.findAll('tbody')[0].findAll('tr')\n", "\n", " for i in range(len(data_rows)):\n", " try:\n", " name = data_rows[i].findAll('td')[0].find('a').string\n", " st_addr = data_rows[i].findAll('td')[0].find('dd').string.split('&')[0].split('near')[0]\n", " city = data_rows[i].findAll('td')[2].find('a').string\n", " addr = st_addr + ', '+ city + ', Ontario'\n", " price = float(data_rows[i].find('a').string)\n", " if addr not in gas_stations:\n", " gas_stations[addr] = {'name':name}\n", " else:\n", " pass\n", " gas_stations[addr]['type_'+gas_type] = price\n", " except:\n", " # The \"other stations in the area\" row can just be ignored\n", " pass\n", " \n", " with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'w') as outfile:\n", " json.dump(gas_stations,outfile)\n", " \n", "def price_main():\n", " # first get postal codes\n", " with open('gas_dict.json', 'r') as infile:\n", " gas_dict = json.load(infile)\n", " \n", " try:\n", " with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'r') as infile:\n", " gas_stations = json.load(infile)\n", " except:\n", " with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'w') as outfile:\n", " gas_stations={}\n", " json.dump(gas_stations,outfile)\n", " \n", " postal_list=[]\n", " for k,v in gas_dict.items():\n", " postal_list.append(k.split(',')[-1].split())\n", " \n", " # If no postal code, just ignore for now...\n", " postal_list = [postal_list[i] for i,p in enumerate(postal_list) if len(p) == 2]\n", "\n", " break_point = 0\n", " while (break_point < len(postal_list)):\n", " postal_code = postal_list[break_point]\n", " try:\n", " for gas_type in list('ABCD'):\n", " get_area_prices(gas_type, postal_code)\n", " sleep(np.random.random(1)*2+0.5) \n", " break_point += 1\n", " except:\n", " # Try again\n", " print('Trying again ...')\n", " sleep(np.random.random(1)*2+2)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Get not only the coordinates but also the Google Maps address, which will likely clean up any messes left over\n", "# from the ontariogasprices.com splits\n", "\n", "def get_coord(addr):\n", " # splitting Toronto addresses into \"Toronto - [area]\" messes up google maps, so get rid of that when querying\n", " if 'Toronto' in addr.split(',')[1]:\n", " addr = addr.split(',')[0] + ', Toronto, '+addr.split(',')[2]\n", " print ('Getting' , addr)\n", " geoloc = geolocator.geocode(addr)\n", " sleep(np.random.random(1)*2+2)\n", " print ('Got ', geoloc, '\\n')\n", " return geoloc\n", " \n", "def coords_main():\n", " # Start at the beginning of gas_dict\n", " break_point = 0\n", "\n", " # open gas_stations, and get the list of all the addresses\n", " with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json', 'r') as infile:\n", " gas_dict = json.load(infile)\n", " addr_list = [k for k,v in gas_dict.items()]\n", " \n", " # do while there are still addresses to process\n", " while(break_point < len(addr_list)):\n", " # run collect_coords starting with the current breakpoint, and gets the\n", " # address that it's stuck on\n", " addr = addr_list[break_point]\n", " try: \n", " geoloc = get_coord(addr)\n", " if geoloc == None:\n", " # GoogleV3 can't get a hold of it. Need to look at it later. Most likely it's a highway. \n", " break_point += 1\n", " sleep(np.random.random(1)*2+2) \n", " else:\n", " gas_dict[addr]['address'] = geoloc.address\n", " gas_dict[addr]['longitude'] = geoloc.longitude\n", " gas_dict[addr]['latitude'] = geoloc.latitude\n", " # If it works, get to the next one. Once break_point gets to the length of the list, the conditions of\n", " # the loop will no longer be satisfied, and the loop ends.\n", " break_point +=1\n", " with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json','w') as outfile:\n", " json.dump(gas_dict,outfile) \n", " sleep(np.random.random(1)*2+2) \n", " except:\n", " # If it doesn't work, don't increase the break_point and try again\n", " print ('Could not get ', addr)\n", " sleep(np.random.random(1)*2+2)\n", "\n", " with open('gas_stations_'+datetime.datetime.now().strftime('%Y-%m-%d')+'.json','w') as outfile:\n", " json.dump(gas_dict,outfile) \n", " #return gas_dict" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def xy_distance(xlat,xlong,ylat,ylong):\n", " # The radius of the planet is approx. 6371.01 km\n", " xlat = math.radians(xlat)\n", " xlong = math.radians(xlong)\n", " ylat = math.radians(ylat)\n", " ylong = math.radians(ylong)\n", " dist = 6371.01 * math.acos(math.sin(xlat)*math.sin(ylat) + math.cos(xlat)*math.cos(ylat)*math.cos(xlong - ylong))\n", " return dist\n", "\n", "def xy_closest(coords, x):\n", " return (min(coords, key=lambda y: xy_distance(x[1],x[0],y[1],y[0])))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def enter_address(df, addr, max_dist):\n", " standard_addr = get_coord(addr)\n", " within_range = pd.DataFrame()\n", " for i in range(len(df)):\n", " dist = xy_distance(standard_addr.latitude,standard_addr.longitude,df.iloc[i]['latitude'],df.iloc[i]['longitude'])\n", " if dist <= max_dist:\n", " within_range = pd.concat((within_range,df.iloc[i]),axis=1)\n", "\n", " within_range = within_range.transpose().sort_values('type_A')\n", " \n", " # Hide the entered address in the dataframe too so we can show it on the map\n", " print ('Gas stations within range: ')\n", " print (within_range[['name','type_A']])\n", " folmap(within_range,15,custom=1,extralat=standard_addr.latitude,extralong=standard_addr.longitude)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def folmap (df,zoom=10,custom=None,extralat=None,extralong=None):\n", " k=7\n", " colorscale = branca.colormap.linear.YlGnBu.scale(0, 60)\n", " cluster_colors = ['red', 'green', 'blue', 'yellow', 'cadetblue', 'purple', 'black']\n", "\n", " if custom == 1:\n", " m = folium.Map(\n", " location=[extralat, extralong],\n", " zoom_start=zoom\n", " )\n", " else:\n", " # center the map on the median coordinates\n", " m = folium.Map(\n", " location=[np.median(df.latitude), np.median(df.longitude)],\n", " zoom_start=zoom\n", " )\n", "\n", " marker_cluster = MarkerCluster().add_to(m)\n", "\n", " for j in range(k):\n", " for r in df[:-1].iterrows():\n", " if r[1]['label'] == j:\n", " pop = [\"Price: \",str(r[1]['type_A'])]\n", " folium.Marker(\n", " location = [r[1]['latitude'],r[1]['longitude']],\n", " popup=''.join(pop),\n", " icon=folium.Icon(color=cluster_colors[j], icon='ok-sign'),\n", " ).add_to(marker_cluster)\n", " \n", " if custom == 1:\n", " folium.Marker(location = [extralat,extralong],\n", " popup = 'You are here',\n", " icon=folium.Icon(color='pink', icon='cloud')).add_to(m)\n", " \n", " display(m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the 4 main functions to get a json file of addresses and prices, and a list of coordinates." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#address_main()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#clean_addresses()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#price_main()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#coords_main()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "712" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with open('gas_stations_2017-12-13.json', 'r') as infile:\n", " gas_stations = json.load(infile)\n", "len(gas_stations)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Clean up\n", "\n", "For some reason, there's a station near Vaughn that only registers in Google Maps if we input it as King City. There are also a few duplicate stations that are listed under different names or address variants, so we get rid of those." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'address': '12001 ON-400 #1S1, King City, ON L7B 1A8, Canada',\n", " 'latitude': 43.8950319,\n", " 'longitude': -79.5574951,\n", " 'name': 'Canadian Tire',\n", " 'type_A': 125.3,\n", " 'type_D': 119.7}" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Special treatment for 12001 hwy 400 nb, since gasbuddy has it as in Vaughan but Google recognizes it as being in King City\n", "gas_stations['12001 Hwy 400 NB , Vaughan, Ontario']" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "geoloc = geolocator.geocode('12001 Hwy 400 NB , King City, Ontario')" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "gas_stations['12001 Hwy 400 NB , Vaughan, Ontario']['address'] = geoloc.address\n", "gas_stations['12001 Hwy 400 NB , Vaughan, Ontario']['longitude'] = geoloc.longitude\n", "gas_stations['12001 Hwy 400 NB , Vaughan, Ontario']['latitude'] =geoloc.latitude" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "gas_df = pd.DataFrame.from_dict(gas_stations, orient='index')" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitude
12731 Hwy 48 , Whitchurch-Stouffville, OntarioUltramar115.6129.6132.6115.912731 ON-48, Whitchurch-Stouffville, ON L4A 7X...-79.28237843.976346
1610 Keele St , Toronto - Central, OntarioShell122.9136.9144.9NaN1610 Keele St, York, ON M6M 3V9, Canada-79.47198843.682320
6897 Finch Ave W , Toronto - West, OntarioEsso121.9132.9142.9118.96897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada-79.61736443.734632
1 Harwood Ave S , Ajax, OntarioPioneer122.6123.4130.9118.61 Harwood Ave S, Ajax, ON L1S 2C1, Canada-79.02510043.861027
1 Thornhill Woods Dr , Vaughan, OntarioEsso122.9NaNNaN118.91 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C...-79.46388643.826660
\n", "
" ], "text/plain": [ " name type_A type_B \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n", " 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n", " 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n", "1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n", "1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n", "\n", " type_C type_D \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n", " 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n", " 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n", "1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n", "1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n", "\n", " address \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n", " 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n", " 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n", "1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n", "\n", " longitude latitude \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 " ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df.head()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitude
995 Pape Ave , East York, OntarioCocomile Service Centre115.9NaNNaNNaN995 Pape Ave, East York, ON M4K 3V8, Canada-79.34783143.687761
995 Pape Ave , Toronto - Central, OntarioNational Gas & Variety115.9NaNNaNNaN995 Pape Ave, East York, ON M4K 3V8, Canada-79.34783143.687761
\n", "
" ], "text/plain": [ " name type_A \\\n", "995 Pape Ave , East York, Ontario Cocomile Service Centre 115.9 \n", "995 Pape Ave , Toronto - Central, Ontario National Gas & Variety 115.9 \n", "\n", " type_B type_C type_D \\\n", "995 Pape Ave , East York, Ontario NaN NaN NaN \n", "995 Pape Ave , Toronto - Central, Ontario NaN NaN NaN \n", "\n", " address \\\n", "995 Pape Ave , East York, Ontario 995 Pape Ave, East York, ON M4K 3V8, Canada \n", "995 Pape Ave , Toronto - Central, Ontario 995 Pape Ave, East York, ON M4K 3V8, Canada \n", "\n", " longitude latitude \n", "995 Pape Ave , East York, Ontario -79.347831 43.687761 \n", "995 Pape Ave , Toronto - Central, Ontario -79.347831 43.687761 " ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df[gas_df['longitude']==-79.347831]" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitude
490 Great Lakes Dr , Brampton, ON, OntarioShell118.9132.9140.9118.9490 Great Lakes Dr, Brampton, ON L6R, Canada-79.77364543.740888
490 Great Lakes Dr , Brampton, OntarioShell118.9132.9140.9118.9490 Great Lakes Dr, Brampton, ON L6R, Canada-79.77364543.740888
\n", "
" ], "text/plain": [ " name type_A type_B type_C \\\n", "490 Great Lakes Dr , Brampton, ON, Ontario Shell 118.9 132.9 140.9 \n", "490 Great Lakes Dr , Brampton, Ontario Shell 118.9 132.9 140.9 \n", "\n", " type_D \\\n", "490 Great Lakes Dr , Brampton, ON, Ontario 118.9 \n", "490 Great Lakes Dr , Brampton, Ontario 118.9 \n", "\n", " address \\\n", "490 Great Lakes Dr , Brampton, ON, Ontario 490 Great Lakes Dr, Brampton, ON L6R, Canada \n", "490 Great Lakes Dr , Brampton, Ontario 490 Great Lakes Dr, Brampton, ON L6R, Canada \n", "\n", " longitude latitude \n", "490 Great Lakes Dr , Brampton, ON, Ontario -79.773645 43.740888 \n", "490 Great Lakes Dr , Brampton, Ontario -79.773645 43.740888 " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df[gas_df['longitude']==-79.7736453]" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Get rid of redundancies\n", "gas_df.drop_duplicates(subset=['longitude','latitude'], inplace=True)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitude
995 Pape Ave , East York, OntarioCocomile Service Centre115.9NaNNaNNaN995 Pape Ave, East York, ON M4K 3V8, Canada-79.34783143.687761
\n", "
" ], "text/plain": [ " name type_A type_B \\\n", "995 Pape Ave , East York, Ontario Cocomile Service Centre 115.9 NaN \n", "\n", " type_C type_D \\\n", "995 Pape Ave , East York, Ontario NaN NaN \n", "\n", " address \\\n", "995 Pape Ave , East York, Ontario 995 Pape Ave, East York, ON M4K 3V8, Canada \n", "\n", " longitude latitude \n", "995 Pape Ave , East York, Ontario -79.347831 43.687761 " ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df[gas_df['longitude']==-79.347831]" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitude
490 Great Lakes Dr , Brampton, ON, OntarioShell118.9132.9140.9118.9490 Great Lakes Dr, Brampton, ON L6R, Canada-79.77364543.740888
\n", "
" ], "text/plain": [ " name type_A type_B type_C \\\n", "490 Great Lakes Dr , Brampton, ON, Ontario Shell 118.9 132.9 140.9 \n", "\n", " type_D \\\n", "490 Great Lakes Dr , Brampton, ON, Ontario 118.9 \n", "\n", " address \\\n", "490 Great Lakes Dr , Brampton, ON, Ontario 490 Great Lakes Dr, Brampton, ON L6R, Canada \n", "\n", " longitude latitude \n", "490 Great Lakes Dr , Brampton, ON, Ontario -79.773645 43.740888 " ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df[gas_df['longitude']==-79.7736453]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mean prices by brand\n", "The top prices may be unreliable, since they are not major chains or might be out of the way, and so user data on their prices is scarce. Not surprisingly, the cheapest brand is Costco. It's a members only service, and in some sense doesn't really count. It's prices are so low that in the K Means clustering below, the price discrepancy largely overpowers any distance effects." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
countmeanstdmin25%50%75%max
name
Speedy Gas1.0125.900000NaN125.9125.9125.9125.9125.9
Good Guys Gas Bar1.0122.900000NaN122.9122.9122.9122.9122.9
Falcon1.0122.900000NaN122.9122.9122.9122.9122.9
Superstore1.0122.600000NaN122.6122.6122.6122.6122.6
Fas Gas2.0121.9000000.000000121.9121.9121.9121.9121.9
Petro-Canada188.0121.8627662.145533113.6122.7122.9122.9123.3
Esso173.0121.7127172.215133113.2121.3122.9122.9123.9
Shell117.0121.3410262.771970112.4119.9122.9122.9123.9
\n", "
" ], "text/plain": [ " count mean std min 25% 50% 75% \\\n", "name \n", "Speedy Gas 1.0 125.900000 NaN 125.9 125.9 125.9 125.9 \n", "Good Guys Gas Bar 1.0 122.900000 NaN 122.9 122.9 122.9 122.9 \n", "Falcon 1.0 122.900000 NaN 122.9 122.9 122.9 122.9 \n", "Superstore 1.0 122.600000 NaN 122.6 122.6 122.6 122.6 \n", "Fas Gas 2.0 121.900000 0.000000 121.9 121.9 121.9 121.9 \n", "Petro-Canada 188.0 121.862766 2.145533 113.6 122.7 122.9 122.9 \n", "Esso 173.0 121.712717 2.215133 113.2 121.3 122.9 122.9 \n", "Shell 117.0 121.341026 2.771970 112.4 119.9 122.9 122.9 \n", "\n", " max \n", "name \n", "Speedy Gas 125.9 \n", "Good Guys Gas Bar 122.9 \n", "Falcon 122.9 \n", "Superstore 122.6 \n", "Fas Gas 121.9 \n", "Petro-Canada 123.3 \n", "Esso 123.9 \n", "Shell 123.9 " ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df.groupby('name')['type_A'].describe().sort_values(by='mean', ascending=False)[:8]" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
countmeanstdmin25%50%75%max
name
Costco7.0111.4714291.511858110.9110.90110.9110.9114.9
Danforth Gas & Wash1.0112.900000NaN112.9112.90112.9112.9112.9
Star Gas1.0113.000000NaN113.0113.00113.0113.0113.0
Mobil1.0113.900000NaN113.9113.90113.9113.9113.9
Cocomile Service Centre1.0115.900000NaN115.9115.90115.9115.9115.9
OP1.0116.500000NaN116.5116.50116.5116.5116.5
Santak Pump1.0116.500000NaN116.5116.50116.5116.5116.5
Race Trac3.0116.6666672.926317114.2115.05115.9117.9119.9
\n", "
" ], "text/plain": [ " count mean std min 25% 50% \\\n", "name \n", "Costco 7.0 111.471429 1.511858 110.9 110.90 110.9 \n", "Danforth Gas & Wash 1.0 112.900000 NaN 112.9 112.90 112.9 \n", "Star Gas 1.0 113.000000 NaN 113.0 113.00 113.0 \n", "Mobil 1.0 113.900000 NaN 113.9 113.90 113.9 \n", "Cocomile Service Centre 1.0 115.900000 NaN 115.9 115.90 115.9 \n", "OP 1.0 116.500000 NaN 116.5 116.50 116.5 \n", "Santak Pump 1.0 116.500000 NaN 116.5 116.50 116.5 \n", "Race Trac 3.0 116.666667 2.926317 114.2 115.05 115.9 \n", "\n", " 75% max \n", "name \n", "Costco 110.9 114.9 \n", "Danforth Gas & Wash 112.9 112.9 \n", "Star Gas 113.0 113.0 \n", "Mobil 113.9 113.9 \n", "Cocomile Service Centre 115.9 115.9 \n", "OP 116.5 116.5 \n", "Santak Pump 116.5 116.5 \n", "Race Trac 117.9 119.9 " ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df.groupby('name')['type_A'].describe().sort_values(by='mean', ascending=True)[:8]" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitude
12731 Hwy 48 , Whitchurch-Stouffville, OntarioUltramar115.6129.6132.6115.912731 ON-48, Whitchurch-Stouffville, ON L4A 7X...-79.28237843.976346
1610 Keele St , Toronto - Central, OntarioShell122.9136.9144.9NaN1610 Keele St, York, ON M6M 3V9, Canada-79.47198843.682320
6897 Finch Ave W , Toronto - West, OntarioEsso121.9132.9142.9118.96897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada-79.61736443.734632
1 Harwood Ave S , Ajax, OntarioPioneer122.6123.4130.9118.61 Harwood Ave S, Ajax, ON L1S 2C1, Canada-79.02510043.861027
1 Thornhill Woods Dr , Vaughan, OntarioEsso122.9NaNNaN118.91 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C...-79.46388643.826660
\n", "
" ], "text/plain": [ " name type_A type_B \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n", " 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n", " 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n", "1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n", "1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n", "\n", " type_C type_D \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n", " 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n", " 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n", "1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n", "1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n", "\n", " address \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n", " 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n", " 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n", "1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n", "\n", " longitude latitude \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 " ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df.iloc[:5]" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": true }, "outputs": [], "source": [ "gas_df.to_csv('gas_stations_2017-12-13.csv')" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": true }, "outputs": [], "source": [ "with open('gas_stations_2017-12-13.json','w') as outfile:\n", " json.dump(gas_stations,outfile) " ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": true }, "outputs": [], "source": [ "coords_df = gas_df[['longitude','latitude']].copy()" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
longitudelatitude
12731 Hwy 48 , Whitchurch-Stouffville, Ontario-79.28237843.976346
1610 Keele St , Toronto - Central, Ontario-79.47198843.682320
6897 Finch Ave W , Toronto - West, Ontario-79.61736443.734632
1 Harwood Ave S , Ajax, Ontario-79.02510043.861027
1 Thornhill Woods Dr , Vaughan, Ontario-79.46388643.826660
\n", "
" ], "text/plain": [ " longitude latitude\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346\n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320\n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632\n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027\n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coords_df.head()" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": true }, "outputs": [], "source": [ "coords_df['long_rad'] = coords_df['longitude'].apply(lambda l:math.radians(l))\n", "coords_df['lat_rad'] = coords_df['latitude'].apply(lambda l:math.radians(l))" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
longitudelatitudelong_radlat_rad
12731 Hwy 48 , Whitchurch-Stouffville, Ontario-79.28237843.976346-1.3837390.767532
1610 Keele St , Toronto - Central, Ontario-79.47198843.682320-1.3870480.762400
6897 Finch Ave W , Toronto - West, Ontario-79.61736443.734632-1.3895850.763313
1 Harwood Ave S , Ajax, Ontario-79.02510043.861027-1.3792480.765519
1 Thornhill Woods Dr , Vaughan, Ontario-79.46388643.826660-1.3869060.764920
\n", "
" ], "text/plain": [ " longitude latitude \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n", "\n", " long_rad lat_rad \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -1.383739 0.767532 \n", " 1610 Keele St , Toronto - Central, Ontario -1.387048 0.762400 \n", " 6897 Finch Ave W , Toronto - West, Ontario -1.389585 0.763313 \n", "1 Harwood Ave S , Ajax, Ontario -1.379248 0.765519 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -1.386906 0.764920 " ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coords_df.head()" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": true }, "outputs": [], "source": [ "coords_df['combined'] = list(zip(coords_df['longitude'],coords_df['latitude']))" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
longitudelatitudelong_radlat_rad
combined
(-79.877743, 43.6495942)1111
(-79.4275737, 43.6709778)1111
(-79.41025359999999, 43.8575578)1111
(-79.4099732, 43.75732610000001)1111
(-79.4098422, 43.8571536)1111
\n", "
" ], "text/plain": [ " longitude latitude long_rad lat_rad\n", "combined \n", "(-79.877743, 43.6495942) 1 1 1 1\n", "(-79.4275737, 43.6709778) 1 1 1 1\n", "(-79.41025359999999, 43.8575578) 1 1 1 1\n", "(-79.4099732, 43.75732610000001) 1 1 1 1\n", "(-79.4098422, 43.8571536) 1 1 1 1" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coords_df.groupby('combined').count().sort_values('longitude',ascending=False)[:5]" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitude
490 Great Lakes Dr , Brampton, ON, OntarioShell118.9132.9140.9118.9490 Great Lakes Dr, Brampton, ON L6R, Canada-79.77364543.740888
\n", "
" ], "text/plain": [ " name type_A type_B type_C \\\n", "490 Great Lakes Dr , Brampton, ON, Ontario Shell 118.9 132.9 140.9 \n", "\n", " type_D \\\n", "490 Great Lakes Dr , Brampton, ON, Ontario 118.9 \n", "\n", " address \\\n", "490 Great Lakes Dr , Brampton, ON, Ontario 490 Great Lakes Dr, Brampton, ON L6R, Canada \n", "\n", " longitude latitude \n", "490 Great Lakes Dr , Brampton, ON, Ontario -79.773645 43.740888 " ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df[coords_df['combined']==(-79.7736453, 43.7408885)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Checking for gas stations with competitors within 100 m\n", "\n", "While competition pushes prices down, especially when there is almost 100% price transparency (all gas stations advertise their regular gas prices prominently), multiple gas stations are usually at very busy intersections. In that case, can the extra demand counteract the competitive effect? We can check." ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": true }, "outputs": [], "source": [ "within_100 = []\n", "for i in range(len(coords_df)):\n", " closest_coord = xy_closest(coords_df['combined'][coords_df.index != coords_df.index[i]], coords_df.iloc[i]['combined'])\n", " within_100.append(xy_distance(coords_df.iloc[i]['combined'][0],coords_df.iloc[i]['combined'][1], closest_coord[0], closest_coord[1]) <= 0.1)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "133" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum(within_100)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": true }, "outputs": [], "source": [ "coords_df['within_100'] = np.asarray(within_100).astype(int)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
longitudelatitudelong_radlat_radcombinedwithin_100
12731 Hwy 48 , Whitchurch-Stouffville, Ontario-79.28237843.976346-1.3837390.767532(-79.282378, 43.9763464)0
1610 Keele St , Toronto - Central, Ontario-79.47198843.682320-1.3870480.762400(-79.4719879, 43.6823199)0
6897 Finch Ave W , Toronto - West, Ontario-79.61736443.734632-1.3895850.763313(-79.6173637, 43.7346325)0
1 Harwood Ave S , Ajax, Ontario-79.02510043.861027-1.3792480.765519(-79.0251004, 43.8610274)1
1 Thornhill Woods Dr , Vaughan, Ontario-79.46388643.826660-1.3869060.764920(-79.46388639999999, 43.8266602)0
\n", "
" ], "text/plain": [ " longitude latitude \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n", "\n", " long_rad lat_rad \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -1.383739 0.767532 \n", " 1610 Keele St , Toronto - Central, Ontario -1.387048 0.762400 \n", " 6897 Finch Ave W , Toronto - West, Ontario -1.389585 0.763313 \n", "1 Harwood Ave S , Ajax, Ontario -1.379248 0.765519 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -1.386906 0.764920 \n", "\n", " combined \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario (-79.282378, 43.9763464) \n", " 1610 Keele St , Toronto - Central, Ontario (-79.4719879, 43.6823199) \n", " 6897 Finch Ave W , Toronto - West, Ontario (-79.6173637, 43.7346325) \n", "1 Harwood Ave S , Ajax, Ontario (-79.0251004, 43.8610274) \n", "1 Thornhill Woods Dr , Vaughan, Ontario (-79.46388639999999, 43.8266602) \n", "\n", " within_100 \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n", " 1610 Keele St , Toronto - Central, Ontario 0 \n", " 6897 Finch Ave W , Toronto - West, Ontario 0 \n", "1 Harwood Ave S , Ajax, Ontario 1 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 0 " ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coords_df.head()" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": true }, "outputs": [], "source": [ "gas_df['within_100'] = coords_df['within_100']" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitudewithin_100
12731 Hwy 48 , Whitchurch-Stouffville, OntarioUltramar115.6129.6132.6115.912731 ON-48, Whitchurch-Stouffville, ON L4A 7X...-79.28237843.9763460
1610 Keele St , Toronto - Central, OntarioShell122.9136.9144.9NaN1610 Keele St, York, ON M6M 3V9, Canada-79.47198843.6823200
6897 Finch Ave W , Toronto - West, OntarioEsso121.9132.9142.9118.96897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada-79.61736443.7346320
1 Harwood Ave S , Ajax, OntarioPioneer122.6123.4130.9118.61 Harwood Ave S, Ajax, ON L1S 2C1, Canada-79.02510043.8610271
1 Thornhill Woods Dr , Vaughan, OntarioEsso122.9NaNNaN118.91 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C...-79.46388643.8266600
\n", "
" ], "text/plain": [ " name type_A type_B \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n", " 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n", " 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n", "1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n", "1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n", "\n", " type_C type_D \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n", " 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n", " 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n", "1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n", "1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n", "\n", " address \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n", " 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n", " 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n", "1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n", "\n", " longitude latitude \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n", "\n", " within_100 \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n", " 1610 Keele St , Toronto - Central, Ontario 0 \n", " 6897 Finch Ave W , Toronto - West, Ontario 0 \n", "1 Harwood Ave S , Ajax, Ontario 1 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 0 " ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For general clustering, use longitude and latitude. For gasoline-specific clustering, also use the price of the specific type of gas. We can just use k means for general clustering, but when looking at only certain types of gas (such as regular type_A)." ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": true }, "outputs": [], "source": [ "cluster_df = gas_df[['name','type_A','longitude','latitude']].copy()" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Alongitudelatitude
12731 Hwy 48 , Whitchurch-Stouffville, OntarioUltramar115.6-79.28237843.976346
1610 Keele St , Toronto - Central, OntarioShell122.9-79.47198843.682320
6897 Finch Ave W , Toronto - West, OntarioEsso121.9-79.61736443.734632
1 Harwood Ave S , Ajax, OntarioPioneer122.6-79.02510043.861027
1 Thornhill Woods Dr , Vaughan, OntarioEsso122.9-79.46388643.826660
\n", "
" ], "text/plain": [ " name type_A \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 \n", " 1610 Keele St , Toronto - Central, Ontario Shell 122.9 \n", " 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 \n", "1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 \n", "1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 \n", "\n", " longitude latitude \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 " ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cluster_df.head()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### Exploring the Starbucks Metric\n", "Use Starbucks' location research to our benefit. As with station proximities to other stations, we will use the \"as the crow flies\" distance.\n", "\n", "Starbucks location data from https://opendata.socrata.com/Business/All-Starbucks-Locations-in-the-World/xy4y-c4mk" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
indexStore IDNameBrandStore NumberPhone NumberOwnership TypeStreet CombinedStreet 1Street 2...Country SubdivisionCountryPostal CodeCoordinatesLatitudeLongitudeTimezoneCurrent Timezone OffsetOlson TimezoneFirst Seen
0108323King & BrantStarbucks4636-99652416-596-0101CO527-529 King Street West527-529 King Street WestNaN...ONCAM5V 1K4(43.6446495056152, -79.3982620239258)43.644650-79.398262Eastern Standard Time-300GMT-05:00 America/Toronto12/08/2013 10:41:59 PM
11223353035 Argentia RdStarbucks4665-100880905-785-7667CO3035 Argentia Road, unIT 73035 Argentia RoadunIT 7...ONCAL5N 8P7(43.5956649780273, -79.7863235473633)43.595665-79.786324Eastern Standard Time-300GMT-05:00 America/Toronto12/08/2013 10:41:59 PM
2124337631 Commissioners Rd EStarbucks4674-100799519-649-1444CO631 Commissioners Road East631 Commissioners Road EastNaN...ONCAN6C 2V1(42.957763671875, -81.2328872680664)42.957764-81.232887Eastern Standard Time-300GMT-05:00 America/Toronto12/08/2013 10:41:59 PM
3133354599 Taylor Kidd BlvdStarbucks4694-102072613-634-1509CO599 Taylor Kidd Blvd, Unit 6599 Taylor Kidd BlvdUnit 6...ONCAK7M 0A2(44.2504081726074, -76.5673141479492)44.250408-76.567314Eastern Standard Time-300GMT-05:00 America/Toronto12/08/2013 10:41:59 PM
4141364King & ShawStarbucks4646-98483416-979-1953CO1005 King St W, Unit 71005 King St WUnit 7...ONCAM6G 1B9(43.6412811279297, -79.4147644042969)43.641281-79.414764Eastern Standard Time-300GMT-05:00 America/Toronto12/08/2013 10:41:59 PM
\n", "

5 rows × 22 columns

\n", "
" ], "text/plain": [ " index Store ID Name Brand Store Number \\\n", "0 108 323 King & Brant Starbucks 4636-99652 \n", "1 122 335 3035 Argentia Rd Starbucks 4665-100880 \n", "2 124 337 631 Commissioners Rd E Starbucks 4674-100799 \n", "3 133 354 599 Taylor Kidd Blvd Starbucks 4694-102072 \n", "4 141 364 King & Shaw Starbucks 4646-98483 \n", "\n", " Phone Number Ownership Type Street Combined \\\n", "0 416-596-0101 CO 527-529 King Street West \n", "1 905-785-7667 CO 3035 Argentia Road, unIT 7 \n", "2 519-649-1444 CO 631 Commissioners Road East \n", "3 613-634-1509 CO 599 Taylor Kidd Blvd, Unit 6 \n", "4 416-979-1953 CO 1005 King St W, Unit 7 \n", "\n", " Street 1 Street 2 ... \\\n", "0 527-529 King Street West NaN ... \n", "1 3035 Argentia Road unIT 7 ... \n", "2 631 Commissioners Road East NaN ... \n", "3 599 Taylor Kidd Blvd Unit 6 ... \n", "4 1005 King St W Unit 7 ... \n", "\n", " Country Subdivision Country Postal Code \\\n", "0 ON CA M5V 1K4 \n", "1 ON CA L5N 8P7 \n", "2 ON CA N6C 2V1 \n", "3 ON CA K7M 0A2 \n", "4 ON CA M6G 1B9 \n", "\n", " Coordinates Latitude Longitude \\\n", "0 (43.6446495056152, -79.3982620239258) 43.644650 -79.398262 \n", "1 (43.5956649780273, -79.7863235473633) 43.595665 -79.786324 \n", "2 (42.957763671875, -81.2328872680664) 42.957764 -81.232887 \n", "3 (44.2504081726074, -76.5673141479492) 44.250408 -76.567314 \n", "4 (43.6412811279297, -79.4147644042969) 43.641281 -79.414764 \n", "\n", " Timezone Current Timezone Offset Olson Timezone \\\n", "0 Eastern Standard Time -300 GMT-05:00 America/Toronto \n", "1 Eastern Standard Time -300 GMT-05:00 America/Toronto \n", "2 Eastern Standard Time -300 GMT-05:00 America/Toronto \n", "3 Eastern Standard Time -300 GMT-05:00 America/Toronto \n", "4 Eastern Standard Time -300 GMT-05:00 America/Toronto \n", "\n", " First Seen \n", "0 12/08/2013 10:41:59 PM \n", "1 12/08/2013 10:41:59 PM \n", "2 12/08/2013 10:41:59 PM \n", "3 12/08/2013 10:41:59 PM \n", "4 12/08/2013 10:41:59 PM \n", "\n", "[5 rows x 22 columns]" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "starbucks_df = pd.read_csv('All_Starbucks_Locations_in_the_World.csv', delimiter=';')\n", "starbucks_df = starbucks_df[(starbucks_df['Country']=='CA') & (starbucks_df['Country Subdivision']=='ON')]\n", "starbucks_df.reset_index(inplace=True)\n", "starbucks_df.head()" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "RangeIndex: 530 entries, 0 to 529\n", "Data columns (total 22 columns):\n", "index 530 non-null int64\n", "Store ID 530 non-null object\n", "Name 530 non-null object\n", "Brand 530 non-null object\n", "Store Number 530 non-null object\n", "Phone Number 486 non-null object\n", "Ownership Type 530 non-null object\n", "Street Combined 530 non-null object\n", "Street 1 530 non-null object\n", "Street 2 191 non-null object\n", "Street 3 110 non-null object\n", "City 530 non-null object\n", "Country Subdivision 530 non-null object\n", "Country 530 non-null object\n", "Postal Code 530 non-null object\n", "Coordinates 530 non-null object\n", "Latitude 530 non-null float64\n", "Longitude 530 non-null float64\n", "Timezone 530 non-null object\n", "Current Timezone Offset 530 non-null int64\n", "Olson Timezone 530 non-null object\n", "First Seen 530 non-null object\n", "dtypes: float64(2), int64(2), object(18)\n", "memory usage: 91.2+ KB\n" ] } ], "source": [ "starbucks_df.info()" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Actually we only need the coordinates.\n", "starbucks_on_coords = starbucks_df['Coordinates'].tolist()" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Convert from a list of strings that look like Python lists to an actual list of lists\n", "# Also need to change it from (lat,long) to (long,lat)\n", "starbucks_on_coords = [ast.literal_eval(starbucks_on_coords[i]) for i in range(len(starbucks_on_coords))]\n", "starbucks_on_coords = [(starbucks_on_coords[i][1], starbucks_on_coords[i][0]) for i in range(len(starbucks_on_coords))]" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": true }, "outputs": [], "source": [ "nearest_starbucks = []\n", "for i in range(len(coords_df)):\n", " closest_coord = xy_closest(starbucks_on_coords, coords_df.iloc[i]['combined'])\n", " nearest_starbucks.append(xy_distance(coords_df.iloc[i]['combined'][1],coords_df.iloc[i]['combined'][0], closest_coord[1], closest_coord[0]))" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": true }, "outputs": [], "source": [ "starbucks_within_100m = []\n", "for i in range(len(coords_df)):\n", " closest_coord = xy_closest(starbucks_on_coords, coords_df.iloc[i]['combined'])\n", " starbucks_within_100m.append(xy_distance(coords_df.iloc[i]['combined'][1],coords_df.iloc[i]['combined'][0], closest_coord[1], closest_coord[0])<=1)" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": true }, "outputs": [], "source": [ "gas_df['nearest_starbucks'] = np.asarray(nearest_starbucks)\n", "gas_df['starbucks_within_100m'] = np.asarray(starbucks_within_100m).astype(int)" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitudewithin_100nearest_starbucksstarbucks_within_100m
12731 Hwy 48 , Whitchurch-Stouffville, OntarioUltramar115.6129.6132.6115.912731 ON-48, Whitchurch-Stouffville, ON L4A 7X...-79.28237843.976346011.2086290
1610 Keele St , Toronto - Central, OntarioShell122.9136.9144.9NaN1610 Keele St, York, ON M6M 3V9, Canada-79.47198843.68232001.1782520
6897 Finch Ave W , Toronto - West, OntarioEsso121.9132.9142.9118.96897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada-79.61736443.73463201.0316030
1 Harwood Ave S , Ajax, OntarioPioneer122.6123.4130.9118.61 Harwood Ave S, Ajax, ON L1S 2C1, Canada-79.02510043.86102710.3816521
1 Thornhill Woods Dr , Vaughan, OntarioEsso122.9NaNNaN118.91 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C...-79.46388643.82666001.4445660
\n", "
" ], "text/plain": [ " name type_A type_B \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n", " 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n", " 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n", "1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n", "1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n", "\n", " type_C type_D \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n", " 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n", " 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n", "1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n", "1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n", "\n", " address \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n", " 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n", " 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n", "1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n", "\n", " longitude latitude \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n", "\n", " within_100 \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n", " 1610 Keele St , Toronto - Central, Ontario 0 \n", " 6897 Finch Ave W , Toronto - West, Ontario 0 \n", "1 Harwood Ave S , Ajax, Ontario 1 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 0 \n", "\n", " nearest_starbucks \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 11.208629 \n", " 1610 Keele St , Toronto - Central, Ontario 1.178252 \n", " 6897 Finch Ave W , Toronto - West, Ontario 1.031603 \n", "1 Harwood Ave S , Ajax, Ontario 0.381652 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1.444566 \n", "\n", " starbucks_within_100m \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n", " 1610 Keele St , Toronto - Central, Ontario 0 \n", " 6897 Finch Ave W , Toronto - West, Ontario 0 \n", "1 Harwood Ave S , Ajax, Ontario 1 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 0 " ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's invent a Starbucks proximity metric, a sigmoid with the log of the inverse of the nearest starbucks in km:
\n", "$\\frac{1}{(1+exp(log(-\\frac{1}{d})^2))}$\n", "(This is not actually used, but keep this here just in case)" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": true }, "outputs": [], "source": [ "gas_df['starbucks_metric'] = gas_df['nearest_starbucks'].apply(lambda d:1/(1+np.exp(-np.log((2/d)**2))))" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Atype_Btype_Ctype_Daddresslongitudelatitudewithin_100nearest_starbucksstarbucks_within_100mstarbucks_metric
12731 Hwy 48 , Whitchurch-Stouffville, OntarioUltramar115.6129.6132.6115.912731 ON-48, Whitchurch-Stouffville, ON L4A 7X...-79.28237843.976346011.20862900.030856
1610 Keele St , Toronto - Central, OntarioShell122.9136.9144.9NaN1610 Keele St, York, ON M6M 3V9, Canada-79.47198843.68232001.17825200.742352
6897 Finch Ave W , Toronto - West, OntarioEsso121.9132.9142.9118.96897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada-79.61736443.73463201.03160300.789857
1 Harwood Ave S , Ajax, OntarioPioneer122.6123.4130.9118.61 Harwood Ave S, Ajax, ON L1S 2C1, Canada-79.02510043.86102710.38165210.964865
1 Thornhill Woods Dr , Vaughan, OntarioEsso122.9NaNNaN118.91 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C...-79.46388643.82666001.44456600.657163
\n", "
" ], "text/plain": [ " name type_A type_B \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 129.6 \n", " 1610 Keele St , Toronto - Central, Ontario Shell 122.9 136.9 \n", " 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 132.9 \n", "1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 123.4 \n", "1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 NaN \n", "\n", " type_C type_D \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 132.6 115.9 \n", " 1610 Keele St , Toronto - Central, Ontario 144.9 NaN \n", " 6897 Finch Ave W , Toronto - West, Ontario 142.9 118.9 \n", "1 Harwood Ave S , Ajax, Ontario 130.9 118.6 \n", "1 Thornhill Woods Dr , Vaughan, Ontario NaN 118.9 \n", "\n", " address \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 12731 ON-48, Whitchurch-Stouffville, ON L4A 7X... \n", " 1610 Keele St , Toronto - Central, Ontario 1610 Keele St, York, ON M6M 3V9, Canada \n", " 6897 Finch Ave W , Toronto - West, Ontario 6897 Finch Ave W, Etobicoke, ON M9W 0A6, Canada \n", "1 Harwood Ave S , Ajax, Ontario 1 Harwood Ave S, Ajax, ON L1S 2C1, Canada \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1 Thornhill Woods Dr, Thornhill, ON L4J 8Y2, C... \n", "\n", " longitude latitude \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n", "\n", " within_100 \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n", " 1610 Keele St , Toronto - Central, Ontario 0 \n", " 6897 Finch Ave W , Toronto - West, Ontario 0 \n", "1 Harwood Ave S , Ajax, Ontario 1 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 0 \n", "\n", " nearest_starbucks \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 11.208629 \n", " 1610 Keele St , Toronto - Central, Ontario 1.178252 \n", " 6897 Finch Ave W , Toronto - West, Ontario 1.031603 \n", "1 Harwood Ave S , Ajax, Ontario 0.381652 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1.444566 \n", "\n", " starbucks_within_100m \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n", " 1610 Keele St , Toronto - Central, Ontario 0 \n", " 6897 Finch Ave W , Toronto - West, Ontario 0 \n", "1 Harwood Ave S , Ajax, Ontario 1 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 0 \n", "\n", " starbucks_metric \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0.030856 \n", " 1610 Keele St , Toronto - Central, Ontario 0.742352 \n", " 6897 Finch Ave W , Toronto - West, Ontario 0.789857 \n", "1 Harwood Ave S , Ajax, Ontario 0.964865 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 0.657163 " ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gas_df.head()" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": true }, "outputs": [], "source": [ "gas_df.to_csv('gas_stations_2017-12-13.csv')" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": true }, "outputs": [], "source": [ "regular_gas = gas_df[['type_A','within_100','starbucks_within_100m']].dropna().copy()" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": true }, "outputs": [], "source": [ "regular_gas.columns=['Price','Competition_Nearby','Starbucks_Nearby']" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": true }, "outputs": [], "source": [ "regular_gas['Relative_Price'] = regular_gas['Price'].apply(lambda p: p/np.mean(regular_gas['Price']))\n", "regular_gas['Price_Difference'] = regular_gas['Price'].apply(lambda p: p-np.mean(regular_gas['Price']))" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
PriceCompetition_NearbyStarbucks_NearbyRelative_PricePrice_Difference
12731 Hwy 48 , Whitchurch-Stouffville, Ontario115.6000.956242-5.289926
1610 Keele St , Toronto - Central, Ontario122.9001.0166272.010074
6897 Finch Ave W , Toronto - West, Ontario121.9001.0083551.010074
1 Harwood Ave S , Ajax, Ontario122.6111.0141461.710074
1 Thornhill Woods Dr , Vaughan, Ontario122.9001.0166272.010074
\n", "
" ], "text/plain": [ " Price Competition_Nearby \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 115.6 0 \n", " 1610 Keele St , Toronto - Central, Ontario 122.9 0 \n", " 6897 Finch Ave W , Toronto - West, Ontario 121.9 0 \n", "1 Harwood Ave S , Ajax, Ontario 122.6 1 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 122.9 0 \n", "\n", " Starbucks_Nearby \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0 \n", " 1610 Keele St , Toronto - Central, Ontario 0 \n", " 6897 Finch Ave W , Toronto - West, Ontario 0 \n", "1 Harwood Ave S , Ajax, Ontario 1 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 0 \n", "\n", " Relative_Price \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0.956242 \n", " 1610 Keele St , Toronto - Central, Ontario 1.016627 \n", " 6897 Finch Ave W , Toronto - West, Ontario 1.008355 \n", "1 Harwood Ave S , Ajax, Ontario 1.014146 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1.016627 \n", "\n", " Price_Difference \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -5.289926 \n", " 1610 Keele St , Toronto - Central, Ontario 2.010074 \n", " 6897 Finch Ave W , Toronto - West, Ontario 1.010074 \n", "1 Harwood Ave S , Ajax, Ontario 1.710074 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 2.010074 " ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "regular_gas.head()" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\statsmodels\\compat\\pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n", " from pandas.core import datetools\n" ] } ], "source": [ "import statsmodels.api as sm" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Y = regular_gas['Price']\n", "X = sm.add_constant(regular_gas[['Competition_Nearby','Starbucks_Nearby']])\n", "lm = sm.OLS(Y,X)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "collapsed": true }, "outputs": [], "source": [ "results = lm.fit()" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "const 120.606227\n", "Competition_Nearby 0.182607\n", "Starbucks_Nearby 0.621326\n", "dtype: float64" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results.params" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: Price R-squared: 0.011
Model: OLS Adj. R-squared: 0.008
Method: Least Squares F-statistic: 3.873
Date: Sat, 06 Jan 2018 Prob (F-statistic): 0.0213
Time: 21:37:31 Log-Likelihood: -1681.8
No. Observations: 675 AIC: 3370.
Df Residuals: 672 BIC: 3383.
Df Model: 2
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
const 120.6062 0.155 777.316 0.000 120.302 120.911
Competition_Nearby 0.1826 0.286 0.638 0.523 -0.379 0.744
Starbucks_Nearby 0.6213 0.230 2.699 0.007 0.169 1.073
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 132.421 Durbin-Watson: 1.985
Prob(Omnibus): 0.000 Jarque-Bera (JB): 211.951
Skew: -1.304 Prob(JB): 9.45e-47
Kurtosis: 3.855 Cond. No. 2.90
" ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: Price R-squared: 0.011\n", "Model: OLS Adj. R-squared: 0.008\n", "Method: Least Squares F-statistic: 3.873\n", "Date: Sat, 06 Jan 2018 Prob (F-statistic): 0.0213\n", "Time: 21:37:31 Log-Likelihood: -1681.8\n", "No. Observations: 675 AIC: 3370.\n", "Df Residuals: 672 BIC: 3383.\n", "Df Model: 2 \n", "Covariance Type: nonrobust \n", "======================================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "--------------------------------------------------------------------------------------\n", "const 120.6062 0.155 777.316 0.000 120.302 120.911\n", "Competition_Nearby 0.1826 0.286 0.638 0.523 -0.379 0.744\n", "Starbucks_Nearby 0.6213 0.230 2.699 0.007 0.169 1.073\n", "==============================================================================\n", "Omnibus: 132.421 Durbin-Watson: 1.985\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 211.951\n", "Skew: -1.304 Prob(JB): 9.45e-47\n", "Kurtosis: 3.855 Cond. No. 2.90\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results.summary()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Gas stations within 100m of a Starbucks location charge 0.62 cents more for regular gas, and this is statistically significant with over 99% confidence.\n", "\n", "Proximity to competition also seems to increase prices, likely due to increased demand at busy intersections (as mentioned above). However, the p value is far too high so we can not be certain of its effect with any confidence." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Clustering and Mapping\n", "\n", "Let's see if we can find any meaningful clusters." ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "collapsed": true }, "outputs": [], "source": [ "cluster_df.dropna(inplace=True)" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Index: 675 entries, 12731 Hwy 48 , Whitchurch-Stouffville, Ontario to Taymall Ave , Toronto - West, Ontario\n", "Data columns (total 4 columns):\n", "name 675 non-null object\n", "type_A 675 non-null float64\n", "longitude 675 non-null float64\n", "latitude 675 non-null float64\n", "dtypes: float64(3), object(1)\n", "memory usage: 26.4+ KB\n" ] } ], "source": [ "cluster_df.info()" ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "collapsed": true }, "outputs": [], "source": [ "cluster_df['type_A_norm'] = cluster_df['type_A']/np.mean(cluster_df['type_A'])\n", "cluster_df['longitude_norm'] = -cluster_df['longitude']/np.mean(cluster_df['longitude'])\n", "cluster_df['latitude_norm'] = cluster_df['latitude']/np.mean(cluster_df['latitude'])" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Alongitudelatitudetype_A_normlongitude_normlatitude_norm
12731 Hwy 48 , Whitchurch-Stouffville, OntarioUltramar115.6-79.28237843.9763460.956242-0.9974351.005572
1610 Keele St , Toronto - Central, OntarioShell122.9-79.47198843.6823201.016627-0.9998200.998849
6897 Finch Ave W , Toronto - West, OntarioEsso121.9-79.61736443.7346321.008355-1.0016491.000045
1 Harwood Ave S , Ajax, OntarioPioneer122.6-79.02510043.8610271.014146-0.9941981.002935
1 Thornhill Woods Dr , Vaughan, OntarioEsso122.9-79.46388643.8266601.016627-0.9997181.002149
\n", "
" ], "text/plain": [ " name type_A \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 \n", " 1610 Keele St , Toronto - Central, Ontario Shell 122.9 \n", " 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 \n", "1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 \n", "1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 \n", "\n", " longitude latitude \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n", "\n", " type_A_norm \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0.956242 \n", " 1610 Keele St , Toronto - Central, Ontario 1.016627 \n", " 6897 Finch Ave W , Toronto - West, Ontario 1.008355 \n", "1 Harwood Ave S , Ajax, Ontario 1.014146 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1.016627 \n", "\n", " longitude_norm \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -0.997435 \n", " 1610 Keele St , Toronto - Central, Ontario -0.999820 \n", " 6897 Finch Ave W , Toronto - West, Ontario -1.001649 \n", "1 Harwood Ave S , Ajax, Ontario -0.994198 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -0.999718 \n", "\n", " latitude_norm \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 1.005572 \n", " 1610 Keele St , Toronto - Central, Ontario 0.998849 \n", " 6897 Finch Ave W , Toronto - West, Ontario 1.000045 \n", "1 Harwood Ave S , Ajax, Ontario 1.002935 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1.002149 " ] }, "execution_count": 63, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cluster_df.head()" ] }, { "cell_type": "code", "execution_count": 64, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from scipy import linalg\n", "from sklearn.cluster import KMeans\n", "from sklearn import svm\n", "from scipy.spatial.distance import cdist\n", "from scipy.spatial import ConvexHull\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5,1,'SSE vs k')" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEWCAYAAABbgYH9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XmcHVWd9/HPNyuEJYQkaCCQZoli\nFAzShs4DuEUQFAgIjsEo+IgEFWbcUEFkRniGGcANGdCZAGqEQKIoEmdEdhhFCDQ7IQJNCBACodlC\nQoBsv+ePU03fdG6vudXV9/b3/Xrd1606dercX41M/3JOVZ2jiMDMzKzSBhQdgJmZ1SYnGDMzy4UT\njJmZ5cIJxszMcuEEY2ZmuXCCMTOzXDjBmNU4SXWSQtKgomOx/sUJxqwHJO0n6W+Slkt6SdJtkt6f\nHRsi6UeSlkhaKekJST8pOXexpNezYy2fC4q7GrN8+F80Zt0kaWvgv4EvA78BhgD7A29mVU4F6oFJ\nwLPAOOADbZo5NCJu6JWAzQriHoxZ970DICKuiIh1EfF6RFwXEQ9kx98PXBURSyNZHBG/7u6PSNo+\n6+lsW1K2l6QXJA2WtJukW7Ne1AuS5nax3SOzXtR7uhuTWXc4wZh136PAOkmzJB0saUSb43cA35D0\nFUl7SFJPfiQilgK3A0eWFH8GuDIi1gD/D7gOGAGMBf6jszYl/V/gHOCjEfFQT+Iy6yonGLNuiohX\ngf2AAC4CmiXNk/S2rMq/k/6ITwcagWckHdummT9IeqXkc3w7P3c5cDRAlqimZWUAa0jDb9tHxBsR\n8ddOQv8a8C3gQxHR1NXrNespebJLs00jaXfgMuCxiDi6zbHNgS8A5wPviYiFkhYDX+zKPRhJ25Du\n4+wKjAcuBcZFREh6O6kX8wngZeBHEfGLMm3UAU8AzwNnRsSFPbxUs25xD8ZsE0XE34FfARvd08ju\nz1xISgATetD2K6RhsH8gDY9dEdm/CiPiuYg4PiK2B04AfiZptw6aOxD4nqQjO6hjVjFOMGbdJGl3\nSd+UNDbb35E0jHVHtv81SR+StLmkQdnw2FbAvT38ycuBY0j3YlqGx5D0qZYYSAksgHUdtLMAOAi4\nUNJhPYzFrMucYMy6bwWwDzBf0mukxPIQ8M3s+OvAj4DngBeAE4EjI2JRSRt/bPMezFUd/N480vDY\nsoi4v6T8/VkMK7M6X42IJzoKPDv/EOAiSQd38XrNesT3YMzMLBfuwZiZWS6cYMzMLBdOMGZmlgsn\nGDMzy0W/nuxy1KhRUVdXV3QYZmZV5e67734hIkZ3Vq9fJ5i6ujoaGxuLDsPMrKpIerIr9TxEZmZm\nuXCCMTOzXDjBmJlZLpxgzMwsF04wZmaWCyeYbpo9G+rqYMCA9D17dtERmZn1Tf36MeXumj0bZsyA\nVavS/pNPpn2A6dOLi8vMrC9yD6YbTjutNbm0WLUqlZuZ2YacYLrhqae6V25m1p85wXTDTjt1r9zM\nrD9zgumGs86CYcM2LBs2LJWbmdmGnGC6Yfp0mDkTxo1L+0OGpH3f4Dcz25gTTDdNnw6LF8PJJ4ME\nRx1VdERmZn1TrglG0kGSHpHUJOmUMseHSpqbHZ8vqS4rP0DS3ZIezL4/UnLO3ll5k6TzJSkr31bS\n9ZIey75H5HltH/0oTJ0Ky5fn+StmZtUrtwQjaSBwIXAwMAE4WtKENtWOA16OiN2AnwDnZOUvAIdG\nxB7AscClJef8HJgBjM8+B2XlpwA3RsR44MZsPzcf+xjMnQvbbZfnr5iZVa88ezCTgKaIWBQRq4E5\nwNQ2daYCs7LtK4EpkhQR90bE0qx8AbBZ1tsZA2wdEbdHRAC/Bg4v09askvJcuQdjZlZenglmB+Dp\nkv0lWVnZOhGxFlgOjGxT50jg3oh4M6u/pJ023xYRz2ZtPQuU7VtImiGpUVJjc3Nzty+q1AknwMSJ\nm9SEmVnNyjPBqExZdKeOpHeThs1O6EabHYqImRFRHxH1o0d3uuJnh975znTD/7nnNqkZM7OalGeC\nWQLsWLI/FljaXh1Jg4DhwEvZ/ljgKuCYiHi8pP7Ydtpclg2hkX0/X7Eracfkyen7jjvy/iUzs+qT\nZ4K5CxgvaWdJQ4BpwLw2deaRbuIDHAXcFBEhaRvgf4BTI+K2lsrZ0NcKSQ3Z02PHAFeXaevYkvLc\n7LUXDB7sBGNmVk5uCSa7p3IScC2wEPhNRCyQdKakw7JqlwAjJTUB36D1ya+TgN2A0yXdl31a7ql8\nGbgYaAIeB67Jys8GDpD0GHBAtp+rzTZLScYJxsxsY0oPY/VP9fX10djYuEltXHklRMCnPlWhoMzM\n+jhJd0dEfWf1vB7MJvKb/GZm5XmqmE0UAfffDwsWFB2JmVnf4gRTAQceCOeeW3QUZmZ9ixPMJpKg\nocE3+s3M2nKCqYDJk+HRR+HFF4uOxMys73CCqYCGhvR9553FxmFm1pc4wVRAfT0MGAC33150JGZm\nfYcfU66ALbeEW2+Fd7+76EjMzPoOJ5gK2W+/oiMwM+tbPERWIUuXwr/+KzzxRNGRmJn1DU4wFbJi\nBZx+Otx0U9GRmJn1DU4wFfKOd8CIEX4fxsyshRNMhfiFSzOzDTnBVFBDQ5qTbPnyoiMxMyueE0wF\nNTTA0KHprX4zs/7OjylX0Ic/DK++mla5NDPr75xgKsiJxcyslYfIKmzOnDR9fz9eKNTMDHCCqbgV\nK+D666GpqehIzMyKlWuCkXSQpEckNUk6pczxoZLmZsfnS6rLykdKulnSSkkXlNTfStJ9JZ8XJJ2X\nHfu8pOaSY1/M89ra0zKzsh9XNrP+LrcEI2kgcCFwMDABOFrShDbVjgNejojdgJ8A52TlbwCnAyeX\nVo6IFRExseUDPAn8vqTK3JLjF1f+qjo3YUKa/NIJxsz6uzx7MJOApohYFBGrgTnA1DZ1pgKzsu0r\ngSmSFBGvRcRfSYmmLEnjge2Av1Q+9J4bOBAmTXKCMTPLM8HsADxdsr8kKytbJyLWAsuBkV1s/2hS\nj6X0dvqRkh6QdKWkHcudJGmGpEZJjc3NzV38qe75+Mdht918o9/M+rc8E4zKlLX9k9uVOu2ZBlxR\nsv9HoC4i9gRuoLVntGHjETMjoj4i6kePHt3Fn+qeb34T5s5N08eYmfVXeSaYJUBpL2IssLS9OpIG\nAcOBlzprWNJ7gUERcXdLWUS8GBFvZrsXAXv3PPTKWL266AjMzIqTZ4K5CxgvaWdJQ0g9jnlt6swD\njs22jwJuajPk1Z6j2bD3gqQxJbuHAQt7FHWFfPSjMG1akRGYmRUrtzf5I2KtpJOAa4GBwC8iYoGk\nM4HGiJgHXAJcKqmJ1HN560+ypMXA1sAQSYcDB0bEw9nhfwA+3uYn/0nSYcDarK3P53VtXTFmDNxw\nQ7oP46EyM+uP1LUOQ22qr6+PxsbGXNq+8EI46SRYvBjGjcvlJ8zMCiHp7oio76ye3+TPiV+4NLP+\nzgkmJ3vuCZtv7gRjZv2XZ1POyeDB8P3vw7vfXXQkZmbFcILJ0be/XXQEZmbF8RBZjtavh4cfhmXL\nio7EzKz3OcHk6Lnn0hDZnDlFR2Jm1vucYHK0/faw446+0W9m/ZMTTM4aGuD224uOwsys9znB5Gzy\nZHjySXj22aIjMTPrXU4wOWt54XL+/GLjMDPrbX5MOWd77QV//CPst1/RkZiZ9S4nmJxtthkcckjR\nUZiZ9T4PkfWCxx6Dc8+FNWuKjsTMrPc4wfSCu++G73wHHnqo6EjMzHqPE0wvaLnR78eVzaw/cYLp\nBePGwdvf7hcuzax/cYLpBVLqxTjBmFl/4gTTSxoa0guXy5cXHYmZWe/INcFIOkjSI5KaJJ1S5vhQ\nSXOz4/Ml1WXlIyXdLGmlpAvanHNL1uZ92We7jtrqK77ylZRchg8vOhIzs96RW4KRNBC4EDgYmAAc\nLWlCm2rHAS9HxG7AT4BzsvI3gNOBk9tpfnpETMw+z3fSVp+w1VbpnRgzs/4izx7MJKApIhZFxGpg\nDjC1TZ2pwKxs+0pgiiRFxGsR8VdSoumqsm31PPzKO/98+OY3i47CzKx35JlgdgCeLtlfkpWVrRMR\na4HlwMgutP3LbHjs9JIk0tO2es3ChXDxxWkhMjOzWpdnginXe4ge1GlrekTsAeyffT7XnbYkzZDU\nKKmxubm5k5+qrMmT4dVXU6IxM6t1eSaYJcCOJftjgaXt1ZE0CBgOvNRRoxHxTPa9AricNBTX5bYi\nYmZE1EdE/ejRo7t5SZum5YVLP65sZv1BngnmLmC8pJ0lDQGmAfPa1JkHHJttHwXcFBHt9mAkDZI0\nKtseDBwCtEzA0q22ijB+PIwY4QRjZv1DbrMpR8RaSScB1wIDgV9ExAJJZwKNETEPuAS4VFITqbcx\nreV8SYuBrYEhkg4HDgSeBK7NkstA4AbgouyUdtvqKyT4xCdg882LjsTMLH/qY//I71X19fXR2NhY\ndBhmZlVF0t0RUd9ZPb/JX5B+nNfNrJ9wgullq1fDhAlw1llFR2Jmli8nmF42ZEi6F+Op+82s1jnB\nFGDy5PQkmYfJzKyWOcEUoKEBXnoJmpqKjsTMLD9OMAXwCpdm1h84wRRgwgQ4/njYeeeiIzEzy09u\nL1pa+wYMgJkzi47CzCxf7sEUJCLdg3nzzaIjMTPLhxNMQa65Js1NNn9+0ZGYmeXDCaYgk7I5oD3x\npZnVKieYgowaBbvt5gRjZrXLCaZADQ3pUWW/cGlmtcgJpkANDfDcc/DUU0VHYmZWeX5MuUCHHgpj\nxsDIkUVHYmZWeU4wBdppp/QxM6tFHiIr2EMPwezZRUdhZlZ5TjAFu/RS+MIX/MKlmdUeJ5iCNTSk\nRcjuvbfoSMzMKivXBCPpIEmPSGqSdEqZ40Mlzc2Oz5dUl5WPlHSzpJWSLiipP0zS/0j6u6QFks4u\nOfZ5Sc2S7ss+X8zz2ipln33St2dWNrNak1uCkTQQuBA4GJgAHC1pQptqxwEvR8RuwE+Ac7LyN4DT\ngZPLNP3DiNgd2AvYV9LBJcfmRsTE7HNxBS8nN9tvn270+4VLM6s1efZgJgFNEbEoIlYDc4CpbepM\nBWZl21cCUyQpIl6LiL+SEs1bImJVRNycba8G7gHG5ngNvWLyZLjrrqKjMDOrrDwTzA7A0yX7S7Ky\nsnUiYi2wHOjSWyGStgEOBW4sKT5S0gOSrpS0YzvnzZDUKKmxubm5a1eSsx//GB58sOgozMwqq8ME\nI2nrDo519gaHypS1nRSlK3XK/fYg4Arg/IhYlBX/EaiLiD2BG2jtGW3YeMTMiKiPiPrRo0d39lO9\nYvvtYYstio7CzKyyOuvB3NKyIenGNsf+0Mm5S4DSXsRYYGl7dbKkMRx4qZN2AWYCj0XEeS0FEfFi\nRLQ87HsRsHcX2ukzzjgDLrqo6CjMzCqnswRT2sPYtoNj5dwFjJe0s6QhwDRgXps684Bjs+2jgJsi\nOp76UdK/khLR19qUjynZPQxY2El8fco118BllxUdhZlZ5XQ2VUy0s11uf8ODEWslnQRcCwwEfhER\nCySdCTRGxDzgEuBSSU2knsu0lvMlLQa2BoZIOhw4EHgVOA34O3CPJIALsifG/knSYcDarK3Pd3Jt\nfUpDQ1pGec0aGDy46GjMzDZdZwlmO0nfIPVWWrbJ9ju9gRERfwL+1Kbsn0u23wA+1c65de00W7bn\nFBGnAqd2FlNf1dAAP/1putn/vvcVHY2Z2abrbIjsImArYMuS7Zb9qnjPpFpMnpy+/T6MmdWKDnsw\nEXFGbwXS3+20E7zrXbByZdGRmJlVRocJRtLxwC0R8ZjSDY9LgCOBJ4FjI8IzaFWIBAsWpG8zs1rQ\n2RDZV4HF2fbRwHuBXYBvAOfnF1b/5ORiZrWkswSzNiLWZNuHAL/O3je5AfCrgRXW1AR77JEeWTYz\nq3adJZj1ksZI2gyYQnpDvsXm+YXVP40ZAwsXwt/+VnQkZmabrrPHlP8ZaCS9xzIvIhYASPogsKij\nE637ttgC9tzTU/ebWW3oLMEsAyYDKyLiZUnHkG7yLwNm5B1cf9TQkN7oX7cOBg4sOhozs57rbIjs\nv4CVWXL5AHA28GtSgvlp3sH1Rw0NsGJFGiozM6tmnfVgBkZEy+STnwZmRsTvgN9Jui/f0Pqn/faD\n6dP9RJmZVb9OE4ykQdlaLVPYcFiss3OtB3bZxZNemllt6CxJXAHcKukF4HXgLwCSdiMtDmY5iICl\nS2GHtsuzmZlVkQ7vwUTEWcA3gV8B+5VMpT8A+Md8Q+u/fvhD2HFHWO4UbmZVrNNhrojYaPrFiHg0\nn3AMYOLE1Iu580444ICiozEz65nOniKzAkyalG7ye2ZlM6tmTjB90PDhMGGCE4yZVTcnmD6qoSEl\nmI4XkDYz67v8qHEfddxx6f7LunUwyP8rmVkV8p+uPmry5NZVLs3MqlGuQ2SSDpL0iKQmSaeUOT5U\n0tzs+HxJdVn5SEk3S1op6YI25+wt6cHsnPOzhdCQtK2k6yU9ln2PyPPaesN998EttxQdhZlZz+SW\nYCQNBC4EDgYmAEdLmtCm2nHAyxGxG/AT4Jys/A3gdODkMk3/nDSjwPjsc1BWfgpwY0SMB27M9qva\nt74FX/960VGYmfVMnj2YSUBTRCyKiNXAHGBqmzpTgVnZ9pXAFEmKiNci4q+kRPMWSWOArSPi9uyl\nz18Dh5dpa1ZJedVqaIAHHoDXXis6EjOz7sszwewAPF2yvyQrK1snm+9sOTCykzaXtNPm2yLi2ayt\nZ4HtyjUgaYakRkmNzc3NXbyUYjQ0wPr10NhYdCRmZt2XZ4IpNx9w24duu1JnU+pvXDliZkTUR0T9\n6NGju3Nqr9tnn/Tt92HMrBrlmWCWADuW7I8FlrZXR9IgYDjwEu1bkrVTrs1l2RBay1Da8z2OvI8Y\nNQrGj/cKl2ZWnfJMMHcB4yXtLGkIMA2Y16bOPODYbPso4KaSCTU3kg19rZDUkD09dgxwdZm2ji0p\nr2p/+AP8+tdFR2Fm1n25JZjsnspJwLXAQuA3EbFA0pmSDsuqXQKMlNQEfIOSJ78kLQZ+DHxe0pKS\nJ9C+DFwMNAGPA9dk5WcDB0h6DDgg2696994Le+4JAwZAXR3Mnl10RGZmXaMOOgw1r76+Phr78B30\n2bPh+OPh9ddby4YNg5kz06qXZmZFkHR3RNR3Vs9zkfVhp522YXIBWLUqlZuZ9XVOMH3YU091r9zM\nrC9xgunDdtqpe+VmZn2JE0wfdtZZ6Z5LqYEDU7mZWV/nBNOHTZ+ebuiPG5dWuNxuOzjvPN/gN7Pq\n4On6+7jp0zdOKOvWwd/+BvvvX0xMZmZd4R5MFfrRj+BDH4I//7noSMzM2ucEU4VOPBH22AM+/Wn4\n+9+LjsbMrDwnmCq0xRZw9dUwdCgcdhi8/HLREZmZbcwJpkqNGwdXXQWLF8MxxxQdjZnZxnyTv4rt\nuy/MmgW77FJ0JGZmG3OCqXJHH926/fjjsOuuxcViZlbKQ2Q1YtYseNe74NZbi47EzCxxgqkRhx+e\nei9HHgmLFhUdjZmZE0zNGD4c5s2D9eth6lRYsaLoiMysv3OCqSHjx8NvfwsLF6a3//vxUj9m1gf4\nJn+NmTIFLrgANt88zV9mZlYUJ5ga9KUvtW4vX56Gz8zMepuHyGrYdddBXR3Mn190JGbWH+WaYCQd\nJOkRSU2STilzfKikudnx+ZLqSo6dmpU/IuljWdk7Jd1X8nlV0teyY9+X9EzJsY/neW3V4H3vgxEj\n0hNmS5YUHY2Z9Te5JRhJA4ELgYOBCcDRkia0qXYc8HJE7Ab8BDgnO3cCMA14N3AQ8DNJAyPikYiY\nGBETgb2BVcBVJe39pOV4RPwpr2urFqNGwR//CCtXpiSzalXREZlZf5JnD2YS0BQRiyJiNTAHmNqm\nzlRgVrZ9JTBFkrLyORHxZkQ8ATRl7ZWaAjweEU/mdgU14N3vhssvh3vugS98wU+WmVnvyTPB7AA8\nXbK/JCsrWyci1gLLgZFdPHcacEWbspMkPSDpF5JGlAtK0gxJjZIam5ubu3M9VevQQ+Hss2GHHdJ7\nMmZmvSHPBFPuIdm2/35ur06H50oaAhwG/Lbk+M+BXYGJwLPAj8oFFREzI6I+IupHjx7dfvQ15tvf\nTguVDRyYVsQ0M8tbnglmCbBjyf5YYGl7dSQNAoYDL3Xh3IOBeyJiWUtBRCyLiHURsR64iI2H1Ay4\n/36YMCF9m5nlKc8EcxcwXtLOWY9jGjCvTZ15wLHZ9lHATRERWfm07CmznYHxwJ0l5x1Nm+ExSWNK\ndo8AHqrYldSQ7baD115LC5U9/3zR0ZhZLcstwWT3VE4CrgUWAr+JiAWSzpR0WFbtEmCkpCbgG8Ap\n2bkLgN8ADwN/Bk6MiHUAkoYBBwC/b/OT50p6UNIDwIeBr+d1bdVszJi0GmZzc5oY8803i47IzGqV\noh8/VlRfXx+NjY1Fh1GIuXNh2rT0ZNnFF3taGTPrOkl3R0R9Z/U8VUw/9elPw0MPwV13pV7MZpsV\nHZGZ1RonmH7sjDPSY8uD/F+BmeXAc5H1YwMGpOTy7LNw0EHwyCNFR2RmtcQJxli9Or3pf+ih8PLL\nRUdjZrXCCcYYNw6uugoWL4b990/7AwakmZhnzy46OjOrVk4wBsC++8LnPw8LFsBTT6U5y558EmbM\ncJIxs55xgrG3XHfdxmWrVsFpp/V+LGZW/Zxg7C1PPdW9cjOzjjjB2Ft22qn9Y6ecAi+91HuxmFn1\nc4Kxt5x1FgwbtmHZZpvB5Mlw7rmwyy7wb//mNWXMrGucYOwt06fDzJnpKTIpfV98Mdx2W5p9+YMf\nhDvvbJ1WxmvLmFlHPBdZP52LrKfefBOGDoVHH4VPfAK+9z347GfTOjNm1j90dS4y92CsW4YOTd8r\nV8Lw4enR5j33TO/R9ON/q5hZGU4w1iPve1+aKPO3v00rZH7yk2kIzatlmlkLJxjrMQmOOirNynzJ\nJWk+s5ahsoULi43NzIrnBGObbNCgtK7Md7+b9m+5JS3LfMQRaWYAM+ufnGCs4vbeG848E266CfbY\nA449Ns1zZmb9ixOMVdxWW8Hpp8OiRXDyyfCb36S5ztasKToyM+tNuSYYSQdJekRSk6RTyhwfKmlu\ndny+pLqSY6dm5Y9I+lhJ+WJJD0q6T1JjSfm2kq6X9Fj2PSLPa7POjRyZXtBsaoJZs2Dw4PQQwA9/\nCK+8kibRrKvzzM1mtSq392AkDQQeBQ4AlgB3AUdHxMMldb4C7BkRX5I0DTgiIj4taQJwBTAJ2B64\nAXhHRKyTtBioj4gX2vzeucBLEXF2lsxGRMR3OorR78H0vltugQ9/GDbfHNau3bBXM2xYetFz+vTC\nwjOzLugL78FMApoiYlFErAbmAFPb1JkKzMq2rwSmSFJWPici3oyIJ4CmrL2OlLY1Czi8AtdgFfah\nD8G996Z3ZtoOmXnmZrPakudq7DsAT5fsLwH2aa9ORKyVtBwYmZXf0ebcHbLtAK6TFMB/RcTMrPxt\nEfFs1tazkrar5MVY5UycmGYEKOfJJ+GQQ+A970mfPfaA3XdvfcHTzKpHnglGZcrajse1V6ejc/eN\niKVZArle0t8j4n+7HJQ0A5gBsFNH0wdbrnbaKSWTtrbaKi0PcN11rT2cKVPghhvS9s9+BmPGpOSz\nyy6eosasL8tziGwJsGPJ/lhgaXt1JA0ChgMvdXRuRLR8Pw9cRevQ2TJJY7K2xgDPlwsqImZGRH1E\n1I8ePbrHF2ebptzMzcOGwc9/Dg88AK+9lt6hmTMHvv71dHzNmrT9yU/CO94BW26ZHomemfVhI+CZ\nZ9qfssYPFZj1rjx7MHcB4yXtDDwDTAM+06bOPOBY4HbgKOCmiAhJ84DLJf2YdJN/PHCnpC2AARGx\nIts+EDizTVtnZ99X53httolabuSfdlrqsey0U0o6LeWDB6eXNSdMaD1n8OC0Js3DD6fZAx58MH23\nJJSlS2HsWBgxYsMhto99DG6/PS3/vGpVqtuyHHRpLGZWYRGR2wf4OOlJsseB07KyM4HDsu3NgN+S\nbuLfCexScu5p2XmPAAdnZbsA92efBS1tZsdGAjcCj2Xf23YW39577x1WO158MeKCCyJOOCFi330j\nhg+PgIjLLosYNy5tt/2MGhWxbFk6/403Il5/vXLxtPyulL4vu6xybZsVCWiMLuQAT9fvx5RrVsuQ\n2dZbwzbbtD909pe/wH77wWWXwec+l4beRo2C0aPT9wUXpPs999+f1sNpKW/5HjEiDbuVmj17wx4T\n+DFsqx1dfUw5zyEys0JJacgM2n+oYMyYNDM0wHvfm4bpXngBmptbvwdl/1/y5z+npaPbWrIEdtgB\n/vM/4fLLU9K5/voNkwu0PobtBGP9hROM9QtnnVW+R/GDH7Q+bLDHHunTnq9+FT7zmZR4SpNQy7Mi\nQ4aknsyjj6b1csp58sl0P6i+Pt1TMqtlHiLzEFm/MXt2+w8VVFpdXfkeU4sttkjDclOnwpe/nE8M\nZnnpC2/ym/Up06enWZ3Xr0/feQ5VdfQY9m9/m1YCffppuPHG1uNf/GLqUTU2euE2qw0eIjPLQWeP\nYR91VPpueZl05Uq47ba0cBukBxM+8AE48cS0kJtZNXIPxiwnXekxtdyH2XLLtAro0qVwxRUwbVq6\nl7NsWTr+6KNpOO288+C++1KbpfwSqfVF7sGY9SFjxqTkMm1a2i99iXTBApg3L+1vuy188INp6QO/\nRGp9lW/y+ya/VZGnn05LHtx8M9x6a3ovZ++9yz9QsOOOaXjOrNK6epPfCcYJxqrcgAHtv0S6fn16\nH+g//iMlobq61s+4cWlyUbPu8ouWZv1Eey+RbrttSi4A8+fD734Hb7zRenz33dN9H4B/+Zc0wWhp\nAqqrS/eGyunNR76tejnBmFW59l4iPf/81v3LLoNLL4Xnn08PHCxe3Jp8IN3H+ctfNkxABx4I116b\ntk88Ma1CWleXzr3wwta6vudj7XGCMatynT0S3UKCt70tffZps/TfddelYbbSBDR8eDoWAX/9a3qS\nrTQBlVq1Kr3H86c/palyRo5MS2Pvv396p+f++1vLhw3bMLl1hXtM1ckJxqwGTJ++6X9w20tAUkoQ\nLQlozJjy93zeeAPuuCNNn/MZJya5AAAJ+0lEQVTqq6nO/vunc/beu7Xe0KEp0ZxxRkpKy5alIbqR\nI1uT0MiRaY64MWNS7+uEE/yUXDVygjGzLmlJQO3d8xk3Dh5/PG2vWdM6G8Hw4XDVVSnxvPhi66dl\nQdnmZvj979NaP6UzGMyaBcccAyefXH7i0JNPTgnmiSfgmmvSPaeRI9P3ttumCUiHDOn+dfZ2b6mW\ne2dOMGbWLe3d8znrrNb9wYNbXyIdNgwOP7z99t7zntTLWb8+9XxaElBdXTr+fNm1aVtfQm1sTPeI\n2mpZhuF3v4Pvfa818bQkoW99K/WQFi1Kn223Teeceiq8/npqI+/eUttlHWqtd+bHlP2Yslm39YWJ\nQ8eNS/eK1qxJvZ8XX0zfLZ9DDklDbjfdlJZSKK3z4otp2G/XXeHcc+E73+k4hoED03tFAwak7YED\n4d57YbPN0suuc+a0lg8cmHpOLfPMnXdeWr6h5diAAWkqoF/+sv1rGzky9fgg9dAi0ppGw4enNjZF\nJf6382PKZpabStzz6arOekyDB7feOyrnIx9Jn/YccwxMnpySzhFHlK+zbl2aG27duvRZv771D/1W\nW6Xfbjm2bt2GDzG8+mrqbbWct25d6/tH7b0I++KLrdvHH7/hpKjDh8OkSenBDIDvfheeey4loBEj\n0vduu8HBB6fjjz+eEuE228Af/tC7PSb3YNyDMevzeqvH1Flvqbd+b+zYNGsDpBkbnngCXnkFXn45\nfY8aBaefno4fcUQaJnzlldZ1iKZMgRtuSNu77pqGADvS3evzm/xd4ARjZqV6e6nrSv/emjWwfDms\nXQtvf3squ/rq1IN6+eXyK7JC6nG1nUC1I31iPRhJB0l6RFKTpI0uTdJQSXOz4/Ml1ZUcOzUrf0TS\nx7KyHSXdLGmhpAWSvlpS//uSnpF0X/b5eJ7XZma1Z/r09Md93Lj0R3fcuPySSx6/N3hw6t20JBdI\ns3DPmJHuM40bV/68lif6Ki23HoykgcCjwAHAEuAu4OiIeLikzleAPSPiS5KmAUdExKclTQCuACYB\n2wM3AO8AtgPGRMQ9krYC7gYOj4iHJX0fWBkRP+xqjO7BmFl/UqkeU1/owUwCmiJiUUSsBuYAU9vU\nmQrMyravBKZIUlY+JyLejIgngCZgUkQ8GxH3AETECmAhsEOO12BmVjN6u4eWZ4LZAXi6ZH8JGyeD\nt+pExFpgOTCyK+dmw2l7AfNLik+S9ICkX0gaUS4oSTMkNUpqbG5u7u41mZlVtd5cOjzPBFNutqG2\n43Ht1enwXElbAr8DvhYRr2bFPwd2BSYCzwI/KhdURMyMiPqIqB89enTHV2BmZj2WZ4JZAuxYsj8W\nWNpeHUmDgOHASx2dK2kwKbnMjojft1SIiGURsS4i1gMXkYbozMysIHkmmLuA8ZJ2ljQEmAbMa1Nn\nHnBstn0UcFOkpw7mAdOyp8x2BsYDd2b3Zy4BFkbEj0sbkjSmZPcI4KGKX5GZmXVZbm/yR8RaSScB\n1wIDgV9ExAJJZwKNETGPlCwuldRE6rlMy85dIOk3wMPAWuDEiFgnaT/gc8CDku7Lfuq7EfEn4FxJ\nE0lDaYuBE/K6NjMz65xftPRjymZm3eI3+btAUjNQZqKGPmkU8ELRQeSklq8Navv6fG3Va1Oub1xE\ndPqUVL9OMNVEUmNX/sVQjWr52qC2r8/XVr164/pynSrGzMz6LycYMzPLhRNM9ZhZdAA5quVrg9q+\nPl9b9cr9+nwPxszMcuEejJmZ5cIJxszMcuEE08d1tMharZA0UNK9kv676FgqSdI2kq6U9Pfsf7/J\nRcdUKZK+nv33+JCkKyRtVnRMmyKbgf15SQ+VlG0r6XpJj2XfZWdo7+vaubYfZP9dPiDpKknb5PHb\nTjB931rgmxHxLqABODFbkK2WfJW0tk+t+Snw54jYHXgvNXKNknYA/gmoj4j3kKaCmlZsVJvsV8BB\nbcpOAW6MiPHAjdl+NfoVG1/b9cB7ImJP0sKQp+bxw04wfVytL7ImaSzwCeDiomOpJElbAx8gzbdH\nRKyOiFeKjaqiBgGbZ7OgD2PjmdKrSkT8L2k+xFKlCyLOAg7v1aAqpNy1RcR12RpcAHeQZqyvOCeY\nKtLOImvV7jzg28D6ogOpsF2AZuCX2fDfxZK2KDqoSoiIZ4AfAk+R1l5aHhHXFRtVLt4WEc9C+oce\nacn2WvQF4Jo8GnaCqRLtLLJW1SQdAjwfEXcXHUsOBgHvA34eEXsBr1G9QywbyO5FTAV2BrYHtpD0\n2WKjsp6QdBppGH52Hu07wVSB9hZZqwH7AodJWgzMAT4i6bJiQ6qYJcCSiGjpbV5JSji14KPAExHR\nHBFrgN8D/6fgmPKwrGWdqez7+YLjqShJxwKHANMjpxcinWD6uI4WWat2EXFqRIyNiDrSTeKbIqIm\n/iUcEc8BT0t6Z1Y0hbS+US14CmiQNCz773MKNfIAQxulCyIeC1xdYCwVJekg4DvAYRGxKq/fcYLp\n+/YlLbL2EUn3ZZ+PFx2Udck/ArMlPQBMBP6t4HgqIuuVXQncAzxI+jtS1dOqSLoCuB14p6Qlko4D\nzgYOkPQYcEC2X3XaubYLgK2A67O/Kf+Zy297qhgzM8uDezBmZpYLJxgzM8uFE4yZmeXCCcbMzHLh\nBGNmZrlwgjHrAUl1pbPTdvPcX0k6qoe/+Zme/KZZEZxgzKpHHdCtBCNpYD6hmHXOCcZsE0naJZvQ\n8v1ljn1b0oOS7pe00Yt6khZLGpVt10u6Jdv+YMmLtfdK2or0ot/+WdnXs3V0fiDprmxdjxOycz+U\nrSF0OfCgpC0k/U8Ww0OSPp3n/z3MWgwqOgCzapZNBTMH+L8RcV+bYweTpnjfJyJWSdq2G02fDJwY\nEbdlE52+QZos8+SIOCRrfwZpJuP3SxoK3CapZVbjSaT1Pp6QdCSwNCI+kZ03vOdXbNZ17sGY9dxo\n0vxUn22bXDIfBX7ZMtdTRLRdb6QjtwE/lvRPwDYla3eUOhA4RtJ9pCUcRgLjs2N3RsQT2faDwEcl\nnSNp/4hY3o04zHrMCcas55YDT5PmiytHQGdzMa2l9f8P31p2OCLOBr4IbA7cIWn3dtr/x4iYmH12\nLlmX5bWSth4F9iYlmn+X9M+dxGRWEU4wZj23mjQEdkw7T3ddB3xB0jBIa7yXqbOY9Mcf4MiWQkm7\nRsSDEXEO0AjsDqwgTVDY4lrgy9lyDkh6R7lFzSRtD6yKiMtIC4XVyrIB1sf5HozZJoiI17KF066X\n9FpEXF1y7M+SJgKNklYDfwK+26aJM4BLJH2XDVcq/ZqkDwPrSNP8X0Na9XOtpPtJ66z/lPRk2T3Z\ntPnNlF/Wdw/gB5LWA2uAL2/iZZt1iWdTNjOzXHiIzMzMcuEEY2ZmuXCCMTOzXDjBmJlZLpxgzMws\nF04wZmaWCycYMzPLxf8H9mSydvgA/5cAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "SSE = []\n", "coords = list(cluster_df[['type_A_norm','longitude_norm','latitude_norm']].values)\n", "for k in range(1,13):\n", " km = KMeans(n_clusters=k)\n", " km.fit(coords)\n", " labels = km.labels_\n", " centroids = km.cluster_centers_\n", " # Get the SSE\n", " SSE.append(sum(np.min(cdist(coords, centroids, 'euclidean'),axis=1))/len(coords))\n", " \n", "plt.plot(range(1,13),SSE, 'bo--')\n", "plt.xlabel('k clusters')\n", "plt.ylabel('SSE')\n", "plt.title('SSE vs k')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4 seems like a decent choice, but interestingly so is 7. Let's try 7." ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAHiCAYAAAD8hSV1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzsvX+cHFWZ7/8+3TPToQkE6JCBGKfa\nuUqIiCJEIVcur+ioCyirq16FbRAC2JuJXsN63b177b0L6Pb+8LoLcTVgrxtUur7hu351V1GCwjAJ\nSMbF+Iu4hgEcp4cQpyETCcKYyWT6fP+o7kz/qOqu7q7+NfO8X696zfSp6lOnqrvrc85znvM8SmuN\nIAiCIAidi6/VDRAEQRAEoT5EzAVBEAShwxExFwRBEIQOR8RcEARBEDocEXNBEARB6HBEzAVBEASh\nwxExF4QGopS6RSmVbHU7qkUpdZ1S6getbocTSqlPKaW+3Op2CEK7IGIuCHWilPpjpdQepdRLSqnf\nKKV2KKUu9rD+sFJKK6W6vKqzHprdHqXUeqXU/vwyrfXfaK1vbMb5BaETEDEXhDpQSn0CuB34G6AX\n6AO2Au9pZbvyaZdOgCAIjUPEXBBqRCm1DPg08FGt9Te11i9rrWe11vdqrf/M5viSEaZSalwp9fbs\n/2/OjvBfVEqllVL/mD3s4ezfF7Kj/3XZ469XSu1TSv1WKfU9pZSRV69WSn1UKfUU8FS27Gyl1ANK\nqUNKqVGl1Afzjg8ppb6dPfdjwH+p8Z4ElFK3K6UOZLfblVKBvP3vUUr9LHueXymlLs2Wb8hey++U\nUmNKqT/Jlp8I7ABWZq/9JaXUyuLpC6XUHyql/lMp9YJSaqdSak3RPf6kUupxpdRhpdT/q5Rakt23\nXCn1nez7DimlHlFKyXNR6DjkSysItbMOWAL8m0f1bQG2aK1PxhLTf82WX5L9e4rWeqnWekQp9V7g\nU8D7gNOBR4DtRfW9F7gQeG1WFB8A/h9gBXAVsFUpdU722C8CR4AzgeuzWy3EgIuA84A3AG8G/hKs\nzgrwNeDPgFOy1zWefd9zwLuBk4ENwG1KqfO11i8DlwEHste+VGt9IP+ESqmzstd+U/Ze3Afcq5Tq\nyTvsg8ClwKuA1wPXZcv/J7A/+75erHsqMa6FjkPEXBBqJwQc1Fof86i+WeDVSqnlWuuXtNY/LHPs\nnwB/q7Xelz3/3wDn5Y/Os/sPaa1/jyWU41rru7TWx7TWPwG+AXxAKeUH3g/8Vda68AvgqzVeQwT4\ntNb6Oa3188CtwDXZfTcA27TWD2itM1rrZ7XWTwBorb+rtf6VttgFfB/4by7P+SHgu9l6Z4HPAScA\n/zXvmM9rrQ9orQ8B92J1NsC652cCRtaq8oiWhBVCByJiLgi1MwUs93BO+gbgLOAJpdSPlFLvLnOs\nAWzJmodfAA4BCnhF3jHPFB1/Ye747HsiwBlYo9KuouNTNV7DyqL3prJlAK8EfmX3JqXUZUqpH2ZN\n3S8AlwPLazmn1jqDdS3592Iy7/9pYGn2//8LPA18P2ve/wuX5xSEtkLEXBBqZwTLNP1el8e/DARz\nL7Ij4tNzr7XWT2mtr8Iyg/898P9lzeN2I8VngD/RWp+St52gtd6dd4wuOn5X0fFLtdaDwPPAMSyx\nzdHn8pqKOYDVccivJ2cWfwabufjsnPo3sEbUvVrrU7BM5crmOiqeUymlsK7l2UqN1Vr/Tmv9P7XW\n/cAVwCeUUgOV3icI7YaIuSDUiNb6MPBXwBeVUu9VSgWVUt3ZUeZnbd7yJLBEKfUupVQ31lxyvnPY\n1Uqp07MjyxeyxXNYYpsB+vPquhP437k5b6XUMqXUfy/T3O8AZymlrsm2sVsp9Sal1Bqt9RzwTeCW\n7DW8FrjWxS0IKKWW5G0+rLnrv1RKna6UWp69PzlHtX8BNiilBpRSPqXUK5RSZwM92fvwPHBMKXUZ\n8M6886SBkLIcDu34V+Bd2Xq7sebBZ4DdDscfRyn1bqXUq7MdgBex7veci2sXhLZCxFwQ6kBr/Y/A\nJ7CE+Xms0efHgH+3OfYwsAn4Mtao8WUs56sclwL/qZR6CcsZ7kqt9RGt9TQQBx7Nmsgv0lr/G9bo\n/R6l1IvAL7AcxZza+TssgbwSayQ7mX1/rjPxMSzT8yTwFeAuF5f/EvD7vO1twF8De4DHgb3AT7Jl\naK0fI+vcBhwGdmHNVf8O+DiWKP8W+GPg23ltfwKrkzCWvf6c2T63fxS4Gvgn4CDWCPsKrfVRF9fw\nGuDB7LWMAFu11jtdvE8Q2golvh6CIAiC0NnIyFwQBEEQOhwRc0EQBEHocETMBUEQBKHDETEXBEEQ\nhA5HxFwQBEEQOpyOyqa0fPlyHQ6HW90MQRAEQWgKP/7xjw9qrU+vdFxHiXk4HGbPnj2tboYgCIIg\nNAWllKvQymJmFwRBEIQOR8RcEARBEDocEXNBEARB6HA6as5cEARBEKpldnaW/fv3c+TIkVY3xZEl\nS5awatUquru7a3q/iLkgCIKwoNm/fz8nnXQS4XAYK0Fee6G1Zmpqiv379/OqV72qpjrEzC4IgiAs\naI4cOUIoFGpLIQdQShEKheqyHIiYC4IgCAuedhXyHPW2T8RcEARBEBrM/fffz+rVq3n1q1/N3/3d\n33levysxV0ptU0o9p5T6hcN+pZT6vFLqaaXU40qp8/P2XauUeiq7XZtXfpVSam/2+PuVUsvrvxxB\nEARBqBPThHAYfD7rr2nWVd3c3Bwf/ehH2bFjB7/85S/Zvn07v/zlLz1pag63I/OvAJeW2X8Z8Jrs\nFgXuAFBKnQbcDFwIvBm4WSl1qlKqC9gCvFVr/XrgceBjtVyAIAiCIHiGaUI0CqkUaG39jUbrEvTH\nHnuMV7/61fT399PT08OVV17Jt771LQ8b7VLMtdYPA4fKHPIe4Gva4ofAKUqpM4E/AB7QWh/SWv8W\neACrU6Cy24nKmig4GThQx3UIgiAIQv3EYjA9XVg2PW2V18izzz7LK1/5yuOvV61axbPPPltzfXZ4\ntTTtFcAzea/3Z8tsy7XWs0qpQWAv8DLwFPBRj9oiCIIgCLUxMVFduQu01iVlXjvkeeUAZ9cq7VSu\nlOoGBoE3AiuxzOz/27ZipaJKqT1KqT3PP/+8R80VBEEQBBv6+qord8GqVat45pn5ce3+/ftZuXJl\nzfXZ4ZWY7wdemfd6FZbZ3Kn8PACt9a+01WX5V+C/2lWstU5orddqrdeefnrFLHCCIAiCUDvxOASD\nhWXBoFVeI29605t46qmn+PWvf83Ro0e55557+MM//MM6G1qIV2L+beDDWa/2i4DDWuvfAN8D3pl1\nejsVeGe27FngtUqpnDq/A9jnUVsEQRAEoTYiEUgkwDBAKetvImGV10hXVxdf+MIX+IM/+APWrFnD\nBz/4Qc455xwPG+1yzlwptR1YDyxXSu3H8lDvBtBa3wncB1wOPA1MAxuy+w4ppT4D/Chb1ae11oey\ndd4KPKyUmgVSwHXeXJIgLAxMEzZvhqkp63UoBFu21PVMEQTBDZGI5z+0yy+/nMsvv9zTOvNxJeZa\n66sq7Nc4OLBprbcB22zK7wTudHN+QVhMmOk0m/eOMXXmDHwhAF/uh6FepqZgwwbrGBF0QRDykQhw\ngtBGmOk00dFRprpmrF/nGTPwyVEYSAMwO1vXChlBEBYoIuaC0EbExsaYzmQKC5dk4C/2HRf0OlbI\nCIKwQJEUqILQRkzMzNjv6MIaoQN9T/c2r0GCIHQEIuaC0Eb0BQKknAR9SQY+Mkb8mIi5IAiFiJld\nENqIeH8/QV+Zn+WKGXF+EwShBBFzQWgjIr29JFavxu+w31gSaGp7BEHwhuuvv54VK1bwute9riH1\ni5gLQpsR6e3lq2vWlIzQFXB5KNSaRgnCIsLcaxK+PYzvVh/h28OYe+tLgQpw3XXXcf/993vQOntE\nzAWhDYn09nLtGWcUJDfQwFcnJzHT6VY1SxAWPOZek+i9UVKHU2g0qcMpovdG6xb0Sy65hNNOO82j\nVpYiYi4Ibcp9U1MU51qazmSIjY21pD2CsBiIDcWYni1MgTo9O01sqL0DPIiYC0Kb4rRMzXH5miAI\ndTNx2D6Qg1N5uyBiLghtSl/A3tnNqVwQhPrpW2af6tSpvF0QMReENsVumVrQ5yPe39+iFgnCwic+\nECfYXZgCNdgdJD5QewrUZiBiLghtSm6ZmhEIoAAjECCxejWR3sYHjTFNCIfB57P+mvU78wpCRxA5\nN0LiigTGMgOFwlhmkLgiQeTc+gI8XHXVVaxbt47R0VFWrVrFv/zLv3jUYgtlJTzrDNauXav37NnT\n6mYIwoLGNCEahek8H6BgsO6Uzm2PudckNhRj4vAEfcv6iA/E636AC+3Bvn37WLNmTaubURG7diql\nfqy1XlvpvTIyF4QOIZ02GRkJs3Onj5GRMOl0Y4bLsVihkIP1eiFna2vUciRBaBYi5oLQAaTTJk88\ncT0zMylAMzOT4oknrm+IoDtlZVvI2do6dTmSIOQQMReEDuCppzaj9dGCMq2P8tRTmz0/V5+D065T\n+UKgU5cjCUIOEXNB6ACOHZuqqjxHLY5s8bg1R55PMGiVL1Q6dTmSIOQQMReEBUrOkS2VAq2tv9Fo\nZUGPRCxnN8MApay/C935rVOXIwlCDhFzQehwnEbf9TiyRSIwPg6ZjPV3IQs5zC9Hes8rTmT7hTB0\nCdx1wRFOOPJoq5smCK4QMReEDsDvt8+WNjcXchx9L0ZHtno44cijbHzVy5yxBHwKVgQyBF+8g2/+\naFOrmyYsAJ555hne+ta3smbNGs455xy2bNniaf0i5oLQAZx11hagu6i0m3/+5y2Oo+/F6MhWD5nf\nJlhSlEh+id8qFxYXjQia1NXVxT/8wz+wb98+fvjDH/LFL36RX/7yl/VXnEXEXBA6gN7eCGvW3EVX\n1/wI3e8/mYMH7Y+fmFicjmz1cFr3XFXlwsKkVl+TSpx55pmcf/75AJx00kmsWbOGZ5991oMWW4iY\nC0IHkcn8/vj/c3NT/NmfRRkYKH3K9PUtTke2ejg066+qXFiYNCNo0vj4OD/96U+58MILPatTxFwQ\nOoSxsRiZTOFTJhCY5iMfKXzK5I++F5sjWz34To1ypGgQfmTOKhcWD432NXnppZd4//vfz+23387J\nJ5/sTaWImAtCxzAzY/80WbFiQkbfHvC+N21l+uRBDh71k9Fw8Kif6ZMHed+btra6aUITaaSvyezs\nLO9///uJRCK8733vq7/CPLo8rU0QhIYRCPRlw7kWsmRJH+PjzW/PQsQSbhHvxUw8bp9oqF5fE601\nN9xwA2vWrOETn/hEfZXZICNzQegQ+vvj+HyFHm0+X5D+fvFoEwSvaJSvyaOPPsrdd9/NQw89xHnn\nncd5553Hfffd502jkZG5IHQMvb3W02RsLMbMzASBQB/9/fHj5YIgeEMk4v1U1cUXX0wjU46LmAtC\nB9HbGxHxFgShBDGzC4IgCEKHI2IuCIIgCB2OiLkgCIIgdDgi5oIgeIJpmoTDYXw+H+FwGNOLgNaC\nILhCHOAEQagb0zSJRqNMZxfnplIpolErclpEItgIQsORkblQFhltLV5M02T58uUopVBKsXz5csfP\nPxaLHRfyHNPT08S8DGgtCB3MkSNHePOb38wb3vAGzjnnHG6++WZP65eRueCIjLYWL6ZpsmHDBmZn\nZ4+XTU1Ncd111wGln/+EQ+Bqp3JBaGfSadPzeA6BQICHHnqIpUuXMjs7y8UXX8xll13GRRdd5Emb\nZWQuOCKjrcVLLBYrEPIcx44dY/PmzSXlfQ6Bq53KBaFdSadNRkej2dDJmpmZFKOjUdLp+qySSimW\nLl0KWDHaZ2dnUUp50GILEXPBERltLV7KfcZTU1MlZfF4nGBR8vRgMEhckqcLHYZddsJMZpqxsfoH\nMXNzc5x33nmsWLGCd7zjHZICVWgOTqOqV7zitCa3RLDDNCEcBp/P+uulO0O1I+pIJEIikcAwDJRS\nGIZBIpGQ6Rih43DKTuhUXg1+v5+f/exn7N+/n8cee4xf/OIXddeZQ8RccCQej3PCCT0FZYEAXHfd\ni3WbnIT6ME0rs1MqBVpbf6NR7wS93Ig6FArZlkciEcbHx8lkMoyPj4uQCx1JIGDfkXUqr4VTTjmF\n9evXc//993tWp4i54EgkEuHP//wkenut7EG9vfDJT8LAwKwnJiehdmKxwhSNYL32yp0hEokwODhY\nUt7d3c2WLVu8OYkgtCGNyk74/PPP88ILLwDw+9//ngcffJCzzz67rjrzEW92oSzr1x9i/frSci9M\nTkLtOE1pe+nOsHXrVt7ylrcQi8WYmJigr6+PeDy+aEbcjfBoFtqfRmUn/M1vfsO1117L3NwcmUyG\nD37wg7z73e/2osmAiLlQgUCgL+vVWVoutI6+Psu0blfuJZFIZNGIdz45j+acI1TOoxkQQV8ENCI7\n4etf/3p++tOfelpnPhXN7EqpbUqp55RStjP1yuLzSqmnlVKPK6XOz9t3rVLqqex2bV55j1IqoZR6\nUin1hFLq/d5cjuA1jTI5CfURj0OR8zjBoFUu1E8jPZoFoRG4mTP/CnBpmf2XAa/JblHgDgCl1GnA\nzcCFwJuBm5VSp2bfEwOe01qfBbwW2FVL44XG09sbYfXqBIGAASgCAYPVqxMyOmkxkQgkEmAYlj+D\nYVivF+EguiE00qNZEBpBRTO71vphpVS4zCHvAb6mtdbAD5VSpyilzgTWAw9orQ8BKKUewOoUbAeu\nB87O1p8BDtZxDUKDaYTJSaifSETEu1HI9JLQaXjhzf4K4Jm81/uzZbblSqlTsq8/o5T6iVLq60qp\nXg/aIQiLkkauN1+syPSS0Gl4IeZ28eh0mfIuYBXwqNb6fGAE+Jxj5UpFlVJ7lFJ7nn/+eQ+aKwjN\nwUynCY+M4Nu5k/DICGY67f05GrzefLEi00tCp+GFN/t+4JV5r1cBB7Ll64vKdwJTwDTwb9nyrwM3\nOFWutU4ACYC1a9dqD9orCA3HTKeJjo4ynckAkJqZITo6CkCk1ztDVLn15mKCrw+ZXhI6CS9G5t8G\nPpz1ar8IOKy1/g3wPeCdSqlTs45v7wS+l51bv5d5oR8AfulBOwShbYiNjR0X8hzTmQyxsTFPz9OM\n9eaCIHjD3Nwcb3zjGz1dX56j4shcKbUdS3iXK6X2Y3modwNore8E7gMuB57GGnFvyO47pJT6DPCj\nbFWfzjnDAf8LuFspdTvwfO49grBQmJiZqaq8Vk47DWzynnCahM8XhJox02liY2NMzMzQFwgQ7+/3\nxKK2ZcsW1qxZw4svvuhBKwtx481+VYX9Gviow75twDab8hRwics2CkLH0RcIkLIR7r5AoAWtEQTB\nLY2aItu/fz/f/e53icVi/OM//qMnbc1HYrMLQgOI9/cT9BX+vII+H/H+fk/Pc+hQdeWCIJSnUVNk\nN910E5/97Gfx+RojuyLmgtAAIr29JFavxggEUIARCJBYvdpT5zdwDt/qdVhXQVgsNGKK7Dvf+Q4r\nVqzgggsuqLmOSkhsdkFoEJHeXs/Fu5h43FqKlu/RLmFdBaF2GjFF9uijj/Ltb3+b++67jyNHjvDi\niy9y9dVXk0wm62lqATIyF4QORsK6Ng9zr0n49jC+W32Ebw9j7pXF/AuRRkyR/e3f/i379+9nfHyc\ne+65h7e97W2eCjmImAtCxxOJwPg4ZDLWXxFy7zH3mkTvjZI6nEKjSR1OEb036rmgSzS/1tOsKTKv\nUZYzemewdu1avWfPnlY3QxCERUb49jCpw6Wx2o1lBuM3jXtyjlw0v+IpE7G01M++fftYs2ZNq5tR\nEbt2KqV+rLVeW+m9MjIXBEGowMRh+yg8TuW1UC6anyBUQsRcEAShAn3L7JcHOJVD9XPsEs1PqAcR\nc0EQhArEB+IEuwuzqAW7g8QH7JcNbPruJq755jVVzbHLMkOhHkTMBUEQKhA5N0LiigTGMgOFwlhm\nkLgiQeTc0slsc6/JnXvuRFPojzQ9O01syNlmHo9bc+T5yDJDwS2yzlwQBMEFkXMjtuJdTGwoViLk\nOcrNseec3GIxy7Te12cJuTi/CW4QMRcEQfCQcoJdbo4dLOEW8RZqQcRcEATBQ/qW9dkuY1Moxzl2\nYXEQDoc56aST8Pv9dHV14eVSa5kzFwRB8BA7ZzmFYuPaja7M9ELrSZtpRsIj7PTtZCQ8QtpMe1b3\n8PAwP/vZzzwVcpCRuSAIgqfkBDs2FGPi8AR9y/qID8RFyDuEtJlmNDpKZtrKnDaTmmE0aqVA7Y20\nbxQ4GZkLQpsiscA7l8i5EcZvGidzc4bxm8ZFyDuIsdjYcSHPkZnOMBarLwUqgFKKd77znVxwwQUk\nEom668tHRuaC0IbkYoFPz1ohwXLrlAERBkFoIDMT9qlOncqr4dFHH2XlypU899xzvOMd7+Dss8/m\nkksuqbtekJG5ILQlsaHYcSHPUWmdcrtjmibhcBifz0c4HMaULCJCGxLos0916lReDStXrgRgxYoV\n/NEf/RGPPfZY3XXmEDEXhDakGbHAm4lpmkSjUVKpFFprUqkU0WhUBF1oO/rj/fiChdLoC/roj9ee\nAhXg5Zdf5ne/+93x/7///e/zute9rq468xExF4Q2pJZY4O1MLBZjuiiLyPT0NDHJIiK0Gb2RXlYn\nVhMwAqAgYARYnVhdt/NbOp3m4osv5g1veANvfvObede73sWll17qUatlzlwQ2pL4QLxgzhzKxwJv\ndyYcsoU4lQtCK+mN9Hruud7f38/Pf/5zT+vMR0bmgtCGVBML3A4znSY8MoJv507CIyOYae/WydZC\nn0O2EKdyQRCqQ0bmgtCmuI0FXoyZThMdHWU6Yy2vSc3MEB211slGeluzTjYejxONRgtM7cFgkLhk\nEREET5CRuSAsMGJjY8eFPMd0JkNsrP51srUSiURIJBIYhoFSCsMwSCQSRCQQuSB4gozMBWGBMTFj\nvx42NXOkyS0p5E1vepSvfGU/oIH9rFz5KCBiLgheICNzQVhg9AUc1sNqzaYffbO5jcny5JObOHDg\nDmAuWzLHgQN38OSTm1rSHkFYaIiYC0Kbk06bjIyE2bnTx8hImHS6/NrseH8/aJt82spH4reZ0vIm\ncOCAfehKp3JBEKpDxFwQ2ph02mR0NMrMTArQzMykGB2NlhV0y8nNRsyBue7TGtPQisxVWS4IC48X\nXniBD3zgA5x99tmsWbOGkZERz+oWMReENmZsLEYmUxhsJZOZZmysfLAV/+yhqsobj7/KckFoHY0K\nPbx582YuvfRSnnjiCX7+85+zZs0aT+oFEXNBaCqmabJ8eRilfCgVZvlyk3LPiZkZ+6AqTuU5oqf6\nYK7I4W3uiFXeAlaujFZVLgitolGhh1988UUefvhhbrjhBgB6eno45ZRTvGgyIGIuCE3DNE2uvz7K\n1JRlMocUU1NRNmxwFvRAwD6oilN5jq1veh+DJ0/jP3oQdAb/0YMMnjzN1je9r76LqJGzztrKypWD\nzI/E/axcOchZZ21tSXsEwYlGhR4eGxvj9NNPZ8OGDbzxjW/kxhtv5OWXX66rznyUtnOUaVPWrl2r\n9+zZ0+pmCEJNhMNhUqmUzR4DwxhnfLx0T27OPN/U7vMFWb06QW+vLOsSBDfs27fPtUnb5/Nhp4tK\nKTKZ2h1I9+zZw0UXXcSjjz7KhRdeyObNmzn55JP5zGc+U7adSqkfa63XVmx3zS0TBKEqnOOQT+C0\nq7c3wurVCQIBA1AEAoZrIW+3kK6C0Ak0KvTwqlWrWLVqFRdeeCEAH/jAB/jJT35SV535iJgLQpNw\nfhj0Ue450dsbYd26cdavz7Bu3bhrIY+OjpKambEM+tmQriLoglCeeDxOMBgsKPMi9PAZZ5zBK1/5\nSkazoZWHhoZ47WtfW1ed+YiYC0KTiMfj9PQEi0qDdHfH8TpEeTuGdBWETqCRoYf/6Z/+iUgkwutf\n/3p+9rOf8alPfcqDFltIOFdBaBK5h8HmzTGmpiaAPkKhOFu2RPA6RLlTSFenckEQ5olEIg3JG3De\neefRKL8vEXNBaCKNekgU0xcIkLIRbsdQr4IgdDRiZheEBUi8v5+gr/DnHfT5rFCvgiAsOETMBWEB\nEuntJbF6NUYggAKMQIDE6tUty2cuCEJjETO7ICxQIr29It6CkEVrjVKq1c1wpN6YLzIyFwRBEBY0\nS5YsYWpqqm7BbBRaa6ampliyZEnNdcjIXBAEQVjQrFq1iv379/P888+3uimOLFmyhFWrVtX8fhFz\nQRCajrnXJDYUY+LwBH3L+ogPxImcu/DC05rpNLGxMSZmZugLBIj398vURwvo7u7mVa96Vaub0VBE\nzAVBaCrmXpPovVGmZ61486nDKaL3WtnTFpKg56Lw5YL35KLwASLogue4mjNXSm1TSj2nlPqFw36l\nlPq8UupppdTjSqnz8/Zdq5R6Krtda/PebzvVKwjCwiM2FDsu5DmmZ6eJDdWXlaqV2MXBlyh8QjNx\n6wD3FeDSMvsvA16T3aLAHQBKqdOAm4ELgTcDNyulTs29SSn1PuClqlstCELHMnHYPquMU3m74xQH\n3y5oD0gUPqExuBJzrfXDwKEyh7wH+Jq2+CFwilLqTOAPgAe01oe01r8FHiDbKVBKLQU+Afx1PRcg\nCEJn0bfMISuVQ3m74zQC9zscL1H4hEbg1dK0VwDP5L3eny1zKgf4DPAPQKG9TRCEBU18IE6wuygr\nVXeQ+IDH2WbqwDQhHAafz/prms7HOo2050Ci8AlNwysxt1uJr53KlVLnAa/WWv9bxYqViiql9iil\n9rTzsgJBECpjmhC7IsL0X76Ef8sz8PgfYywzSFyRaBvnN9OEaBRSKdDa+huNOgu600g7F3VPovAJ\nzUC5XUSvlAoD39Fav85m35eAnVrr7dnXo8D63Ka1/pP844BTgP8DHMXyqF8B7NZary/XhrVr1+pG\nZZwRBKGx5ERyOs8WFwxCIoHnWePqIRy2BLwYw4Dx8dLyYq91sEbgItyCFyilfqy1XlvpOK9G5t8G\nPpz1ar8IOKy1/g3wPeCdSqlTs45v7wS+p7W+Q2u9UmsdBi4Gnqwk5ELnUI2JUlg8xGKFQg7W61ib\nObFPOPjhOZVLHHyhHXC1zlyePfYbAAAgAElEQVQptR1rlL1cKbUfy0O9G0BrfSdwH3A58DTWHPiG\n7L5DSqnPAD/KVvVprXU5RzqhwykefeVMlNBeoy/Be0zTEuaJCejrg3i88DOvViRbRV+f/ci8r4x/\nnsTBF1qNazN7OyBm9vanWhNlO1FJjARn3JjQO+W7YZqwYQPMzs6XdXfDXXfJ90FoPs02swsC0Dmj\nr2KqdXoSCnFjQo/HLYHPJxi0ytuN4uRabpJtmXtNwreH8d3qI3x7GHOvfHmE5iFiLniKkymynImy\nHfByPncx+gy46cRFItZI3TAscTSM9nN+A+szP3q0sOzo0fLfhVyI2tThFBpN6nCKa755DZu+u6mx\njRWqZqH+PsXMLnhKp3gsF+PzWSPyYpSConggZenU66+XTjGhu6GW70L49jCpw6U3QKG4+313t82y\nu8VOJ/4+xcwutIROGX0V45VFoVM8tr2mk0zolajlu+AUilajOzrm/EJjIf8+RczbENM0CYfD+Hw+\nwuEwZofZgSIRazSWyVh/213IwTsx6lSfgXrp1E6cHbV8F8qFou3UmPMLkYX8+xQxbzNM0yQajZJK\npdBak0qliEajHSfonYZXYtQJPgONmjPsxE6cHbnvQig0X3bCCeXfEx+Io2wDXnZuzPmFSCf8PmtF\nxLzNiMViTBfZgaanp4ktBDtQm+OFGDXT3GzuNVn+2eWoWxXqVsXyzy6v6EHdTl777e6I9Pvfz/8/\nNVX+PkXOjbBx7cYSQW+3mPOLnYU0HVSC1rpjtgsuuEAvdJRSGiuufcGmlGp10wSXJJNaG4bWSll/\nk8kGnOPxpO7+dLfmFgq2ns/06OTjzic0DK0tGS/cDMP7NpYjmdQ6GCxsQzDYmHtVC7Xep+TjSW3c\nZmh1i9LGbUbZz0JoDfX+PpOTk9rYvVur4WFt7N6tk5OTjWjmcYA92oU+ijd7mxEOh0nZuAUbhsF4\np7kFCw3DyXsawFhmMH7TuO0+r7z266Xdvd/b5T4J7UUr4vCLN3uHEo/HCRbZgYLBIPEFYQcSvKKc\nU1W5fe0yZ9jujkjtcp+E1mOm04RHRvDt3Mm1+/bZ5q6PjY21qHXziJi3GZFIhEQigWEYKKUwDINE\nIkGkU72JhIZQzqmq3L52mTNsd7Fsl/sktJbcSDw1M4PGylFvh1NO+2YiYt6GRCIRxsfHyWQyjI+P\ni5A3iXZ3yMonPhCn29ddUt7j7ynrcFWL1/6mTdY9UcraTjqp/nvT7mK5kJbaCeXJH3mHR0Yw0+nj\n+2JjYyUjcTuccto3ExFzQaC9vLzdEDk3wl3vvYvQCfPrp0InhNj2nm0Vo41V47W/aRPccUfh/PFL\nL8F119V3bzpBLOtZ3bBp0ya6urpQStHV1cWmTRLWtR0pHnmnZmaIjo4eF3Q3I+6gz0e8v7/BLa2M\nOMAJAu3vkNUqurpgzsG2uJjuTaWMeqZpEovFmJiYIBgM8vLLL5fUMTg4yNatW5vYaqES4ZERUjaC\nbQQCjK9b57jfD2SwRuTx/v6Gpr916wAnYi4IiPeyE+WyhS2We1Mpnncu0FNxfIhi/H4/x44da3Br\nhWrw7dyJnQIqILN+fUu810vaIt7sguCednfIahV+v/O+RtybdvRbqBTP2y7Qkx1zTiYOoWU4zXXn\nyiO9vSRWr8YIBFBYI/ZmCnk1iJgLAu4csjo9Zn4tRKP25V1d3jir5Yv38uVw/fXt57dQaRndhMv1\ndP5yPSOhJcT7+wn6CmWweA480tvL+Lp1ZNavZ3zdurYUckAiwAlCjnKRoZLJpA4GgwVR+YLBoE62\nS8iyBjI4aN2TXBS0pUu9idRmFwWuHaLTFVMpGpxhGLZRG4u3wcHBVl6G4ECzI7pVCy4jwLVcoKvZ\nRMyFVuH0wDZarTR10Iyws+VwEsnirdGRjCvdh0qhZ5PJpGMYZkD7/X4RcqFmRMwFwUOaGTO/KbHd\nmxAbvdJ15I/2WzUyL7Y6ON2HStcyODhY8h1ZLJabxUD+5/+h0KT+ztLd+iGG9XZ26w+FJhvaERYx\nFwQPadbIvFkJSBqdcMXNdbgZmTcy+Uoy6dyhqOU+JJNJbRiGVkppwzBEyBcIyaTWPT3W92KASb2D\nXXqY4ePbjsCwHvjYw3rwG40xz7sVc1maJggusFt+FAwGPQ+126z17sqnQdusO1ManSmzHs0lbq7D\nbslXdzecfDIcOmS/nttLnNoIi2fZnWBhptPExsaYmJkpWTu+fLmVAhdgOyOcQem688le+OOv+rj7\nPO893WVpmnCcxeiF7TWNiplfvBTLSVy8TkDiP+XZqspzbNpkebIrZf11CmzmJpGKXRS4u+6Cgwfr\nyynvlnL3dLEvSVxMbHrySa7Zt88xClxOyAFW2Ag5wIrnQAdam3BFxHyBkxtRplIptNakUimi0agI\neg14HTPfLoSsU5AWr8Vl7q3/C7qLopR1v2yVO5AL7ZpbLj03Z722E3S36/brCZlaL6ed5ryvXWLE\nC43FTKe588CBksAxTpnQnsN+XfpzK6y/rUy4ImK+wLELaDE9PU0sF/FCaBl2wUi0LhX0RiQgMf7b\no3DFR2DZOJCx/l7xEavcgUTCfXm7J1Ipx9Kl7RUjXmgcsbEx2whwMC/Mofn0B3yZfo4UyeaRAHz5\nRut/PRloWbAjEfMFjlNAC7eBLgR31BK5zOkj0LrxCUjiA3GCF3wL/vRVcIsf/vRVBC/4VtmMa04B\nzOzKOyGRyqFD9uU2YdVrxjStOddcxrnly1sfBEeYp9xIOhcFbssWy5eDgTRD28f43P/OMHk6ZJQ1\nV/65T8LQ24EjPvhyf8uCHYkD3AInHA6TspmINQyD8cWSJcNjTBM2b56fSzvxRDh6FGZn54/Jj93t\nRKuTu5h7TWJDMSYOT9C3rI/4QLxsxjWnpCt+P3gZcrxSUhOvaPT9N00rot3Ro4Xl3d2Wb0A7dWwW\nK06JVBRw95o1x53ZNn0zzZ0njqID816R3cDJXV1MzR6DdAC+3A9D885vXn2PxAFOACAejxMssncG\ng0HinWDvbENyD+h8p5iXXy4UciiM3e1Eq03RkXMjjN80TubmDOM3jVdMneoU2tWpvBaamYq2lvtf\njQUmFisVcrC+KzLL1R7YhXNVwMaVKwu80u87c6xAyAFmgaV+P2pgPVy1rkDIwXun1Yq4Wb/WLpus\nM68NWf/qHW6jlrmNXNbqKGzVMjiotd9vXZ/fb732kkavfy+mmvtfbQyAckFxGh3VbrFRz+/ITThX\nHhrWDJduanh4/js7MKnZvlszNKzZvluHPuTNunNknbkgVEfOvJtKWabjuTnLVJZv5nVKlWpHo83l\nzTJHe0mlNpe7v8lka6+vWrN8uaWGiykXfKOplKLWi/qv6RpB99rnPY+PrWODmWb246OwZH703pPx\nse2c+tedi5ldEKog37wL83PDxWZet0vEGm0ub6Y52u7ctaQpddPmcve3EddXzbW4WTufTzwOPT2l\n5X5/Z3j1dwqVUtR6Ub/+537LwS0PNWNlV4tE4ORPjBUIOcBRX5PXnbsZvrfLJmZ2oVFUMp/nzLz5\noR3zN79f61CoeeZyr83RyceT2rjN0OoWpY3bDJ183P4C3JianUyebtpcKZOal+b2as3mtdzzwcHS\n47u72386pZMoN53hxX0+Xn+RGZ2BeTO6sjHB58zw9YKY2QXBPZXM5/nhPYu92UMha/lKM03ATu2t\nJQypudckem+U6dm8ULXdQRJXJEqc4iqZmsuZPK+5xl2bTROuvtq+rV6GWa3WbF6LObfVKxYWA/nh\nVoup1txuptNs3jvGlH8G0gFC/94PD/ba1p//GTp5xRuBAOPr1rk7uQNiZheEKqhkPs/fH4lYIUdz\n/f+DB909LMy9JuHbw/hu9RG+PYy5t3absdsIa26IDcWY/vF74LZfwy1zcNuvmf7xe9j89/9RYoKu\nZGouZ/KsJiqcYbg7th6qNZvXsna+2nMI3lKNud1Mp7n+P0eZ6pqxXNrPmGHqulFeWJu21pnnUTyN\nZucVH/RZZvhmIWIuCNgvU8rhxfy3udfk+m9dT+pwCo0mdTjF1d+8GnWrqknYvVzWlnrkLXDvP8Ph\nMOCz/n7rLqa2f7ZkftspBGpOZMuJVzVtbsayvVo6RNWGn/Wy0yXY4xT8J4fbjlNsbIyjviKzz5IM\nc9eNcfLJ5Ttxkd5eEqtXYwQCKKwReWK190lXyiFiLhynVsemhUD+qAssJyXwLnLZ5h2bOTpns+gY\nSB1OEb03WpWgexlhzT/89zB7YmHhXAAySwqKciPuciJbTryqaXMzIsg1o8PQ6lgCi4FqrGrlcIwG\nt2KGQ4cqd+Iivb2Mr1tHZv16xteta6qQg0SAE7I0ennHYiKdNhkbizEzM0Eg0Ed/f5wz7nSYBM7D\nWGYwftN44xtYhGM6VLtjFdx9t/Pysk77HjVjeV8nLiFsJyr5qNh95xhIw41jsGKGUCbAlnP7K4qr\n07w3kwGMv1jXMh8Ht3PmIuYCII46XvHkk5s4cOBOyEvf4PMF+fQvphl6vvx7FYrMzc1Pol1uPXQx\nbr4PXoqXCOHixjRhw4bSCIs9PbBtW6Gg52JEsPlJ+MMDBXbnoM9X0eydmzMvMLUf8dH9+dXcFelt\n2fdOHOCEqhBHHfek0yYjI2F27vQxMhImnTaPlxcLOUAmM030v1T+qfUta81Eqp0puKeHik4/TniV\n1rSVa+kr4Tavu1AfsVipkIMVJjffsS33nUtOplHvPVCibE4pTfOJ9Pay7ZzVhI4FrJ/wZIDQV1or\n5NUgYi4A4qjjlnTa5IknrmdmJgVoZmZSPPHE9cdN68VCnuP0QIZuX7ftPrCWgpXLWNZIfwa7+elt\n26xkIK3MetboYCC1Uk1ed6E+8gcTA6TZzghD7GQ7I7w6lS453k1K03JEens5+PZ16LeuR1+5joP3\ndIaQg4i5kEUcdQoxTZNwOIzP5yMcDmNm1fOppzajdaEjm9ZHeeqxDzNzxNlWvSRgcNd778JYZqBQ\nhE4IETohhEJhLDNs13TPt6X8CNULobcbTXs1wq6VVluLzL0my6/+OOqUcZTKsHzlS5hmdXndhfrI\nDSYGSPNJRjmDGXzAGczwZ2qUtFko6G5Smi5UZM5cOI7MT1ps2rSJO++8k/zfRjAYJJFI8IpXODiy\naQg8p5jptfs9KdasuZve3tpuZjl/hni8sxzOqqGVfhzmXpMNn3mQ2X//QoGnf8+SYxw90uX4vg56\nnHYEuTnzr82OcAalQh0wAqwbnw/K4jalaSchc+ZC1bR6JNYOmKZZIuQA09PTxCrYd/v/WeObKfYK\nV6xcubFmIYfyI9R2NUV7QSutRbGhGLPfv7lkyV45Ic8tZ6yEk8/FYqWcZSkSsaZ7VtgIOcDMRGG5\n25SmCxERc0HIIxaLlQh5jomJCfz+kO0+/2HoHYLV/1cTCBiAIhAwWLPmbs46a2tdbSrnz9BqU3Qj\nacZacycmDk/A4eocRtzkdd+6dRNveMM1vOUtKa68UvPd76YYHY0uWkF34+QYicAJxryJ/EEe5Equ\n5G28jSt9Vx6fAgP74C13r1nD1rPOauJVtYhKwduBbcBzwC8c9ivg88DTwOPA+Xn7rgWeym7XZsuC\nwHeBJ4D/BP7OTRB5LYlWhCaglNJYXmwlm2EYenIyqYeHu/XwMPPb99CTA95kArHLPV8uIUiz838v\nFozbDM2yXzve21ryuieTSR0IFH6/AgF0LIbevdto9CW1JW6/v5PJSb0ruEvHiOkAgYJ7GAwGdXIB\nZ67BZaIVN2J+CXB+GTG/HNiRFfWLgP/Ilp8GjGX/npr9/9SsmL81e0wP8AhwmZvGipgLjcYwDFsh\nV0odf2BMTib17gdCengIvXt7npCXS7lVxORkUu/ebejhYaV377Y6CclkUgeDQdsHlVMmsmozfy0G\nahHaYpKPJ3X3f79O0/1Swb3tWTJb8711+m719qKHh1VtlXYgo6ODenjYr4eH0Q884Ncf//hgiZgr\nm9sxmZzUvf5ex472QsUzMbfqIlxGzL8EXJX3ehQ4E7gK+JLTcXnlW4CPuGmHiLnQaOwEVSmlB+0U\nwUlhKzA5mdS7dgULRve7dgX1qlWhmh5UNTZjQWKXchScBb3cvUs+ntShyP/IjtDndOjM39V1b52s\nPkotnpG5JeQUbA89RImgO33lne/hwu0MNVPMvwNcnPd6CFgLfBL4y7zy/wN8sui9p2RH7P1lzh0F\n9gB7+vr6GnvXBEHbm7q9xBqRU7IpZW/eb/aDqpM7B7kRefHm95ce22yrhvPIXOnJyQ66yXWQG5EX\nbw884Hf1GTjdQxmZa08c4OyCOusy5dablOoCtgOf11o7hubRWie01mu11mtPP/30uhsrtA9Oa7lb\nTSQSYXx8nEwmw/j4OBGPPa5mZuy901assD++r4mRe9o16prbiGu5QC5uypu9EiAejxMscs8PBBR/\n9Vf1rXboLOw/IL9/zpWTo909DAaDxBdrQIw8vBDz/cAr816vAg6UKc+RAJ7SWt/uQRuEDsM0TaLR\nKKlUCq01qVSKaDTaNoLeSAIBe3HeuDHU8gdVuyx1y1+utHSp+4hrTsvD7MqbvRIgEomQSCQwDAOl\nFKGQwdKld/Oxj21d0FkK8z/LuTn7D0gpv6slscX30DAMEomE5x3ujsTN8J3yZvZ3UegA91i2/DTg\n11hOb6dm/z8tu++vgW8APjfnz20yZ75wWIzmshxOc+Y5J7hGmvgroZS9mbqZln4787fd5vOVTgdU\nM2feypUAi8Vxsfg6P/7xQf3QQ6Vm9tHRGrwUFwl46M2+HfgNMIs12r4B2AhszO5XwBeBXwF7gbV5\n770ea8na08CGbNmq7IN7H/Cz7Hajm8aKmC8c2t2RpdHzxnbe7F5Qb7vbYambUxsqbTkxdOvN3kpB\nbYf73EhynVJQGgwNyQJBf/DB3Ny5X4S8Ap6JeTttIuYLh3YemXfqqMmLdtvVoVRty7tqxck64Gar\n9uvTKme/drCA5OOlRchuRQgECwS9TfrsHYGIudDWlFtT3WraadRkJzZOAuRVuwcHS8Wm1s5MLWJZ\n68i8k0Sivb5j3v4WnTrq1gi9ddfZqYiYC21Pq+eHnWiXUZPdKLmnR+vubnuh9ardtQpNsXAPDtZm\nKbC7bp/P2nKm8xNPtG9jKFTdtbaKdrL+eG0lc46iqFp6nZ2KiLnQEDp5DbJb2mXUVM0I1TCcjw+F\nqvvMaukUOJnna72Plb5nyaTVsSmuu7u7c76Tjf4tua3fa/+VciPzhfrMaCQi5oLntNNoopG0y3VW\nM3eslH27u7tLRa/StdTSmamm4+GVhSMU8r7T5UUo2FaSH6/f7VRJPSPzasMSC9UjYi54TruMWJtB\nO1ggqh2Z27W7FsFzWt41MOD8nmo6Hvnnruc+ez0dUm0o2HbDzZI+u8+9FvGdnEzqhx8OlSwxa5cl\nlgsJEXPBc9plLnmxUO2cuR21fGZOnYjc6L/a9zi1tV4LiNedy0qhYNuhg1cON50/p8+9GvG14qsr\n27Csw8OLJ858s3Ar5pLPXHBNubzagj31hKy1y+e9bRvcdZf7HN+1fGZOEdC0do4EF49DUfA6gkHY\nuNG5rfVGm3M6Z60B88qFgm1mmNv8iGnVRIZzE7nO6XN3G8I4nTY5cOBOmI/MXYJTuGKhwbhR/HbZ\nZGTeWtplLrldKfXmbv3cYS2fWbkRXiUnuEY72pU7ZyhkbbWOnMuNzBs1xZRMFk6FnHhi9T4OOSqN\nzL34rTolCZKReeNAzOxCI2iUqbFREdGahb03t1Eg5Lmt2YFx8p2icoJV7rMrt8yt2uVp9TjaVVNX\nMlk6/VCtZ3u5OfNGTDE5eeTX2mkot6Kgnt9q/ufw0EPO5vX8OXPBO0TMO4DBwUHt9/s1oP1+v33O\n7EWA3Rxcpz0U7IWpfULW2glHT4/zA76WwDHVWgHKHV9tXU6OftWuO3fyZvdqZJ4vjE6WgHo6DV53\ntos/h+3bnUfmjzwS6qjfbKcgYt7mDA4O2j7oF5ugT04mHZ1p2sFcNzmZ1I88Mu+1+/DD9g8s+5Gb\n0RYjc61rE7tqhaEWwfMqml05IfSCWjoqdpH73CSQ8brTUI+oF38OAwNJvWNHsOi3qiS+egNxK+bK\nOrYzWLt2rd6zZ0+rm+EJXV1dzNl43Pj9fo4dO9aCFrWGkZEwMzMph72K9eszTW1PPum0yRNPXI/W\nR4v2dLNmzV0FOajDYcspqhATpaJoPe/lFQwGW5KyUSnnfV49Anw++7qUgkyVH6NTXeB8Die8uj7T\ntJzzJiYsR7J43N7xMOcsl+/cpxSceCK89FL15w0Gyzs5ujl/tXXksPscBgZMbrzxU/T2TjDnD3Hu\n6i2LKB9781FK/VhrvbbSceLN3iLshLxc+UKlnOdrV9dpDT9/Oc/hsbFYoZA/OABXboe33c++NwRI\nm+nju+w9qyNs3NiY3Mu1ejw3Ei9XOzi9Ryn7aw2F7I93Kq+FSMTKt10p77adl77W1Qu5m9UKbs9f\na056u89haCjCVRt38baH4T0/+D0PPld9vYL3iJi3CL/fX1X5QiUQcH7SN9JoZJqwfDlcfbXzcqOC\njsaDA/C5T0L6DNA+SC9nNDp6XNDtlpElErB1q7slP9W2vdplUs0QOy+XisXj9qNtre1FacsW6Okp\nLOvpscqbjZslYm6o1Gmo9vy1tMvuM6X7ZRj4FADTs9PEhmroJQieI2LeIqLRaFXlC5X+fucn/dzc\noYacMyeGU1Ol+/JHMAUdjS/fCDNLCo7NTGcYi40df+125OaGdNpkZCTMzp0+RkbCpNPzSl3LyGvL\nFujuLizr7vZW7Jw6NLXch0jEuTNnJ0qRiLUGv3hNfr19p1osIF7EXSg3bVCJ0xwMWk7l5l6T8O1h\nfLf6CN8extw7f5H5nylkYNk4XPEReP3248dMHJZ15e2AiHmL2Lp1K4ODg8dH4n6/n8HBQbZu3dri\nljWX3t4Ifr/98LDcqL0e7MQwn5xY9PfHUSo73Htuhe2xMxMzHrfOEvLR0WjWl0AzM5NidDR6XNBr\nGXlFInDjjZAz/Pj91muvp+697NBYAlKKc+CT2s9tmibLl4dRyodSYZYvN9m0qbZAMU5WBbAsIfkd\njoEB++P8/uZMnZh7TaL3RkkdTqHRpA6niN4bLRH08XEwbuuHP31VgZAD9C2TqFHtgIh5C9m6dSvH\njh1Da82xY8cWnZDnOOusLfh8hbY8ny9YdtReD5XMjTmx6O2NcPbZ2+jqCsEK+4nBQF/A49ZZc/WZ\nTGFvI5OZZmwsVtC+YsqNCE0TvvrV+Shnc3PWa7eC0eg5erv6vY7wVny+pUstUVXK5Oqro0xNWZ0n\nSDE1FeWOO8ya5p4jESvyXbGgB4OWJSS/w/Hgg1Y7ijl2rLY5boBDDgYtu/LYUIzp2cKLdDKdxwfi\nBLsLP5Bgd5D4QGN+p0J1iJgL3lHjE7+3N8Lq1QkCAQNQBAIGq1cnXHnI1hIutZzoFYtFb2+Eiy8+\nyJp/eCu+YOHPxRf00R/vr3i+anFyCsyV1yJy9ThFNTqUqVP94J3Zvvh8H/4wvPxyriQGFJtqprPl\npbiZe966Fe6+213b59tRSCpV2z2uprPnZCK3K4+cGyFxRQJjmYFCYSwzSFyRIHKueLK3BW7Wr7XL\ntpDWmS84Ghzr1ctUi07rfUOh8s2dTE7q3cZuPayG9W5jt55MTnpybcU4hczMX3ffzNCpjc6W57b+\nctdcX+Q5++A+VnnjrrvS9UP5wD5OVPNT/NBdIb19B3roIfT2HeiBL6K5BW3cZnhxaYIHIEFjhKbS\ngCd+cnJSG7t3azX8kO4dvkfHhgcKIsStWhWyfQi7CcrSzhmwJieTeteuoKcR8er5eBqdLc9N/dVG\ni+vudo7TXno+w0HMjaqj4LmlOKZ8ubCu1UaxK67f6fs9OZnUQ8M9Bd+zHQ+iL/tSt04+3kY/iEWO\nWzGXoDGCN9QRMcQuGAdvTxMdHWU6770BjvBJPsfbGQLgbW9zOqUiU22UkjYjnTYZG4sxMzNBINBH\nf3+8rsAcToFMNm60TMLlsA+IY5mOx8drblJV9Zc7Buz35ZMfNKW0LhOIUmhqD9LdneDGGyPcd1/l\nQDFuyH3PUynr3ud/d7u7YXbW+b2NeEw7BWw65gvx9ksOen9CoSYkaEwn0Y4RQKqlxogh8/OlJlqH\nSaV8XHNNmD+5/Y4CIQeYYQlf5sbjr1fYO5jTtwBysvb2Rli3bpz16zOsWzded4StSASuvbbQKUtr\nd05wjXREc1t/OQ9+N3PY+f4B8bj1U5snAiQAy2cDDEKhBHfdFWHrVm+88zdtgmuume9EFItzOSFv\nFE6+GV2ZxiwJFRqLiHmraWai5Ebi4omcNtOMhEfY6dvJSHiEtJnOOmblRkaWN7HWKV6+7e8sV98i\nnmNewTduDBEsOmcwGCTulcosMO67r1RE3Hpne+mIlgvYo5SJUmGuvtqHzxcmFDId6y/XV3Tbd8uJ\nfiQCX/uaFV41h88XYXBwHK0zaD3OwYORkutz2+cuPm7TJrjzztpH1/UG9nGKWeC09LNRS0KFBuPG\nFt8u24KcM69xMjOZTGrDMLRSShuG0dQc2WUa5ThRN5mc1LuCu/Qww8e3XcFdeoBJ5znL3l7N8HDB\n1ju8vWAOuS3vQxOpZu6/0XPfbpjP3pbUUOi82NPj7LxY7Zy5145rbp3KyqUhrbSFQvWncS2mnP9F\nI3wzBO9BHOA6hBqesLV6cbeS3cbuAiHPbV/379aO3sRKFQh5YPh+HRse6Mh8542g2rSmjfZKL8Yu\nneh8G+w7cOWcF916s9s5lFWb4awYt/eunGd6uS2/Y+KlY+YDu0JlV0bYrRIR2gsR806hhiesYVT/\nIPSKWkfCw6pUyIcZ1g8xrJWyv57QqlVZb/ZhbezerZOTtS8Fa2fv9VqpNq1pg1cPFjA4WEnAGpvr\n3e3n7faelLuWfNyOwtYD4p0AACAASURBVIt/6o34DJKPJ/XQQ/a5x4eHm2iOEepCxLzNcOwB1/CE\nVaqxD0In6rEIOI3Mdxu79eBgUivVOEtDM0WsmbgVmHya1anJjcidt+Z1SMtds9u+tNP1+P2FxznV\nVyzySlkdnkZi3Gbo7TvsxTw/ZoHQ3oiYtxEV56aqfMKWG5k3cg65HouA05x5LvBKY9vt7oHdCuoR\n11rEvFmUa1stc+a1Uqkj53aWy+29djpfboqhmZYhdYvSA1+01o4XryUXc3rnIGLeRriJ6FUNTiPk\nwcHBhs6l12sRaFYEtWLawfHLjnotBtWa2b1qsxtRKjeSTSZzbU9mR+hKh0KNcV6s1JHL3z/ApN7O\nbj2U9eXI/35W0yH0yvpRbz3GbYbmFiuqW36Utw/d1cAviOA5IuZtxPCw8nzeanBwUPv9fg1ov9+v\nBwcHGz6X3sq5+npo1cjczgHMTbtCIffzvV57P5cjmdS6q6vwfF1d9udzmjNvtGm5mHJz2Pme8ANM\n6h2Usxw1dqqmWLgHB+s/X/LxpA7Gg5pbOL4F40GJ7tZhiJi3Ec0amdt6hFcxcq71vO3sRa91a+bM\n3YiZW2epej2xvWLpUvv2LV1qf3ylzowTXl5TOe/yfA/yr/udfToa0a58qlnOVm0HNPl4Uhu3GVrd\norRxmyFC3oGImLcRlebMq10e4jRCzo3Ui7dVq0JV1V+uPZ26rrvZ3uxuHKaqWcbktRWhlvtRrn1e\ntsvLjlelNei5++q02mJYDXt1aQVtyr/3TtMldlurp4aE5uNWzCUCXBMol+IznTYZHY1mYyRrZmZS\njI5Gj0dpsmPCIX7l3NxcSUS0E07o4brrXnRdfzptsm/fhoLj9+3bcPz4SCTC+Pg4mUyG8fFxIvXm\no2wSkUj9YTmrSbeayxtertwuaJ4TTiFLa4kE3M5BB+tJ1WpHLnqdE7n76pSX3ut89Xb3fmrK/fsX\nQKRioUGImDcJp1jbY2MxMpnCp1cmM83YmPPTyyn2uGEYJBIJDMNAKYVhGPz5n5/EwEBh4Ody9T/5\n5GagOFD0bLZ88WKaJtFolFQqhdaaVCpFNBp1FHS/376e/HK7MKlOoTvtPnInUX7726Gry6qzq8sK\nJ5pPrYKZH9fdTXktOCVMqZRIpRyRyHxClmJy97U/3u9JvvpKnSu7e+9E8X31Mh6+sPAQMW8xTskO\nnMoB4vG4Y0zy4pHz+vX2SROc6p+bsx8mOJUvFmKxGNNFT+Hp6WliDgoYjdrXU1xebDHYssV9UhMn\nUR4amrcAzM3BHXcUCnq5pCXl2LixuvJacNMJqoVKqQN6I72sTqwmYARAQcAIsDqxmt5Ir2198/Hl\nrW35cuseV7J4uEkKk2vbxo3excMXFgFubPHtsnXqnHk5nJzjHnmk/Dz38blr0Ibfr5O5CcCiycVq\nne/sve6tbTFTy7K8RjuAVRNtzM1cvZt5+VqvyS2NnJf3cslY8SoCt/4O9a5gaNQ1Ce0Lks+8M8jN\nmeeb2pXqwfpc5s3dPl/w+Dz7ceySVOcnbnao37auLD/4wXKOHSsdhXd1hbj44sWb4zgcDpOysfUa\nhsG4F0m9a8Apx7cTuZ+6i69Ny2h07nQvqPa+K2VZXsDdvTf3msSGYkwcnqBvWR/xgTiRc0s/mHb+\nHAXvkHzmHYKdc5zPdxLF89a289xOdtarrz4+YVfO+c6O17xmC0r1FJQp1cNrXrOlruusmxbnfC83\ntdEI3FxuNQ50lebqqxWARn0cjc6dXg5zr0n49jC+W32Ebw9j7rW/KLem8hz5/g6V7r251yR6b5TU\n4RQaTepwiui9Udu2eO0sKHQ4bobv7bItRDO7Ha6DzFSys5ZZ0zM5mdQPPxwqMOvXulSu4VRar9Qk\nW2OzluVVszyr+NIHBuy/Cl6aw5sdQKUhSUiKg7T87SOuA6yUW1JY/JOs9r7korYVb8ZtRsmx7RrZ\nUPAWxMzePqTTJmNjMWZmJggE+ujvjzuOjAFGRsLZpWGFBAIG69aNzxe4sffZ2Cdzy8+KR/9K9XD2\n2dvKtq0llLO9xuMLztZYr6l50ybr8ufmrBF5NApbt7ZP+5qBaVoj1IkJOO00q+zQIWuUfPnl8NWv\nFn5lVPc0+oob4fXbC+oxlhmM3zReUveGDTBbtOijpwduuAHuu886b1+f9fWs5mvou9WHpvSZrFBk\nbs4UlHXC5yDUj5jZ24Ry68jTaZORkTA7d/oYGQkfX8vd3x/H5yu0Nfp8Qfr7i2yNbuysNjZBy1xf\nvPwMtD5adklcyyjnfl3J1thi83yOappRq7d5jq1b4dgxa5x27Ji3Ql6uHdWanxtF8ZK9qSlr01kP\n8zvvLP3K6NkgDP1NSV0Th0svKhKBu+4qXEYYCsG2bda9rieeQd8y+2WnduWtnJIQ2g8R8wbjtI78\nySc3O4q863nu/Ak4J2wWKJdb9lZuX8twipTR11deWdokOkq1zSh3uY1oW7V9nWa2rxYqreV2NEYe\nLr0AJ3GNRODgwXnj9sGDtQYiKrz/lx9JEuwu8s3oDhIfKFVoL3wfhIWDiHkDSKdNfvCD5ezcqWzN\n5WCt2y4XLMYpyEwJuYXKyaTrbnog4PzULbevZZQbgpRTlhZ4CNlFiau2Gc0acdXa12n1iLBSB8S9\nhcAEwliPwTCc8IWCvU4iWgt2bba7/1/9zMVc6/sexjIDhcJYZpC4ImHrzQ7eRDYUFghuJtaBbcBz\nwC8c9ivg88DTwOPA+Xn7rgWeym7X5pVfAOzNvufzYM3fl9s8d4BrgKfN5GRS79zZU3a9dvmtDu8V\nl9czOZnUw8PdJefeubOn9c5uOezSSNldWzlvrCo8hLxwbnNOgJOs2lGpYU5geRUb/mfKrosu14ZW\nOKnlZzkr52TmLuZ9aT71ru6ADkVCniclGfzGpFb37NYMDWu279YMTGrQ2uezb1ubJyEUmgxeJloB\nLgHOLyPmlwM7sqJ+EfAf2fLTgLHs31Oz/5+a3fcYsC77nh3AZZXa4amYJ5Na9/To5Llo4ya0uhlt\n/Ck6ubU+t1+nIC3FSVYeeSRku6/WTGrVUs6bveVU6y7tpCwuo6N4lQ3OOQGOUfNDu9pORnLwEW34\nn9GKOW34n9HJwUfyKyu4r4o5R7FrRaa5gutwOL9TUpL8e1kpuUowqPXSpfafVa3pfJ2ykyUnJ7W6\nf5dmeHh+27HruKCLN7pQCU/F3KqPcBkx/xJwVd7rUeBM4CrgS8XHZfc9kVdecJzT5qmYh0I6eS46\n+KnCJSDBGHX1yJ2Xlc2L9eRk0sqkNlQ4gt811EYj41biVQJyl4rkVZ52pyhxoGoSxmo7GcnBR3QP\nRwrO08OReUEvuq8Gv3YUEyfRDIWquiU1U01GOTsBzO/fhULWlt/XcxvRz01n6huPDep7dig99BB6\n+w70wBfnl7UZu3cXCnlu277b8VpyX7uBe0c1Dw5rHhrWPDisB+4dbcStFtoct2Lu1Zz5K4Bn8l7v\nz5aVK99vU948pqaIDcB0YXwUprshNlT7nGr5+Wjj+Px374Ow+nOawCSQgcCk9br3QdrGA7tleOUu\n7dJDyCkLnVO5E84JcPpqclRyEw8+/6vy4TvWcZTCLF9HCbD5S6/NXVDBvjifQlG43AksSXHK5DU1\n1ZyvYz2BWaBwLvngQWvLn1d2+qzyy90k10mnTU783Z30LtH4FJyxBD65GtadMm1FcZuZsW/wCvvy\nnO/B27/zJEMnHgA/lu3SD0MnHuDt33nS/U0RFhVeibld3iRdQ3lpxUpFlVJ7lFJ7nn/++TqaWMrE\nModym+UodtgtLevvj5dEULPoLlxaFovRu2OWdVfB+gFYdxX07piFzZu988Du1E6Bl+7SLjyE3DzY\n3VA+AU71jkqVOhnFDlQZ7LORTGVOtf4pup4I2+1/dBVoRoQxp1sfCnnjfOcmop+bztTYWIyAr/Au\nLvHDjf3Wc6Qv4JBC9bn5cr+/tJM3dMKB0qekypYLgg1eifl+4JV5r1cBByqUr7IpL0FrndBar9Va\nrz399NM9ai4QCtF32H6XT/kqhnR0Wj8OcPbZ2+jqml+E6veHWLPmrkKPdKehx9SUNx7YbbIsqyaa\n7C7tVajWSCRSkoI2kUjUnPPdqTNx4olWeTXpNAHb+2qo/baHOqViheasJ3f6CmzZUvtyrPy+bSwW\n4dprSz8riBw/JpWqbLFxWsq5ImAta4v39xP0FT1mj/jgy/3Hr+mrX7Xp5Dk9mWX9keCEG1u8rjxn\n/i4KHeAey5afBvway/nt1Oz/p2X3/Sh7bM4B7vJKbfDaAS75Bl/JnHnx5hTSsdpsZCVUOymYm0xz\n64Hk1bxzq2hyOqhmhWqttk1dXYVz5pYHdlIPDrrPmhY68ff5lRbc1+TgI47z+W6czRp7/d59Bdy4\nT5QeU9mXwuk5cM8OVeAEZ+zerdXwsA49sFuHPjRZ8Zp40GaefdiaOxcWF3jszb4d+A1W2LD9wA3A\nRmBjdr8Cvgj8Cmu52dq8916PtfzsaWBDXvla4BfZ93yBFi1NS64PHfdm99/icx0X2XX89DLnrspd\n16XnVEFqVLBSo5bzFKoFybtYE7XcNqWSWVFR2b/WMje/321/MFMxLrtTu1rt0V6Jau6nm75t6TGl\ny9eKHRAnJ5N6165gwTPg/iGlv/FYfatiBu4dtRzf8oX8IXGCW4x4KubtsjU60Yq6RdmKubqlVADr\nHplr7X4hrcuhka33c7Gg1zusqvUJ36AOQDuOqO3IroTMEwlLoEOh8m0u9zWw+yi6uqru/1Vsdzv2\n26r9GroJOWB/jPVZlft+NSoxkXizC1qLmNdENRmL7Hrku3YF9eRksmTN6eB3Bm3XoDqSe4I6PcUd\nRteOS6y8HFbVYr73eoiXvT9J0MGiJUa1rA+v9fzVKNy8waXyaC8fv9/+dvv99k2p1TRebFYPhdpH\nuO2o9mvodPyHPjQvxF//uqEHBkoD/HTKzJSwMBExr4Hk40nXaRC1tu+R29Xhdh6+hCqfWI5rZ3Pv\ncfl0LjvarSXvopfz93kdA8PmWovnND2nxo7J/OEOHS6HNg8O2t86J9N5LR9PodVgfuvubl9Br/Y6\n7T62yy5L6qGhwg75jh3BAkFXyrrX1fbfnILICEK1iJjXiOOP0OWv2Wl072a0X9qY6oTDi+AnFQOV\n1CLMlZ68g4PzQ1C/v3zy7bzzKwcxLw78URWVPucaOybzh7oLVpJPjbfHdfPKGYHadVRaq4Eo/6N9\n4AH7qbLt242COnt6rI6N2/5btYMCQSiHiLmXVCGqTvPububhHc/tckjgRVjSih2CWkam5Z68dQw9\nDQcxz7W1nIXBdp+ba6tl6KvzTdj1d7jKYXcJSmn92tc6dwjKecW3a2hRL2ZunJxYh4aU4/2w6ziU\nTHVE/kftHXhBKELE3EuqGAZ4OjKvgXodwpxDkjJfl0MHo6xVw+nJW2lSuJi8zyKJ5eBn13lJJpP/\nf3vvHubIWd/5fn7SdGuQLwPWeNrh0iU6hFmctQ8bJuTxc06Obdoh2ByH6xoTGYyNo3XPyXpM1jlL\nopNgONHB7PoJ9nIYE2WwcVAxDixnwT6PTbDbMxN2p1kwC/YAToPTaTUOGdluhwl2Z3pmWu/5o0o9\nulSVSlLp1v37PI8etd4qVb31qrq+7+V3MWNjY3XbxsbG1reNj483Xd+M11xz4+/c4ci8WKwm1mhe\nMx8fj3adP6zbWlXQ+zkyj9rdrJtjPXzIOz9C48g8qKPjaa869oLhne/tvAOvKDWomEdJO9m3olwz\nHwB+I/NWo/yWU4t+T96gp6XnieqfnkUwlojjilfTeUmlUp7XkEqlfLeBhytf4+/cxZCw3gjOMlV3\ns1Qq2nshbAiDWiO6fqyZe51nfHww6/LFJ4rm8j8bMw89Ui/kD33tJZ5GcH4dHd+23vZ3OjJXIkHF\nPEr8/mNTqVAj1Lat2QeI11R9mOngdjwB6mh3ZO5UsuWQzK/+rV4pv6d24/kbzcZDmH93OEPfNmGD\nydT2l/phzd7L5C3tjtKr9+v0p53kKNUkKW/7ozc39dWC1sz923rNs2M7Kq6UyvCgYh4lXqOxdq1i\nRohisegrdn6GWmF89D39cdtdMw9Jp2JOYz28ftMOh7JRGvUH0e7IvF+E6VR0QieTJU336zvf64ym\nWfPMsubXWfDt5//Cz5s68FGl2lU2FyrmURPWoTceH74IGx3QrmV8q5F5kF9+W+baIel0mp2qugb9\nhh0uMkftbh90Hr/l/wj7S23TKzHvyIK/9n5953udde4Ofpewv2mxWDTxeLyt/ylFMUbFvPeEmcvs\nw0g9smm7hs5KcWamvVzaLdbMI4mY19blNBu5jY+PrxvA+U6zn3lm64N3Yf4dpQFYEEFRgSPqL0VW\npzDT7EHt1pFvfe39uu3v2u4MhK2bsz146aorV0plw6Ni3mvCzmWGCb3V6unuZz0e1bSdz/CiODPT\nVkchKFBG17HsOyCoozMzM9P0UB3fsiVc242AY3a/1ufboVhsXpkKY2TXavTb6fJF9X6FtUjayn8q\n3vIVch2ZK61QMe81nj4pbT4Rwpj3BjzJ2p0K96UPi7n9Hpl74jH70NGsRr/Mv7sg6Cft1+yA13k6\nOXer27Pb5Ysobv+gOgS5e+qaudIKFfN+UH0yBYl50Bp6mHnHgCeNb/jWhg5Ey6n4PgzjAtfM+0HU\nC9Zdmn/3KjlHbfW8LndmxiOozNiLhnf+dqSeFlE2d5jbs5sOShR1DeoQ+HW64/G4CrnSEhXzsEQx\nTAk75d74hAjat0rAkyzMyDzUVHyfzKx7LWCB9MuUPASeHZuviTn65WgXsr1u7VZ+0WFjIMzMzKwb\ndMXjcTPTsAgf5M3ZLu38dJ3aUnb7GAjqcKgVu9INKuZh6KRL7jUi8xru+L1qn0BhxDzgSRbmIRFq\nKr5fZtZebdmP+V5jupt9iLievksO90nP2zyMX3Sr2ABe9gZAnaAH2QhGZNLRdJweeTmGovVSgPqX\nK52hYh6GdkdrXhY8tUOA2ge+35OsVjzCTLN7PcmqT8oQ676BmdRq695PYfW7rl52IDq2koq+nr7G\ngLMh6tMlrUbmYcKO+rlYxWsc16O2EQxze3YSf6gbikePGuvwYSMHDpjUw4fN2OVH+94fVjY+KuZh\naHe0FvSEapw/DCMeYc17a+dHG+vc4onhOzIfxPClvmLRPenD0Kko96CeviPz/f73XlQju2LRmPGt\nJ+svpyGWeKuReaCPfs15wvRno8TvfDXViozi0aMmeeiQ4cCB9df47CGTes/RjRBmQhkiwop5jM3M\n5GR75UtL/sdaXq7/nM9DMllflkw65VUyGbjnHrAsEHHe77nHKa8lk4HFRWe7MfXbVlYgl/OtVj6f\nJ9lQjySQb9zxz/7M9xg9wa8tg9q4GzIZKBTq27pQaG7rsPXpop5TU3liq1JXFjsOU/vwvPds2yab\nzVIqlTDGUCqVyGaz2Lbd/skvtDFX/g5sWwQqzvuVvwMX7gcgOZYkP910d9QRj8dblmcykEp5f9/v\n36tbfKrlWx4W24btV5eR++aQAwfZ/sgce378Y1Yqlbr9TsQqnLlngUrF+XdtdWspSpRsbjEPI7i1\ntPMUCiserlDbn6+QZpHY+zKk084DpIkOhCWTyVAoFLAsCxHBAgpA03OmUvE5aY9otyPVLbYNe/ZA\nqeR0iF54Idz3elDPiYkMO1+8kURZoAKJo7DzdpiY8773crkcKysrdWUrKyvkAjpxfuRmc5z85c/B\nh14Nt8add1fIrW0WhSsLZC4IVqFsNhuq/M472/v3CqJctpmbS3PwYIy5uTTlcvO96lMt3/Iw2DZc\nZ5dZ/sA8nLcKAstbVlk+ecpz/6XV1c5PpijdEGb4PiyvgVuzB80ddpEtIvQMcBRTvn4Li72c4vai\nn2vm3fiF97KeIe+9sC6IYWgVQ7828M977km5aUKbvQ9aWbO3eYmBeFn/f21WzJe/1XzOqCMDW5Yx\n7D9cN50e9LIOH+7uhIrSACGn2cU0TtsOMbt27TKPPfZYf09q28409tKSMxp7zWtgdrZ+n/FxuPvu\njufV0mlnwNiIZTnTdXV1yWadqfUqyWS46eIqu3fDXXd5bxNxRuj9orFt8/ng69i927nWtTVn7jSb\nhb17W5/Hr4HBo5EjqGeE2Edsrr34Wtb+ca1pm2VZLLaqewPpO9KUjjW3hbXNIj+dJ/tAlpWTK0yf\nC7fshK01U9SxWJKdOwtMTPTp2t1mv+22NOed11zn8nHh1Cs+z2U7YGEhx+rqEonEJFNT+ZZ1tMtl\ncgsLLK2uMplIkJ+aIjMx0bRfLAZm9iBI8zEaScZiFHbu9DyOonSKiHzHGLOr5X4q5gHYNlx/PZw4\ncbpsfBw++EF48MHIHu6xWPNSOPhoaxTCcuaZ8OKLzeVhhG1Q+HVCZmZaC7pfA0P/OzBtYB+xHXH9\nzgo8AJw8vS2ZTFIoFMi0+duvH/Pk6Q5hcixJ4coCudncutDv/zU4b2vz9xMJi4suWuzgatqjtt86\nOxsjFmv+/SoGCqUUu6f+mUrl9PW06nTY5TLZ+fm6NW8/Id5+dZnlf/Okp5gLMJlItOwQKEo3qJh3\nSblss3Do/axur5B4xjFOmqgOyFMpeO65yM4VemQeFVGM8PvNli3OiLyReBxOea9frtPtyHxA1I2i\nnwBmgWMQf1mcez91b9tCXsU+YpObzbF0bInJbZPkp/NkLsgQ+2gMg/M8mP1fIeY5GhUuuaT3nZ/a\nn2z/fmdk/sgjsG8fPPMM7NgB734//C9var/TkZ6bo+Sxtm0lEixedFFd2fZH5lje4r8Obi65JOQV\nKUpnhBXzzW0A50O5bDM/n2V1RwVisHoezN8C5Wl3h+XlSI3F2rXD65pOLbsHiZeQB5XXks87MyqN\njI31sJG7Z+lYjWHjhcCHgFuhctMamWuuARHsM88kvX07sViMdDodysI9c0GGxZsXqXykwuLNi+sG\nb5PbThv1PeOjX4lEjwwUG6i16dy3L8+DD45x++1QLjuTLOUy3HUHHPmG9/ePH1/y/Rf1M1LzKn8+\nQMjjOKN8RRkGVMw9WFjI1U3bAVS2wsINNQXZbGSC3rG22rYzhInF8DeBDzjp4iJt+9F0c85u6Mbv\nKJNxbBpqfaVSKW83wCGiVlzryo857zaQffFFSsvL3busAfnpPMkxp1e5bwGON/STYrEkU1PhOj/d\n3ia1zgKzsxk+9amzadTayknY91nvR1i5POn7LzqZSHif06Pcb1+ANSA7P6+CrgwFKuYerK56u3qt\n7qj50MK/u10atRVs0um0/4irOlVedbUqlSLtYHgyiHNW6dbvKJNxlkaqNunPPTfUQg714loleQLy\n7nJPDlhp+E6nLmvgjNgLVxawtlk8+qzwuZ+kOBVLAUIiYYU2foviNmmcrTp+/HnP/Z59pkIsVt9G\nx48n2bcv7/svmp+aIhmrf/QlYzHyU1Oh9q1lpVIht7AQcCWK0ifCmLwPy6tfWdMCI3R1EcoqbBSv\nYUqOUsegk5VE7Xc0AtTliL8ZU7zgdLuLTxS2TlzWOqqbb/7uaG6T2uPH45bntVqW4zK3f79lZmfF\n7N9vmenpYst/0dpQrNbhw6Z49Kh/Pdx9/dzR5MCB9i5MUdoADefaOZ5ZrR7CHJ3u/OnUTuakUMlR\ngjJZ9CKW5CDic/abfsenb5cGlbR8xNzqQ+cqOH939LdJq/+ffvQz/QRdfcuVXqJi3iV16TofTpmj\nlzfEUG8zcEgogXYJFSTEfXoV3Ye6uO/FDusXiNeTexAj814yqMxx7dAQ/KYIJtlwj/QqtWZjP8cv\nR5CTv7v6uWjAMiAGLJNKdVevoJmtfvx8XvHYk4cOBY7qFaVbVMyjpstRWztRvMKmLS2OjTU/zGsF\nPSqRDUowM2yC1ymDXkIIS0MK3uIZZxgrlYoktaafWLbqyzWOvp38QUUD9SPp8fHe5vDux8RKO9Pz\nihIFKuZDRjsj87BT8lYq5X3MKOY1a4kyOfWw0ou54X4QkYIF3XNBfTm/vk8qFf5+VxTFHxXzIaOd\nNfPq/nWjpMZ86cVi61zlvR6Zb6QH8xBdY+jRX5tzy3VLRw2x1oM6m0F9Ob9TRxlPvp3rUJSNhop5\nEAMydOo4J7XPQztwZN7rNfONMr1eZUiu0WtdVma/Zma+9eXmnf06IB5JfzyNOg8l14UwSHyDTuP3\nb+TXOUilrLrvePRRA2l1HYqy0VAx92NIHtpt4fM0LYJJNjyEk2CKqVRvrNmDnrrDbgkehiG4Bj+L\naXnoPlN8oqE+bSx/+LpbHrac8waMzDv5l/GaiRofT7pr6f7VbnXcVtehKBuNsGK++YLG5HL1Mckh\n8gAwkeORr9zGDRpiDNUYaJZlUSgWyfQiIEpQxDjbxrZt0rfdRuyRR0jfdpsT5Kaf+dGjoNOoeBHi\nF2rUJM4ld/+e+sKgnOoN97NvICS3PJ/Pk2yIKZxMJsnn8x1FKMxkMhQKBSzLQkSwLIuzzipw8mRw\nm7b6V2x1HYqyWdl8Yu4hjIHlw0DDQ9sGskA1dcgakATyV1zRcfKNbrAfeIDsTTdROu88TCxG6bzz\nyN50E/YDD/S9LqOOb/jQ1WdYOrlcXxYUV77hfvaLqV4t9xLf2qxsnfRzMpkMi4uLVCoVFhcXef75\ncPdm0L9iq+tQlM3K5hNzv9FM0Chn0DTEtvQM4wnkPvOZgYyGc29/Oytb61NXrWzdSu7tb+97XYaS\nNgKV56emkFPH6wvXjsPCvvWY7OtkMvXx5mtpuJ+npvJNYU8bY603im/UHcOw/2JB+1Wv45FH4Oqr\n4U1vgquvFr773SsGlzdAUYaAzSfmfU9RFgG185yA38BlyZiBLBcs+eRw9ivfVLQZqDwzMcGNX/0M\n8s9HwVTg+FGYv53k38+S/56HcN95Z6j7eWIiw86dBRIJi3ZjrUeF179eI63+FScmMnz/+9dy++1S\nk0HN8Psf+iz2IWtM6QAAIABJREFUddcNJm+AogwDYRbWh+U16tbskVAsGsvH8tiC9vyiI2oH6+GH\nPY22Jh79S3Uf6sTlrVg0xTeMGetmjHwEJyb7G8b8f58Rup8bw+tPT7dfdV9jvSFxLVSUKEGt2Tcu\nxZkZ/8hvYR9eDaFBDTifOxCC4tGjJjk7WyfkiUcfMrkD0+vWxl+bFfPlb238xChNdBqMplOBHmJh\nb2UVH7bqfm50IKbIe9trZ0UZclTMNzjFmRljidTHZG/Hxc4vuLaHj3Ko+tQEOpk48KU6Ic8dmDYT\nB/YbHp01sa/9pbfPdNsnHF7RqqOfwWiG3O0yKJJcKtXct/Srut/IHCyT5IV6QdeRuTLiqJhvBroR\ntCBn3y45cEDqhDxx4KH6KfhHHupO0IdVtLx+j37WNaDjMAxR08JGkmvV5/HyYXfiwDs+7BZ/Nzz3\nhKJ0SVgx33wGcBuJMP5CA7DwrXUT2scNrFJv6U58K4V/rHR+gmGMFeBn6AbtO2l3io9PV/k1Jebn\ns6yulgDD6mqJJ598HwcPCnNzacrl/hiJdeIw4nVJVTc6cIz5nPcC4LTpEpO9bWdFGUJUzDcyQZbU\nfi5NfuUhKZdtTp16AWOcz8+ww3O/tbFzOj/JMMYKCOpgdBqMpt2OmI9aLtwYp1JpdGZ0fqDVVUfo\n+yHoYazZG/HrAGQyGSxrEagAi1SFHGDSig0s6I+iDIpQYi4ibxGReRF5SkQ+7LHdEpFZEXlCRA6K\nyCtrtn1CRL7vvt5TUz4tIv9DRL4nIv9VRF4TzSWNKL0YQfsJzLXXwlVXwdhY/baxMbjzTuxymfTc\nHLGDB0nPzWGXy6FOVy7bzM9nWVtbRsQpO9c847lv/OTz7V7NaYYxVkDUHYw2XdoAX7fL1dRa4Kkq\nlRUWFno/q9HgYdmSVm5qo+hlqii9oqWYi0gc+DRwOXA+8F4ROb9ht9uBvzDGXAh8DPi4+923Ar8C\nvB74NeD3ReRs9zt3ARljzOuBLwD/Z/eXM6J08uAOg5+QrK3Bvn1www3107/33IN92WVk5+cpra5i\ngNLqKtn5+VCCvrCQaxoB/o7sY9w0B0HJvqyLSaFhfIpH3cHoZCnBJ+5qYmtr9QwbDtW2bdLpNLFY\njHQ67YTtrd8hsFNanaQoFluP0luHjO3fCoaiDDthnqhvBJ4yxiwYY04A9wFva9jnfGDW/ftAzfbz\ngUPGmFPGmBeBx4G3uNsMUBX2bcBPO7uEDUCv1oCDhOTkSfjiF5umf3MLC6xU6tezVyoVcgsLLU/n\nJQiXMcstcjvxE8+BqRA/8RwzZ6+w91ff2ebF1NDJU7zXtgNRdzA6Hel7TOl7RX9rJEw4VNu2yWaz\nlEoljDGUSiWy2expQbdtuP76+k7p9dd7tnX1J4zHmzYBzk8aLmTswMPpK8pQEEbMXwH8pObz025Z\nLY8D73L/fgdwloik3PLLRSQpItuBS4FXufvdADwoIk8D7wNu6+wSNgC9WgNutUi5vNxU5Jfow6+8\nFj9B+N8ST3Hqze/GXPomTr353d0JeZV2nuK9mvlorE+Uw8QIR/r10d/AMRo7TWNYVz9ye/aw0tDp\nXFlZIVftdO7ZAydO1H/pxAmn3INMBu69d/gmWRRlFAkj5uJRZho+3wJcLCLfBS4G/h44ZYz5OvAg\ncBjYD8wBp9zvfAi4whjzSuAe4E89Ty6SFZHHROSxZ599NkR1R5BerQFXBaadqvgk+vBNAFJDmPjf\nfnS6Th+Kflm/RzlMjHikPzGR4aKLFrnkEsPrXvf59sO62jZLHp0/gKVqp9Nnu285OlWuKFEhxjTq\ncsMOIhcBtxpjftP9/AcAxpiP++x/JvA3rkg3bvsCUAS+DXzTGPOLbvkk8DVjTONafB27du0yjz32\nWMuLGjmqI8dawUkmo3uqbd/u/UBNpeC55+qrUi6TnZ+vm2pPxmIUdu4kEyLWerlss7CQY3V1iURi\nkqmpfEuh6Pac3ge1HbFeWgK/e1zEEd5hpfYaJicdIR+UyqXTpEul9Ux9tViWxeLiIutWj160eM4o\niuKNiHzHGLOr1X5hRubfBn5JRF4tIuPA1cD9DSfbLiLVY/0BcLdbHnen2xGRC4ELga8D/whsE5HX\nut/5DeDJEHXZmEQ4PLGP2KTz24ndKqQ/JNiXbncs18fH63ccH3eSdDRWZWKCws6dWImE48GbSLQl\nqqdHgBUuumgxVCIPv3X6a598srMReuO0uh+dzHz002/fa6RfPb8IbNnivPcjfsDSEnmcVLu1VFPv\nkk77f7dLd0dFUUIQJrIMcAXwI+BvgZxb9jHgt9y/3w382N1nH5Bwy7cCP3Rf3wReX3PMdwBHcNbV\nDwJTreqhEeCCKT5RNMmPjhtuZf2V/EM3ScfMzNCGPxWPJC3VV3J21hSPHm3vgEFxQ7uJwuYRza04\nNmasVMqIiLEsyxR72a5e0eT6FQHPbdOiGz54PYzwGWf41wmMGQtIEKMoSksIGQGu5TT7MLFhp9kj\nIn1HmtKx5olQ62ew+F8sZ3Q3hKTn5igFGNhZp06xeNll4Q8YiwVPrXc6ZZ1OO6N9FxvIUp9bPplM\nUigUIs8F7nX+Jqwe/sZ+S0EveYn/mrhlDXZpQFE2AFFOsytR0eMp2qVj3tbvS9sYbHS0FuSnpkjG\n/G/FJT//JT/8ps8tqzPjtOrv1iCkOeqFHBqsu6Om1W/Yy9/YbynoeZ/gPyLqK6YofUTFvF/0wT1q\ncpu3iE0eY7DR0VpQXaePnzrluX2y3XXzqCzBbdsxHrzmGs8RsZ90LkUhql4dv1a/Ya9/Y681/GGM\nxqcomxAV815TfShfc03P3aPy03mSUm/oljwB+W+MDb3jbmZignv37SN5vD5aXPL4cfJf+UqbB4vA\noLDa+Qpwq/KTq3POceLO20ds0nekiX00RvqONPaRkB03v47fFVf4xw0YlHP2MEbjU5TNSJiF9WF5\njZwBXJDBUvUlEu0pnyga609SRj6CsW7GFC9JdWyAVCwWjWVZ4Y272k3J2rj/zIwpXn65sfbvNzI7\na6z9+03x8ssHY0AVwoiumEqZ8S1bGlJxYuJxzJ6PTZtkPllvjJhPmuITDdfi1WZBOdBrt8fj9eWD\nYlRyyyvKCILmMx8CwlhVeyVsjoouHrLFYtEkx8frRCo5Pu4v6O3m7fbbv49W98Unisb6pGXkVjHW\nJ616oW2VfNv93VKpVJOYA2bHDsz0p6kTc27FWJ+0WrdBnzp+iqIMPyrmw0ArQeilO1G74tqA5SNS\nVirl8wWrvc5Ku/s3XpuP4B89WjSHD1vmwAExhw9b5uhR7+stPlEMHjkHdcRq2lFEPNtJBLP/IUfQ\n9z+EmX3Ueb/p2mnzpZd9yTwqj5rD8S+Zo0w3H7864g4anSuKsikIK+bqmtZLWrkSpVJO4JZ+ujGF\ndF+KiTTF7AUntm/F657xcwfzi7LW7v5Vqsk8amOAj4/D3XdTvgzm57N1mdtisaRnuFJfN75tFos3\nL3q7YkHTb5ZOpyl5tPPEBOzfD6sV2Fo1xn9kGnP7Lcjq1tP14zg7uZ2J9TxFLslk87lrt2nMU0XZ\nFKhr2jAQJtFJ1Ak/qnSZvMXPuGsSvN3r2rVq7tQKOiCZh1cKVr9c3b5ufNVyLyO6YtEJf1sjovl8\nnkSiPoxpIuFkl62YGiEH2HdDnZADVNjKAjfUV6JqsOeX+LsXceUVRRlpVMx7SaMgePlL9+rB3KXL\nUD6V8g7dCd7ude1aNXdqBR2QzMMvJ7dXua8bX215iMQpmUyGP/3TG5mYEEScEfktt8Cb35ykyXX+\nmR3e9aOmvNoG1XP7xTsf4rgBiqL0HxXzXlMrCH7Tx714MPuIpX3FFaTTaWKxGOl0ej0XtW3bdeVc\ndRWFsTEsnKl1CygATXJW7YwEuYNV/bVFnNf27c53I4pHv/ty2PJHcPS495KRV2rW/HSe5Fh9+yTH\nkuSnPToTtk356u3M3SccPCDMPbKdcvn0bMru3Xt5/PHP89/+m8V99wlvfauTiWxromFkveMZ7/rF\nn/dvA/XjVhQlBLpm3k+6XMdum4asW/YVV5C99966nNTJZJJrr72Wez3KC9deS+bBB09n7fJb/w9a\n57ZtuO46OHmyvtxd525bvBsywO2+HO56IyAwfS7csrN+attvzRwcP/DcbI6lY0tMbpskP50nc0HD\nfrZN2b6O+ZtOUqmZIY9Vxtn5y3c3H7emzctXncN89udUYu6ygNeaeTLGzsJOJjI+iWx6nVFPUZSh\nJuyauYp5Pxnwg9nPWCsej7O2ttZUvp7a8vQB/AXdLw53q++024lp6Bxs+SNYqxHv6XPhhinYkYCX\nbLVCpWANJJ1m7rYSq+c1b0okLC66qKb+Hr9v+fIxFn7vbI5vWWb5RJxvfOESfv3+G0n9LMXWya1M\n5af8hbz2uMOSClVRlL6iYj6sDPDBHIvFaOf3FhEqtSNuPwvvKl4dk1ZJTzrJJ17ThvLHxlkH8MB8\nJIJ7Oxbj4CMGYvDII7BvHzzzDOzY4Ri5/cmf1Jyj3zMviqJseNSafVgJYVTVKyZ91ln90pg07V+7\nLu6FlzFf0Nqu17YwyWhq2jAe8659XOrLy3aZufQcB2MHmUvPUbZDxHu3bYjFSDzjCPntt0O57PRN\nymW4/XZZtzkAuvYgUBRF6RQV801EPp8n2WAUl8RJ49lkuZ5MkveyLG/Xyjqfh7Gx5v3Gx5st1ztI\nRpN9Q9a/3O0YlOUy5t/3PVZLq2BgtbTKfHY+WNCrdVlbY2of7PtzaMzSurpq6jOkqbGaoigDQsV8\nE5HJZCgUCliWVWehvtd9X7dct6zWObn9BOqcc+pH1gD33OMEW6mSSnkbv+VybSej2fvWvczsmlkf\nicclzsyuGf7nJUh/533EPlDiG2fdQMUk6r5XWamwkFvwv76aukzMOlPrXtRlSPPyIBBxEqQoiqL0\nEF0z36x0u77rtX4+NuaIV21Ql3YM/DqNCldbrSM2ex7aw/LK8vpa+uyts8S8+q0Cl1QuCVWXNOBl\nxlc1ErRtm1wux1KpxCSOP/76Fav1uaIoHaJr5kow3aau9PIrP/vs5uhs7QTF6XKa2j5ik30gy/I/\nL9cZxT2zzce/ezLhWe51zjz+SxG2bZPNZimVShgc0c8C64sDGrFNUZQeo2K+WYki53ejMd/zz3vv\nF9YArMsORm42x8rJZkv7fdP7OD5Wnyc9lowxlZ8KXZcMDUsRqdT6UkQul6vz0QdYAerk26sNwhj7\nKYqihECn2ZXoiMI1qwvXvdhHYxjP9DAw/cQ0Nz58A9t/PkHCasO/+/3v957iT6WcOO34u/wJsP7N\nxjbQYDCKooRAp9kHwbCNtPpdn26n7qEr1z2/eOsAc+c/yuon/ppLzKVctHhRayGv1sVvrb4mCp2f\ny996qVcbdGDspyiK4oeKeVR4uVW9732we3fwd3olttVUobX1uf76SM9hH7FJ35Em9tEY6TvS2BcS\nWbz1TvCKtw6QekmKwlWfJzOzt+Nj2zhGcDH3vbYVPV3+RJykNH5toD7piqJESZik58PyesMb3tBt\nnvfeYVnGOLJZ/xIxplhs3r9YNCaZrN83mfTetxNSKe/6pFKRHL74RNEk80nDray/kvmkKT4RUf27\nqJf1ScvIrWKsT1rd1yeVMkUwSTDUvJJgijW/VbFYNJZlGRExlmXVbfPE736xrO7qqyjKhgJ4zITQ\nR10zj4qgsKVea8a9Dv3pF9QF/OvZBuk70pSONdff2maxePNi8JdHKda4bZO+5ppAt7ROj6tr5oqi\ntELXzPtNkPuU19TpqE6zuulMl37mnTxl6ViL+ncQ5W2gZDL4XdFSN79VFN4EiqIoLirmUZHP+4+G\nvYQ+wtCf5bLN3FyagwdjzM2lnVzbtRHXavErD0N1HX55mclj3rsEGaEBgzH86tI2YdInFr2f4Vto\nBhinX1GUjYWKeVRkMnDjjc2C7mfNHYXlN46Qz//gelZXS4BhdbXE/A+up1y4qjkm+tgY3HlnW8ev\nI5dbDwqTn4VkQ3yY5FiS/HSL+vuNZkulrowAm4zxjrjHiWAmwNPAzS92vaIoygBQMY+SvXvh858P\nN3Ua0TTrwpE9VGL1qlqJnWDh7C86MdFrj3/PPd2N/mqEOHMECg+A9TMQ47wXriyQuaDF8YNGs62E\n1meEXY38VjpWwmAoHSuRfSDrCHoEMwF1Me1FwsWuVxRF6SNqADfs2Dbs2XParzmVckbXrpAcPCDe\n+bwNXHJpxL+tn9EedBfTPcyxAgzG0s/m/I3xfm+p43jv6/HWl5aYnJwkn8+rgCuK0lfUAG4jULNG\nvc7yMlx33fqoNOGTxdOvvCvyeSd1aSNjY53FdPfDayo+YITtZ3S3dGypY9uEunjrxlAqlchms/X5\nyxVFUYYEFfNBEcYoq2aNuo6TJ9eniae+kiJWH3ac2HGnPHIyGSd1aWM603an76uGX36C7iW0Adb/\nfkZ3k9smO7ZN8Iy3vrJSn79cURRlSFAxHwRhjbKCXJ/cbRNX3snO/zRG4ihQgcRR2Pmfxpi4sgtD\ntyAyGScmeTXMyXPPdb4O347QBoywvSK/rRvjdWib4Od21pU7mqIoSq8IE1lmWF5DHQGuHfyif8Xj\n9RHg/PZrjBRWLDqfRZz3qKLIdUKxWB99LpUKrk/YureImBd15DfLsuoivlVflkZoUxSljxAyAtzA\nBbqd1yiI+dGjRXP4sGUOHBBz+LBljh71EBURf5GuDelaLBozPt68TyzmiGSvxLvTzoFffcfGoqlj\nHzstxWLRJJPJ+hCuyWTrMK2KoigRomI+AI4eLZpDh5LmwAHWX4cOJZsFPWjE7TXqrh3pnnFGs2CG\njekeRgy7iRkf5rpGSAzbjreuKIoSMWHFXF3TImRuLu0Gb6knkbC46KLF0wWt3LOC3KY6jekeNhZ4\nNzHjg+LTB51TURRF8URd0wbA6qq3cVRTedUoKx73PlC7cd6DyquEDZ7STcz4MOFNow7dOmw55BVF\nUQaAinmEJBLeYuZZnsnAvfe27zblJ5ixWLCghRXpdv2ya8X0hRdgyxbv/cLUpV1GLWmLoihKj1Ax\nj5CpqTyxWL04x2JJpqZ8xDmM21TjyPOKK5o7AABra8GCFlak23EXaxTT5WXnOs480/tcrerSLoNI\n2qIoijKMhFlYH5bXsBvAGRPSmj0sfsZoMzOnDdni8dZGdEHH8jOCC2M17mfwVj13N8Z0YfDzChCJ\n5viKoigDBrVm3wC0Ektj2hO0qF27wpzb75xR1CVM+3igVuqKoowKKuYbgTBi2aGgRUIn5250tetm\nxN7ByF/9xxVFGSXCirmumQ8zYda5I8qL3hHtnru6xl6bOKZKJ2vdHYRq1ZjriqJsRNTPfJjx8A23\n3zDGA78zztutF5lIwFo8xQXPXcXE7z3oWIlPuslF+uXHbbs5w8OcOyiFKoRKS9otsVgMr3teRKj0\n+NyKoijtEqmfuYi8RUTmReQpEfmwx3ZLRGZF5AkROSgir6zZ9gkR+b77ek9NuYhIXkR+JCJPishN\nYS9u09Aw8rQvSWF/sMIHXvMi5211tG9LZZkfnPNZyv897wjh4mJ/A7JUM6CFOXcrl7SorNwDfM8n\nfc7hV64oijIKtBRzEYkDnwYuB84H3isi5zfsdjvwF8aYC4GPAR93v/tW4FeA1wO/Bvy+iJztfucD\nwKuAf2GMeR1wX9dXsxGpEcvc287k/a9eY2tDrJkYJ1hYCDlNPMggK0GCGdXSQAvf83w+T7JhaSCZ\nTJLvx7KEoihKjwgzMn8j8JQxZsEYcwJHdN/WsM/5wKz794Ga7ecDh4wxp4wxLwKPA29xt80AHzPG\nVACMMc90fhmbg6VjS+xIeG/ziz5Xh21Ttq9j7rYSBx8xzN1Womxf11tBbwwqMz7evE8qFV2I1wbf\n8/I0zH12hYMvv4a5uTSXXQaFQgHLshARLMuiUCiQ0fCyiqKMMGHE/BXAT2o+P+2W1fI48C7373cA\nZ4lIyi2/XESSIrIduBRnNA7wi8B7ROQxEXlIRH6p04vYLExum+SZVe9tftHnaik/sIf5m06yeh4Q\ng9XzYP6mk5Qf2BNtRat4BZUxxhHvqsFasdhdTvRGaqbyy9Mwf4tznQisrpaYn89y2WWwuLhIpVJh\ncXFRhVxRlJEnjJiLR1mjBdEtwMUi8l3gYuDvgVPGmK8DDwKHgf3AHHDK/U4COO4u7P85cLfnyUWy\nruA/9uyzz4ao7sYlP53nL5bGOL5WX15hvCnKnH3EJn1HmthHY6TvSGMfsVl4+zKVrQ3f3QoLb/ew\nLo8CrwhtJ086EeJ6tb5fM5W/cAPN11tZCb8koSiKMiKEEfOnOT2aBngl8NPaHYwxPzXGvNMY86+A\nnFt2zH3PG2Neb4z5DZyOwY9rjvtl9+//AlzodXJjTMEYs8sYs+vcc88NeVkbk8wFGTIX3cPnfpLi\n6HFnkHsqluKXX3c3ExOnRdE+YpN9IEvpWAmDoXSsRPaBLMcnvI+7uoPerKN3k7Sl07X9Gne51R3e\nu4RaklAURRkhQmTF4NvAL4nIq3FG3FcDv127gzuF/ry7/v0HuKNs13jupcaYZRG5EEewv+5+7SvA\nm9x9LwZ+1P3ljD72EZvcbI6lY0tMbpskP50nc8Fpoc5ckKn77EVuNsfKyQZf6pMrLK/G2L612f0q\n8SynXcaqBmPQ/ah5ctLbFa2V5XijS147dapuz+VIPFNyptgbCLMkoSiKMkq0HJkbY04Bvwv8FfAk\n8EVjzA9E5GMi8lvubpcA8yLyI2ACqM75jgHfEJEfAgXgGvd4ALcB7xKRIzjW7zdEdE0ji9+I2j7S\n3kh56Zj3yPMzCxVilXoDtNhxmPrzhh1XVmDPnu5H650GtOk2gYrrATB1abG9xDeKoigjigaNGSLS\nd6QpHWseyVrbLBZvXgx9nO3/YTvL/9y8Dp56SYofXHsnCws5VleXSCQmmfqjEhOzHgdpJJnszOK8\nnaAyVWIxZw2hkQ6CypTLdv31TuXrliQURVGGmbBBY1TMh4jYR2OYJttCEITKR8KL2Jn/95m8ePLF\npvIzxs7ghT98ob6wVVS2WizLMVrrNX516tf5FUVRhoRII8Ap/WFym090Mp9yP7yE3LfcayrcjzCG\na1EwyHjziqIoI4iK+RCRn86THGuITjaWJD/dQxHzSlaSSjXtZgPpWIxYLEY6ncbuZaCZDhKoKIqi\nbGZ0mn3IaGXNHga/NfOYOElGWh63wZrcBrJArUlaMpnUyGmKoig9RtfMNzH2EZvrv3o9J9ZO+O6T\nHEtSuLLQLOhVg7VSCeJxWFsjHY9TWltrOoZlWSzqGraiKErP0DXzAWLbNul0uj9T0h5kLshw99vu\nxtpmIQhxiTfts3Jyhdxsg6tXbfhVgLU1SCZZ8hBygKVSCbZsgd27o74ERVEUpQ10ZB4xtm2TzWZZ\nqfGTHvSUdGgreR8rct+RObBY/TAzA3v3RlFdRVEUxUVH5gMil8vVCTnAysoKubABT3pAaCt5H2v1\n/Npac9pQTkcGAhwDtR5SLtvMzaU5eDDG3Fyacrm/sx2KoijDjIp5xCz5CKJfeT9oZSW/npTljw3p\nm8G+oP77GTdNqGVZCM6IvADUzTP4TMVHQblsMz+fZXW1BJj17Gcq6IqiKA4q5hEz6RN33K+8H2Qu\nyFC4srC+hm5ts9aN3+pCyAqUXgrZK2sE3fXvzmQyTtrQeJxFGoQcHGO5HrGwkKNSqZ/t0OxniqIo\np1Exj5h8Pt88JZ1Mkh9wwJPMBRkWb16k8pEKizcvrluxeyZlGYfcNN7+3dWEJ41Uy20btm93/MNF\nnL+7NAD0y3Km2c8URVEcVMwjJpPJnJ6SFsFyp6iH1R/bLylL6aUQu26J9LO5+kQve/c6xm7VkXg8\nftr4zbbh+uthucbHfXkZrruuK0H3y3Km2c8URVEc1Jp9k+OX3KUWX5/0poOl/eO8dxFXvbpmXjvV\nHosl2bmzoElTFEXZ0Kg1uxIKL+O4RlZOrpC7f0/rgwUZ+XVhADgxkWHnzgKJhAUIiYSlQq4oilLD\nlkFXQBkstWvnQSP0pZPLzlR50HLB5KT/yLxLA8CJiYyKt6Ioig86MlfWjeME8d1n8hhOmNcg8nkY\nH28uHxtzttm2MxUfiznvfY6MpyiKslFRMVfWOecl53iWi4H8LK2nyjMZuPvu+qxrqRTcc4/zdzVU\nrDHOezZbL+gq9oqiKB2hYq4ATuCYn5/4efMGAzd+CzJHaDlVbts26VyO2PPPk7Ys7GIRnnvOEflc\nbj0L2zorK6dH+7Vx4f3EXlEURfFExVwBnDVzryxrqRXY+xDrwWP8qMakL5VKGGMolUpks9nTSWb8\nRvWl0ulMbUFiryiKoviirmkKEJCMxUDlc5Yj5AHGb+l0mpKH8dt6mtQgt7VkslnI1ysgUKl4b1MU\nRdngqGua0ha+yVhe6vqHtwh60zImfT7viLYXKyv+4WAHGAZXURRlVFAxV4DWyVha0TImfSYTnFnN\nzZ1eX4HgqX1FURTFQcVcAYKTsQTiWqDnSyWSUu/a1hSTPpNxIsF5UY0Db1nO1LpXXHhFURTFE10z\nVzqnaoHurnfbQE6EJWOYtCzybra1oO8AzghchVtRFKWJsGvmGgFO6ZwGC/QMkDEmOA57VbBzOcfC\nfXKypXGdoiiKEoyOzJXOicUcn/BG1AJdURQlEtSaXek9fpbmaoGuKIrSV1TMlc7xcjdTC3RFUZS+\no2KudE7V3Uwt0BVFUQaKGsAp3ZHJqHgriqIMGB2ZK4qiKMqIo2KuKIqiKCOOivlGZPdu2LLFWcfe\nssX5rCiKomxYVMw3ErYNZ54Jd93lxDoH5/2uu/or6G6IV2Ix511zkiuKovQUFfONgm3D9dfDiy96\nbw9KchJ1PbJZJ92pMc57NquCriiK0kM0AtxGYft2WF4O3qcfv7Vf3vKgEK+KoiiKJxoBbrPRSsj9\n8oVHjU8kTXwJAAAMgklEQVRec99yRVEUpWtUzDcL2Wx/zqMhXhVFUfqOivlGIZXy3zYzA3v3NpdH\naahWPVap5FjR16IhXhVFUXqKivlG4c47YWysvmxsDIpFfyGPylCt9ljgHK8q6BriVVEUpeeoAdxG\nwrbD5wmP0lBNjd4URVF6QlgDOBXzzUqUucg1r7miKEpPUGt2xRPbtkmn08SMIQ00Tap3YqimRm+K\noigDRcV8E2HbNtlsllKphAFKQJYaQe/UUE3zmiuKogyUUGIuIm8RkXkReUpEPuyx3RKRWRF5QkQO\nisgra7Z9QkS+777e4/HdT4nIC91dhhKGXC7HyspKXdkKkIPuDNU0r7miKMpAaZnPXETiwKeB3wCe\nBr4tIvcbY35Ys9vtwF8YY+4VkTcBHwfeJyJvBX4FeD2QAA6JyEPGmH9yj70LeGmkV6T4suQTuGVJ\npHtDNc1rriiKMjDCjMzfCDxljFkwxpwA7gPe1rDP+cCs+/eBmu3nA4eMMaeMMS8CjwNvgfVOwn8E\n/o/uLkEJy6TPGrZfuaIoijIahBHzVwA/qfn8tFtWy+PAu9y/3wGcJSIpt/xyEUmKyHbgUuBV7n6/\nC9xvjPmHoJOLSFZEHhORx5599tkQ1VX8yOfzJBvWtpPJJHld21YURRlpwoi5eJQ1+iHdAlwsIt8F\nLgb+HjhljPk68CBwGNgPzAGnROTlwL8GPtXq5MaYgjFmlzFm17nnnhuiuhuIiPOSZzIZCoUClmUh\nIliWRaFQIFMzPb5u7R6LkU6nsTXbmaIoytDT0s9cRC4CbjXG/Kb7+Q8AjDEf99n/TOBvjDGv9Nj2\nBaCI00H4LHDc3TQJLBhjXhNUl03lZ757t5OHvBG/0KwRULV2rzWSSyaTTYKvKIqi9IfIgsaIyBbg\nR8A0zoj728BvG2N+ULPPduB5Y0xFRPLAmjHmj9118ZcaY5ZF5ELgC8DrjTGnGs7xgjHmzFaV3VRi\nvmULrK01l8fjcOpUc3kEpNNpSh6R3CzLYlEjuSmKovSdsGLe0prdGHNKRH4X+CsgDtxtjPmBiHwM\neMwYcz9wCfBxETHAXwP/u/v1MeAb4sTp/ifgmkYhV3zwEvKg8gjwtXbX9KWKoihDjYZzHVaGaWQe\nj7N4773qeqYoitJnNJzrqOOXf7yHeck9rd2B/Npa5xnVFEVRlJ6jYj6s7N3rGLvF487neLynxm9Q\nY+0ejyOABRSADMDKipORTVEURRk6dJpdaUazoCmKogwFOs2udI5mQVMURRkpVMwHjW1DOu2MhtPp\n3q5L155r+3bn5XVezYKmKIoyUrR0TVN6iG07hmXVIC2l0mkDt6gtxxvPtbx8elvjeavnzuVgackZ\nkefzas2uKIoypOia+SBJpx0hbcSyus9iFvZcvT6voiiK0jG6Zj4K+AVj6UWQljDH1OAwiqIoI4mK\n+SDpp6FZmGOqgZuiKMpIomI+SPppaOZ1rn6cV1EURek5KuZ9oFy2mZtLc/BgjLm5NOWyazmeyUCh\n4KxVizjvhUJvDM0az5VKOa9en1dRFEXpOWoA12PKZZv5+SyVyum0orFYkp07C0xMhBBP21arckVR\nlE2KGsANCQsLuTohB6hUVlhYCBEatepOVio5EdmqLmQaI11RFEWpQcW8x6yueluI+5XXkcud9guv\nojHSFUVRlAZUzHtMIuFtIe5XXkc/XdcURVGUkUXFvMdMTeWJxeqtyGOxJFNTISzHNUa6oiiKEgIV\n8x4zMZFh584CiYQFCImEFd74TWOkK4qiKCHQ2Ox9YGIiE068G9EY6YqiKEoIdGQ+7GQyTrz0SsV5\nrxXyfmZcUxRFUYYWHZmPKv3MuKYoiqIMNToyH1XUbU1RFEVxUTEfVdRtTVEURXFRMR9VfNzTyled\n4x0HXlEURdmwqJiPKh5ua+XLx/ibf3OM1dUSYFhdLfHkk9fwve9dNpg6KoqiKH1BxXxU8ci49uN/\nl8DIqaZdf/azWX70o90DqKSiKIrSD1TMR5kGt7VT8Rd8d/3pTwv9q5eiKIrSV1TMNw1rg66AoiiK\n0iNUzDcQ8XgqaGvf6qEoiqL0FxXzDcRrX3snIJ7bXv7ybH8royiKovQNFfMNxMREhte97vOInFFT\nGuPlL5/hta/dO7B6KYqiKL1Fw7luMDpO6qIoiqKMLDoyVxRFUZQRR8VcURRFUUYcFXNFURRFGXFU\nzBVFURRlxFExVxRFUZQRR8VcURRFUUYcFXNFURRFGXFUzBVFURRlxFExVxRFUZQRR8VcURRFUUYc\nFXNFURRFGXFUzBVFURRlxFExVxRFUZQRR8VcURRFUUYcFXNFURRFGXHEGDPoOoRGRJ4FSgOuxnbg\nuQHXYRTRdusMbbfO0bbrDG23zulF21nGmHNb7TRSYj4MiMhjxphdg67HqKHt1hnabp2jbdcZ2m6d\nM8i202l2RVEURRlxVMwVRVEUZcRRMW+fwqArMKJou3WGtlvnaNt1hrZb5wys7XTNXFEURVFGHB2Z\nK4qiKMqIo2IOiMi/FpEfiEhFRHwtEUXkLSIyLyJPiciHa8pfLSL/XUR+LCJ/KSLjDd97t4iYoGOP\nKr1qOxH5PRH5oYg8ISKzImL143r6RQ/bLeF+fsrdnu791fQPETlHRB52r/thEXmZz36fEJHvu6/3\n1JS/SUT+h1t+r4hsccu3icgDIvK4+7tc169r6he9ajt32yUi8j237Q7143r6RS/bzd3+qyKyJiLv\n7qqixphN/wJeB+wEDgK7fPaJA38LTAHjwOPA+e62LwJXu39/Bpip+d5ZwF8D3/Q79ii/etV2wKVA\n0v17BvjLQV/riLTbbuAz7t9Xb8B2+w/Ah92/Pwx8wmOftwIPA1uAM4DHgLNxBi8/AV7r7vcx4IPu\n339YPRZwLvA8MD7o6x2Rtnsp8ENg0v28Y9DXOgrt5n6OA48CDwLv7qaeOjIHjDFPGmPmW+z2RuAp\nY8yCMeYEcB/wNhER4E3Af3b3uxd4e833/i+cm+F4xNUeCnrVdsaYA8aYFbf8m8Aro6/94OjhPfc2\n9zPu9ml3/41C7fU1/q9VOR84ZIw5ZYx5EacT9BYgBawaY37k7vcw8C73bwOc5bbVmThifqo3lzAw\netV2vw38v8aYJQBjzDM9qv+g6FW7Afxb4MtA122mYh6eV+D0sKo87ZalgJ8ZY041lCMi/wp4lTHm\n/+tnRYeQttuugQ8CD/W0hsNJJ+22/h13+zF3/43ChDHmHwDc9x0e+zwOXC4iSRHZjjPL8yqcyFxj\nNcsa73bLAf4fnNmSnwJHgD3GmErvLmMg9KrtXgu8TEQOish3ROT9Pb2K/tOTdhORVwDvwJlZ65ot\nrXfZGIjII8B5HptyxpivhjmER5nxKxeRGPBJ4AOhKzmk9LvtGs59DbALuDjEeYaKAbVbyzYddoLa\nLcz3jTFfF5FfBQ4DzwJzwCljjBGRq4FPikgC+DqnR9+/CXwPZ8bjF4GHReQbxph/6u5q+suA2m4L\n8AZgGngJMCci36wZjQ49A2q3O4B/b4xZi2LybNOIuTHmsi4P8TSne6LgTPv+FKfn9VIR2eKOhKrl\nZwH/Ejjo/lDnAfeLyG8ZYx7rsi59ZQBtB4CIXIbzz3SxMWa1yzr0nQG1W/U7T7uGNttwpoxHhqB2\nE5GyiPyCMeYfROQX8JmeNMbkgbz7nS8AP3bL54Bfd8vfjDOqBLgOuM04C5lPicjfAf8C+FY0V9Uf\nBtR2TwPPudPLL4rIXwP/EzAyYj6gdtsF3Ofqw3bgChE5ZYz5SifXoNPs4fk28EviWBGP4xgX3e/+\n8x/AmT4BuBb4qjHmmDFmuzEmbYxJ46z7jpyQR0RbbQfrSxR/htNmG20NLixttxtwv/sZd/uj7v4b\nhdrrq73udUQkLiIp9+8LgQtxRkSIyA73PQH8e05PcS7hjCwRkQkc48SFnl3FYOhV230V+HUR2SIi\nSeDXgCd7eB39piftZox5dY0+/Gdgd6dCjnvATf/CWbd4GlgFysBfueUvBx6s2e8KnN7m3+JMlVbL\np3B68E8BXwISHuc4yMa0Zu9J2wGPuMf7nvu6f9DXOiLtttX9/JS7fWrQ1xpxu6WAWZxRzyxwjlu+\nC9hX0wY/dF/fBF5f8/3/iCM088DNNeUvx3n4HgG+D1wz6GsdlbZzt/2++53vN24b9Vcv261mn8/R\npTW7RoBTFEVRlBFHp9kVRVEUZcRRMVcURVGUEUfFXFEURVFGHBVzRVEURRlxVMwVRVEUZcRRMVcU\nRVGUEUfFXFEURVFGHBVzRVEURRlx/n/JfDP+qCmi+AAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "k = 7\n", "km = KMeans(n_clusters=k)\n", "km.fit(coords)\n", "labels = km.labels_\n", "centroids = km.cluster_centers_\n", "colors = ['r', 'g', 'b', 'y', 'c', 'm', 'k']\n", "plt.figure(figsize=(8,8))\n", "plt.title('Clustered Locations')\n", "for j in range(k):\n", " current_coords = np.asarray(coords)[np.where(labels==j)]\n", " j_longs = current_coords[:,1]\n", " j_lats = current_coords[:,2]\n", " plt.plot(j_longs, j_lats, '%so' % colors[j], label=j)\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "7 is interesting indeed. Mississauga and Toronto are separated, and all but one of the Costco stations are grouped together regardless of distance." ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "collapsed": true }, "outputs": [], "source": [ "cluster_df['label'] = km.labels_" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "label\n", "2 122.913025\n", "0 122.780663\n", "6 120.150000\n", "3 118.368056\n", "1 116.464407\n", "4 113.891667\n", "5 111.114286\n", "Name: mean, dtype: float64" ] }, "execution_count": 68, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cluster_df.groupby('label')['type_A'].describe().sort_values('mean',ascending=False)['mean']" ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Alongitudelatitudetype_A_normlongitude_normlatitude_normlabel
1411 Warden Ave , Toronto - East, OntarioCostco110.9-79.29773643.7594290.917363-0.9976281.0006125
150 Kingston Rd E , Ajax, OntarioCostco110.9-79.01845443.8629500.917363-0.9941141.0029795
1570 Dundas St E , Mississauga, OntarioCostco110.9-79.57577643.6103040.917363-1.0011260.9972025
55 New Huntington Rd , Vaughan, OntarioCostco110.9-79.63918943.7699370.917363-1.0019241.0008525
71 Colossus Dr , Vaughan, OntarioCostco110.9-79.54159043.7856310.917363-1.0006961.0012115
935 Liverpool Rd , Pickering, OntarioShell112.4-79.08691243.8285480.929771-0.9949761.0021925
Taymall Ave , Toronto - West, OntarioCostco110.9-79.50661643.6238870.917363-1.0002560.9975125
\n", "
" ], "text/plain": [ " name type_A longitude \\\n", "1411 Warden Ave , Toronto - East, Ontario Costco 110.9 -79.297736 \n", "150 Kingston Rd E , Ajax, Ontario Costco 110.9 -79.018454 \n", "1570 Dundas St E , Mississauga, Ontario Costco 110.9 -79.575776 \n", "55 New Huntington Rd , Vaughan, Ontario Costco 110.9 -79.639189 \n", "71 Colossus Dr , Vaughan, Ontario Costco 110.9 -79.541590 \n", "935 Liverpool Rd , Pickering, Ontario Shell 112.4 -79.086912 \n", "Taymall Ave , Toronto - West, Ontario Costco 110.9 -79.506616 \n", "\n", " latitude type_A_norm \\\n", "1411 Warden Ave , Toronto - East, Ontario 43.759429 0.917363 \n", "150 Kingston Rd E , Ajax, Ontario 43.862950 0.917363 \n", "1570 Dundas St E , Mississauga, Ontario 43.610304 0.917363 \n", "55 New Huntington Rd , Vaughan, Ontario 43.769937 0.917363 \n", "71 Colossus Dr , Vaughan, Ontario 43.785631 0.917363 \n", "935 Liverpool Rd , Pickering, Ontario 43.828548 0.929771 \n", "Taymall Ave , Toronto - West, Ontario 43.623887 0.917363 \n", "\n", " longitude_norm latitude_norm \\\n", "1411 Warden Ave , Toronto - East, Ontario -0.997628 1.000612 \n", "150 Kingston Rd E , Ajax, Ontario -0.994114 1.002979 \n", "1570 Dundas St E , Mississauga, Ontario -1.001126 0.997202 \n", "55 New Huntington Rd , Vaughan, Ontario -1.001924 1.000852 \n", "71 Colossus Dr , Vaughan, Ontario -1.000696 1.001211 \n", "935 Liverpool Rd , Pickering, Ontario -0.994976 1.002192 \n", "Taymall Ave , Toronto - West, Ontario -1.000256 0.997512 \n", "\n", " label \n", "1411 Warden Ave , Toronto - East, Ontario 5 \n", "150 Kingston Rd E , Ajax, Ontario 5 \n", "1570 Dundas St E , Mississauga, Ontario 5 \n", "55 New Huntington Rd , Vaughan, Ontario 5 \n", "71 Colossus Dr , Vaughan, Ontario 5 \n", "935 Liverpool Rd , Pickering, Ontario 5 \n", "Taymall Ave , Toronto - West, Ontario 5 " ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cluster_df[cluster_df['label']==5]" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### Folium\n", "\n", "Let's put it on a Toronto map with Folium.\n", "\n", "Then run our function for getting the nearest gas stations within a specified distance in kilometers." ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import branca\n", "import folium\n", "from folium.plugins import MarkerCluster" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nametype_Alongitudelatitudetype_A_normlongitude_normlatitude_normlabel
12731 Hwy 48 , Whitchurch-Stouffville, OntarioUltramar115.6-79.28237843.9763460.956242-0.9974351.0055721
1610 Keele St , Toronto - Central, OntarioShell122.9-79.47198843.6823201.016627-0.9998200.9988490
6897 Finch Ave W , Toronto - West, OntarioEsso121.9-79.61736443.7346321.008355-1.0016491.0000450
1 Harwood Ave S , Ajax, OntarioPioneer122.6-79.02510043.8610271.014146-0.9941981.0029352
1 Thornhill Woods Dr , Vaughan, OntarioEsso122.9-79.46388643.8266601.016627-0.9997181.0021492
\n", "
" ], "text/plain": [ " name type_A \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario Ultramar 115.6 \n", " 1610 Keele St , Toronto - Central, Ontario Shell 122.9 \n", " 6897 Finch Ave W , Toronto - West, Ontario Esso 121.9 \n", "1 Harwood Ave S , Ajax, Ontario Pioneer 122.6 \n", "1 Thornhill Woods Dr , Vaughan, Ontario Esso 122.9 \n", "\n", " longitude latitude \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -79.282378 43.976346 \n", " 1610 Keele St , Toronto - Central, Ontario -79.471988 43.682320 \n", " 6897 Finch Ave W , Toronto - West, Ontario -79.617364 43.734632 \n", "1 Harwood Ave S , Ajax, Ontario -79.025100 43.861027 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -79.463886 43.826660 \n", "\n", " type_A_norm \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 0.956242 \n", " 1610 Keele St , Toronto - Central, Ontario 1.016627 \n", " 6897 Finch Ave W , Toronto - West, Ontario 1.008355 \n", "1 Harwood Ave S , Ajax, Ontario 1.014146 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1.016627 \n", "\n", " longitude_norm \\\n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario -0.997435 \n", " 1610 Keele St , Toronto - Central, Ontario -0.999820 \n", " 6897 Finch Ave W , Toronto - West, Ontario -1.001649 \n", "1 Harwood Ave S , Ajax, Ontario -0.994198 \n", "1 Thornhill Woods Dr , Vaughan, Ontario -0.999718 \n", "\n", " latitude_norm label \n", " 12731 Hwy 48 , Whitchurch-Stouffville, Ontario 1.005572 1 \n", " 1610 Keele St , Toronto - Central, Ontario 0.998849 0 \n", " 6897 Finch Ave W , Toronto - West, Ontario 1.000045 0 \n", "1 Harwood Ave S , Ajax, Ontario 1.002935 2 \n", "1 Thornhill Woods Dr , Vaughan, Ontario 1.002149 2 " ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cluster_df.head()" ] }, { "cell_type": "code", "execution_count": 72, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "folmap(cluster_df,zoom=10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Enter the address and distance threshold to see nearby stations." ] }, { "cell_type": "code", "execution_count": 73, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Getting 80 Bloor St W, Toronto, ON M5S 2V1\n", "Got 80 Bloor St W, Toronto, ON M5S 2V1, Canada \n", "\n", "Gas stations within range: \n", " name type_A\n", "835 Yonge St , Toronto - Central, Ontario Canadian Tire 119.3\n", "1077 Yonge St , Toronto - Central, Ontario Shell 122.9\n", "150 Dupont St , Toronto - Central, Ontario Esso 122.9\n", "241 Church St , Toronto - South, Ontario Esso 122.9\n", "505 Jarvis St , Toronto - South, Ontario Petro-Canada 122.9\n", "581 Parliament St , Toronto - South, Ontario Esso 122.9\n" ] }, { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "enter_address(cluster_df,'80 Bloor St W, Toronto, ON M5S 2V1',2)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }