notebooks/01_explore_cas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploration of CAS (Crash analysis system) data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import seaborn as sb\n",
    "import contextily as cx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from crash_prediction import cas_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set seaborn default style\n",
    "sb.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data retrieval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need to retrieve the dataset from the [Open Data portal](https://opendata-nzta.opendata.arcgis.com/datasets/crash-analysis-system-cas-data-1).\n",
    "Multiple file formats are available (csv, kml, geojson, ...), the most compact\n",
    "being the .csv one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dset_path = Path(\"..\") / \"data\" / \"cas_dataset.csv\"\n",
    "if not dset_path.exists():\n",
    "    dset_path.parent.mkdir(parents=True, exist_ok=True)\n",
    "    cas_data.download(dset_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we load the data and have a quick look to check if there no obvious\n",
    "loading error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dset = pd.read_csv(dset_path)\n",
    "dset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset contains 72 columns, describing various aspects of the recorded\n",
    "car crashes. The full description of the fields is available online, see\n",
    "https://opendata-nzta.opendata.arcgis.com/pages/cas-data-field-descriptions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dset.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `X` and `Y` are geographical coordinates using the WGS84 coordinate\n",
    "system (see [EPSG:4326](https://epsg.io/4326))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Spatial features\n",
    "\n",
    "First, we will look at the location of the crashes. More accidents happen in\n",
    "densier areas and it would be good to compare with population density.\n",
    "\n",
    "*Note: We removed Chatham island data here to ease plotting.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_hexmap(dset, ax=None):\n",
    "    if ax is None:\n",
    "        _, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "    hb = ax.hexbin(\n",
    "        dset[\"X\"], dset[\"Y\"], gridsize=500, cmap=\"BuPu\", mincnt=1, bins=\"log\"\n",
    "    )\n",
    "    ax.set_xlabel(\"Longitude\")\n",
    "    ax.set_ylabel(\"Latitude\")\n",
    "    cx.add_basemap(ax, crs=4326, source=cx.providers.CartoDB.Positron)\n",
    "\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "    ax.figure.colorbar(hb, cax=cax)\n",
    "\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_hexmap(dset[dset.X > 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In dense aread, like in Auckland, there are enough crashes events to map the\n",
    "local road network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dset_auckland = dset[dset[\"X\"].between(174.7, 174.9) & dset[\"Y\"].between(-37, -36.8)]\n",
    "plot_hexmap(dset_auckland)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At a coarser level, there is also the region information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_perc = dset[\"region\"].value_counts(normalize=True)\n",
    "ax = region_perc.plot.bar(ylabel=\"fraction of crashes\", figsize=(10, 5))\n",
    "_ = ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha=\"right\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f\"The top 4 regions account for {region_perc.nlargest(4).sum() * 100:0.1f}% \"\n",
    "    \"of the crashes.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Temporal features\n",
    "\n",
    "The dataset contains few temporal features:\n",
    "\n",
    "- `crashYear` and `crashFinancialYear`, respectively the year and final year\n",
    "  of each crash,\n",
    "- `holiday`, whether it occurs during a holiday period.\n",
    "\n",
    "So we won't be able to study daily, weekly and yearly patterns with these data.\n",
    "\n",
    "If we look at the yearly counts, we can see some fluctuations, mostly driven\n",
    "by Auckland region but still noticeable in other parts of the country.\n",
    "Year 2020 is much lower as it's the current year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "year_counts = dset[\"crashYear\"].value_counts(sort=False)\n",
    "_ = year_counts.plot.bar(ylabel=\"# crashes\", figsize=(10, 5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "year_region_counts = (\n",
    "    dset.groupby([\"crashYear\", \"region\"]).size().reset_index(name=\"# crashes\")\n",
    ")\n",
    "_, ax = plt.subplots(figsize=(10, 5))\n",
    "sb.pointplot(data=year_region_counts, x=\"crashYear\", y=\"# crashes\", hue=\"region\", ax=ax)\n",
    "ax.set_xticklabels(ax.get_xticklabels(), rotation=90)\n",
    "_ = plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also explore the spatio-temporal patterns too. Here we focus on\n",
    "Auckland (excluding 2020)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = sb.FacetGrid(\n",
    "    dset_auckland[dset_auckland.crashYear < 2020], col=\"crashYear\", col_wrap=5\n",
    ")\n",
    "grid.map(plt.hexbin, \"X\", \"Y\", gridsize=500, cmap=\"BuPu\", mincnt=1, bins=\"log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other temporal attribute is the holiday. Christmas is the holiday period\n",
    "with most of the accidents. How the period is computed is not clear, so the\n",
    "larger amount of accident could be partly due to the time extent. Easter,\n",
    "Queens Birthday and Labour weekend are 3 to 4 days periods. Christmas & New Year\n",
    "is probably 1 to 2 weeks period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "holiday_counts = dset[\"holiday\"].fillna(\"Normal day\").value_counts()\n",
    "ax = holiday_counts.plot.bar(ylabel=\"# crashes\", rot=0, figsize=(10, 5))\n",
    "_ = ax.set(yscale=\"log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Road Features\n",
    "\n",
    "From the dataset fields description, the following features seem specific to\n",
    "the type of road:\n",
    "\n",
    "- `crashSHDescription`, whether the crash happened on a state highway,\n",
    "- `flatHill`, whether the road is flat or sloped,\n",
    "- `junctionType`, type of junction the crash happened at (may also be *unknown*\n",
    "  & crashes not occurring at a junction are also *unknown*),\n",
    "- `NumberOfLanes`, number of lanes on the crash road,\n",
    "- `roadCharacter`, general nature of the road,\n",
    "- `roadCurvature`, simplified curvature of the road,\n",
    "- `roadLane`, lane configuration of the road (' ' for unknown or invalid\n",
    "  configurations),\n",
    "- `roadMarkings`, road markings at the crash site,\n",
    "- `roadSurface`, road surface description applying at the crash site,\n",
    "- `speedLimit`,  speed limit in force at the crash site at the time of the\n",
    "  crash (number, or 'LSZ' for a limited speed zone),\n",
    "- `streetLight`, street lighting at the time of the crash (this is also\n",
    "  a sort of **temporal information**),\n",
    "- `urban`, whether the road is in an urban area (derived from speed limit).\n",
    "\n",
    "Unfortunately, not all fields are actually available in the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "road_features = set(\n",
    "    [\n",
    "        \"crashSHDescription\",\n",
    "        \"flatHill\",\n",
    "        \"junctionType\",\n",
    "        \"NumberOfLanes\",\n",
    "        \"roadCharacter\",\n",
    "        \"roadCurvature\",\n",
    "        \"roadLane\",\n",
    "        \"roadMarkings\",\n",
    "        \"roadSurface\",\n",
    "        \"speedLimit\",\n",
    "        \"streetLight\",\n",
    "        \"urban\",\n",
    "    ]\n",
    ")\n",
    "missing_features = road_features - set(dset.columns)\n",
    "road_features -= missing_features\n",
    "print(\"The following features are not found in the dataset:\", missing_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, 3, figsize=(15, 12))\n",
    "for ax, feat in zip(axes.flat, sorted(road_features)):\n",
    "    counts = dset[feat].value_counts(dropna=False)\n",
    "    counts.plot.bar(ylabel=\"# crashes\", title=feat, ax=ax)\n",
    "    ax.set(yscale=\"log\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `urban` feature is derived from `speedLimit`, so we can probably remove it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Environmental features\n",
    "\n",
    "The environmental features are weather and sunhsine:\n",
    "\n",
    "- `light`, light at the time and place of the crash (this is also a sort of\n",
    "  **temporal information**),\n",
    "- `weatherA` and `weatherB`,  weather at the crash time/place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "env_features = [\"light\", \"weatherA\", \"weatherB\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(13, 4))\n",
    "for ax, feat in zip(axes.flat, env_features):\n",
    "    counts = dset[feat].value_counts(dropna=False)\n",
    "    counts.plot.bar(ylabel=\"# crashes\", title=feat, ax=ax)\n",
    "    ax.set(yscale=\"log\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Possible next steps\n",
    "\n",
    "We have checked the spatial, temporal, road and environmental features related\n",
    "to the accidents.\n",
    "\n",
    "If these features inform us in which conditions there are more accidents\n",
    "relatively, we will need additional baseline information if we want to create\n",
    "a predictive model.\n",
    "\n",
    "For the road features we could use a [LINZ dataset](https://data.linz.govt.nz/layer/50329-nz-road-centrelines-topo-150k/)\n",
    "or another [NZTA dataset](https://opendata-nzta.opendata.arcgis.com/datasets/NZTA::national-road-centreline-road-controlling-authority-data)\n",
    "that brings more information about the road type and traffic. But then we need\n",
    "to attribute each crash to a road.\n",
    "\n",
    "Another option would be to regrid the data, and for each cell containing at\n",
    "least one crash event we associate road features from the crash events. With\n",
    "this option, we don't make any prediction for cells of the grid where we don't\n",
    "have information about.\n",
    "\n",
    "For the environmental information, we need weather information for all days in\n",
    "a year:\n",
    "\n",
    "- for example from a climate reanalysis, e.g. [ERA 5](https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview),\n",
    "- or we can use the [station data from NIWA](https://cliflo.niwa.co.nz/),\n",
    "- or some simple [climate summaries](https://niwa.co.nz/education-and-training/schools/resources/climate/summary).\n",
    "\n",
    "The prediction task can be formulated in different ways:\n",
    "\n",
    "1. exclude weather & holiday features, and fit a regression model with count\n",
    "   data using year & location features (and accounting for traffic volume to\n",
    "   compare the number of crashes per car on the road),\n",
    "\n",
    "2. group data by location, time, weather type (e.g. rain vs. no rain), and\n",
    "   perform a binomial regression using the total number of days in each\n",
    "   category (e.g. number of rain days for a particular location & year),\n",
    "\n",
    "3. predict crash severity from the whole dataset, assuming the non-severe\n",
    "   crashes are a good proxy for normal conditions (weather, holidays, etc.)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Original computing environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!date -R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!uname -a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip freeze"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}