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)