1 Import libs
In [1]:
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import os
import pandas as pd
import sqlite3 as dbi
2 Get and clean the data
In [14]:
crashes = pd.read_csv('../data/crash_data.csv')
crashes = crashes.dropna()
# crashes = crashes[crashes['road_inv'] == '005']
3 Prepare the parameters
In [15]:
lons = crashes['longitude']
lats = crashes['latitude']
counts = crashes['tot_acc_ct']
rates = (counts * 1000000)/ (crashes['avg_aadt'] * 365 * 6 * crashes['seg_lng'])
In [16]:
lons = lons.tolist()
lats = lats.tolist()
counts = counts.tolist()
rates = rates.tolist()
4 Prepare the hot map plotting tool
In [17]:
def draw_crash_map(shpurl, name, llon, llat, rlon, rlat, lons, lats, data):
'''
A function to plot hot spot map for crash rate and
crash severity
Parameters:
@shpurl {string} the shapefile url, without suffix,
e.g. '../data/highway/wgs84'
@name {float} either 'rate' or 'severity'
@llon {float} longitude of the top left corner
@llat {float} latitude of the top left corner
@rlon {float} longitude of the bottom right corner
@rlat {float} latitude of the bottom right corner
@lons {a list of float} longitudes of the crash spots
@lats {a list of float} latitudes of the crash spots
@data {a list of float} data for crash rates or severity
Return: Nothing
'''
# set up a map canvas
map = Basemap(llcrnrlon=llon, llcrnrlat=llat,
urcrnrlon=rlon, urcrnrlat=rlat, resolution='i',
projection='tmerc', lat_0=(llat+rlat)/2,
lon_0=(llon+rlon)/2)
# map.drawmapboundary(fill_color='#9999FF')
# map.fillcontinents(color='#ddaa66',lake_color='#9999FF')
# map.drawstates(color='0.5')
# read the shapefile
map.readshapefile(shpurl, 'highway')
# get the max value from the data to scale the size of the marker
max_val = max(data)
min_val = 10000
# plot the hot spots one by one with the marker
# size corresponding to the data
for index in range(len(data)):
x, y = map(lons[index], lats[index])
map.plot(x, y, marker='o', color='r', markersize=((data[index]*15/max_val)+5))
if data[index] > 0 and data[index] < min_val:
min_val = data[index]
red_dot1, = plt.plot([], "ro", markersize=20)
red_dot2, = plt.plot([], "ro", markersize=5)
str1 = 'Max crash rate:' + str(round(max_val,2))
str2 = 'Min crash rate:' + str(round(min_val,2))
plt.legend([red_dot1, red_dot2], [str1, str2])
plt.title('WA Highway Crash Rate Hot Spot Map')
# show the map
plt.show()
5 Invoke the tool
In [18]:
draw_crash_map('../data/highway/wgs84', 'rate', -125.0, 49.50, -116.5, 45.0, lons,
lats, rates)