Mapping Queensland's Fatal Crash Locations
Visualization of Queensland Road Crash Sites¶
Introduction¶
Background¶
Recently I attended (mostly as an observer) the 2019 GovHack event, which describes itself as the "THE LARGEST OPEN DATA HACKATHON IN THE SOUTHERN HEMISPHERE" (which is a pretty bold claim). I couldn't commit to the time needed to participate fully (which was why I was mostly observing), but found the process very interesting. The event was very well run at the Peregian Digital Hub, just north of where I live.
Govhack provides access to a large array of Government [datasets])(https://hackerspace.govhack.org/data_sets), including road crash locations for Queensland.
Goal¶
The road crash locations databse is one that I have been looking for, as I have a dream of a mobile phone app that would warnimg drivers about dangerous sections of the road, so that they can pay extra care to driving. This notebook is about the first stage of achieving that dream, in visualizing crash locations in various ways. I've decided to show the initial data exploration and display, as well as the more polished final results.
Implementation¶
This sections contains the housekeeping commands related to implementation details.
Notebook magic commands¶
Show all graphics inline in notebook.
%matplotlib inline
Watermark is used to show environmental details for reproducability
%load_ext watermark
I use black as my python code formatter, tailored for Pelican webpage production.
%load_ext lab_black
Notebook imports¶
Note that some imports are not actually used by the visualization code, but are there to support display of reproducability information.
# all imports should go here
import sys
import os
import subprocess
import datetime
import platform
import datetime
import seaborn as sns
import pandas as pd
import geopandas as gpd
from cartopy import crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.img_tiles import GoogleTiles
import cartopy.io.img_tiles as cimgt
import six
from PIL import Image
import matplotlib.pyplot as plt
Data Load¶
Initial Exploration¶
So the first thing to do is to load and explore the data (held in a downloaded csv file. We get a somewhat opaque warning message: it turns out that the auto-discovery of data column types can get confused by things like PostCodes that (in Australia in any case) mostly look like integers.
crash_path = 'D:\\QLDCrashLocations\\data\\locations.csv'
crashes = pd.read_csv(crash_path)
As usual, we do a .head() to see what we have.
crashes.head()
This all looks promising: lets see what types of crashes we have recorded.
crashes["Crash_Severity"].unique()
Lets see what the pesky column 15 is all about.
crashes.columns[15]
A little Googling on Stack Overflow reveals the source of the problem. We declare PostCode to be strings, and reload data (not really necessary, but just checking this fixes the warning message)
crash_path = 'D:\\QLDCrashLocations\\data\\locations.csv'
crashes = pd.read_csv(
crash_path, dtype={'Loc_Post_Code': str}
)
Initial Geodata Exploration¶
At this stage I decide that fatal crashes were the only ones I was interested in, so I create a Pandas view to the crashes dataframe. This is a touch unfeeling of me, in that I am using data representing human tragedy and suffering for a blog post to amuse my readers, but my long term goal is an app that might reduce this tragedy.
fatal = crashes[crashes["Crash_Severity"] == 'Fatal']
Let's check what this this gives us.
fatal.head()
Let's see how many fatal crashes we have. Turns out that we have quite a few.
len(fatal)
Now, let's explore the geographic data, and check that it looks reasonable.
fatal["Crash_Longitude_GDA94"].min()
fatal["Crash_Longitude_GDA94"].max()
The range of Longitude values looks reasonable for Queensland. If we look at the dataframe for all crashes, we see that there are some crashes with a suspicious (wrong) value. Clearly, some crashes have a missing Longtitude value, that gets coded as 0.0 in the csv file.
crashes["Crash_Longitude_GDA94"].min()
Let's see how many crashes have no Lat/Lon values.
len(crashes[crashes["Crash_Longitude_GDA94"] < 0.1])
Compared to the total numbers of crashes, the fraction with Lat/Lon missing values is quite small.
len(crashes[crashes["Crash_Longitude_GDA94"] > 0.1])
It turns out that all the fatal crashes have a recorded Lat/Lon pair.
len(fatal[fatal["Crash_Longitude_GDA94"] > 0.1])
At this stage, I decided to use GeoPandas, to support my preliminary visualizations. We show all the columns again, to refresh my memory of what data we have. I could have chosen to extract just the Lat/Lon pairs, but decided to keep all the data columns in my GeoPandas data structures: the numbers of rows is not huge, and I may choose to use this data somewhere in my visualizations.
fatal.columns
Creating a GeoPandas DataFrame¶
fatal_gdf = gpd.GeoDataFrame(
fatal,
geometry=gpd.points_from_xy(
fatal.Crash_Longitude_GDA94,
fatal.Crash_Latitude_GDA94,
),
)
The Lat/Lon pairs in the geometry column look reasonable.
fatal_gdf.geometry.head()
We now do the almost the most minimal visualization. It turns out that the Natural Earth world dataset doesn't think that Australia is a continent, so I had to use Oceania, which includes French Polynesia on the far side of the Date Line.
world = gpd.read_file(
gpd.datasets.get_path('naturalearth_lowres')
)
# We restrict to Oceania
ax = world[world.continent == 'Oceania'].plot(
color='white', edgecolor='black'
)
# We can now plot our ``GeoDataFrame``.
fatal_gdf.plot(ax=ax, color='red')
plt.show()
A listing of what Natural Earth thinks are the continents(?!).
world.continent.unique()
The absolutely minimal GeoPandas visualization. You can make out the shape of Queensland, and even what I suspect is Birdsville in the South West corner.
fatal_gdf.plot()
Define an imagery class for Stamen, should we decide to use it.
class StamenTerrain(GoogleTiles):
def _image_url(self, tile):
x, y, z = tile
url = 'http://tile.stamen.com/terrain/{}/{}/{}.png'.format(
z, x, y
)
return url
# end _image_url
# end StamenTerrain
We have to chose a Coordinate Reference System (CRS) for our maps. To make things easy, I have chosen PlateCarree: Queensland looks OK in this projection (unlike, say continential USA). We re-create the GeoPandas dataframe.
# define the CRS we will use for all maps
my_crs = ccrs.PlateCarree()
fatal_gdf = gpd.GeoDataFrame(
fatal,
geometry=gpd.points_from_xy(
fatal.Crash_Longitude_GDA94,
fatal.Crash_Latitude_GDA94,
),
)
Plotting Queensland Fatal Crashes¶
In this figure, I use the higher resolution coastline and state borders (1:10M) to provide geographic context. I manually set the extent of the map.
Each fatal crash is shown as a dot with a very low alpha (i.e. very translucent), so that areas with high fatal crash density are easier to perceive. Mount Isa is unexpectedly very prominent.
fig = plt.figure(figsize=(20, 20))
ax = fig.add_subplot(1, 1, 1, projection=my_crs)
ax.set_extent((137.4, 154.2, -30, -9.8), crs=my_crs)
ax.coastlines(resolution='10m')
# Create a feature for States/Admin 1 regions at 1:10m from Natural Earth
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='10m',
facecolor='none',
)
# plot state border
SOURCE = 'Natural Earth'
LICENSE = 'public domain'
ax.add_feature(states_provinces, edgecolor='gray')
# plot fatal crashes
fatal_gdf.plot(ax=ax, color='red', alpha=0.1)
# set title
ax.set_title('Highways of Blood', fontsize=20)
plt.show()
Showing Local Crashes¶
If we want to show the local area (for, say a Sunday drive), then we need more geographic context, showing roads and towns. The example below shows fatal crashes in the context of Google tiles. I have chosen to use the CRS of the image tiles for this map, as re-projecting images can result in a fuzzy looking basemap.
As before, the map extent is set manually to my local patch, and the zoom level is set manually to a value that looks good, and has the right amount of detail. On the recommendation of Stack Overflow, I have the interpolation='spline36' parameter on the add_image call, as it might improve image quality if I choose a different CRS.
fig = plt.figure(figsize=(20, 20))
imagery = cimgt.GoogleTiles()
ax = plt.axes(projection=imagery.crs)
ax.set_extent(
(153, 153.2, -26.6, -26.4), crs=ccrs.PlateCarree()
)
# Add the imagery to the map.
zoom = 12
ax.add_image(imagery, zoom, interpolation='spline36')
# Here's what the plot looks like in GeoPandas
crs_proj4 = imagery.crs.proj4_init
fatal_gdf.crs = {'init': 'epsg:4326'}
fatal_gdf2 = fatal_gdf.to_crs(crs_proj4)
fatal_gdf2.plot(ax=ax, color='red', alpha=0.8, zorder=4)
plt.title(
"Local Fatal Crashes: Sunshine Coast",
{"fontsize": 30},
pad=40,
)
plt.show()
The diagram above represents a mock-up of what might be my final app UI. If you use Google Maps to plan a trip, then the red dots representing fatal crashes might allow the user to avoid 'Black Spots', and show where extra care might be needed.
Conclusions / Summary¶
This is a work in progress: the nxet steps will be to try Folium, and get a fata crash webserver going.
Reproducibility¶
This section of the notebook provides information to support reproducability.
Notebook version status¶
theNotebook = (
'2019-10-21-dc-QLDCrashLocationsNotebook.ipynb'
)
# show info to support reproducibility
def python_env_name():
envs = subprocess.check_output(
'conda env list'
).splitlines()
# get unicode version of binary subprocess output
envu = [x.decode('ascii') for x in envs]
active_env = list(
filter(lambda s: '*' in str(s), envu)
)[0]
env_name = str(active_env).split()[0]
return env_name
# end python_env_name
print('python version : ' + sys.version)
print('python environment :', python_env_name())
print('pandas version : ' + pd.__version__)
print('current wkg dir: ' + os.getcwd())
print('Notebook name: ' + theNotebook)
print(
'Notebook run at: '
+ str(datetime.datetime.now())
+ ' local time'
)
print(
'Notebook run at: '
+ str(datetime.datetime.utcnow())
+ ' UTC'
)
print('Notebook run on: ' + platform.platform())
%watermark
%watermark -h -iv