SI671-Project.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, math, random, re\n",
    "from datetime import datetime, timedelta, date\n",
    "\n",
    "import pandas as pd\n",
    "from pandas import DataFrame, Series\n",
    "import numpy as np\n",
    "\n",
    "from scipy.sparse import lil_matrix\n",
    "\n",
    "import pysal\n",
    "import pysal.esda as esda\n",
    "import pysal.esda.mapclassify as mc\n",
    "\n",
    "import geopandas as gpd\n",
    "from geopandas.tools import sjoin\n",
    "import geoplot as gplt\n",
    "from shapely.geometry import Point, Polygon\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colors\n",
    "import seaborn as sns\n",
    "\n",
    "from util import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 23.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "detail_data=DataFrame()\n",
    "loaded_files = []\n",
    "for filename in os.listdir(\"data/details\"):\n",
    "    if filename.endswith(\".csv\"):\n",
    "        loaded_files.append(pd.read_csv(\"data/details/\"+filename))\n",
    "detail_data = pd.concat(loaded_files)\n",
    "detail_data.reset_index(drop=True,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>BEGIN_YEARMONTH</th>\n",
       "      <th>BEGIN_DAY</th>\n",
       "      <th>BEGIN_TIME</th>\n",
       "      <th>END_YEARMONTH</th>\n",
       "      <th>END_DAY</th>\n",
       "      <th>END_TIME</th>\n",
       "      <th>EPISODE_ID</th>\n",
       "      <th>EVENT_ID</th>\n",
       "      <th>STATE</th>\n",
       "      <th>STATE_FIPS</th>\n",
       "      <th>...</th>\n",
       "      <th>END_RANGE</th>\n",
       "      <th>END_AZIMUTH</th>\n",
       "      <th>END_LOCATION</th>\n",
       "      <th>BEGIN_LAT</th>\n",
       "      <th>BEGIN_LON</th>\n",
       "      <th>END_LAT</th>\n",
       "      <th>END_LON</th>\n",
       "      <th>EPISODE_NARRATIVE</th>\n",
       "      <th>EVENT_NARRATIVE</th>\n",
       "      <th>DATA_SOURCE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>195004</td>\n",
       "      <td>28</td>\n",
       "      <td>1445</td>\n",
       "      <td>195004</td>\n",
       "      <td>28</td>\n",
       "      <td>1445</td>\n",
       "      <td>NaN</td>\n",
       "      <td>10096222</td>\n",
       "      <td>OKLAHOMA</td>\n",
       "      <td>40.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>35.12</td>\n",
       "      <td>-99.20</td>\n",
       "      <td>35.17</td>\n",
       "      <td>-99.20</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>PUB</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>195004</td>\n",
       "      <td>29</td>\n",
       "      <td>1530</td>\n",
       "      <td>195004</td>\n",
       "      <td>29</td>\n",
       "      <td>1530</td>\n",
       "      <td>NaN</td>\n",
       "      <td>10120412</td>\n",
       "      <td>TEXAS</td>\n",
       "      <td>48.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>31.90</td>\n",
       "      <td>-98.60</td>\n",
       "      <td>31.73</td>\n",
       "      <td>-98.60</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>PUB</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>195007</td>\n",
       "      <td>5</td>\n",
       "      <td>1800</td>\n",
       "      <td>195007</td>\n",
       "      <td>5</td>\n",
       "      <td>1800</td>\n",
       "      <td>NaN</td>\n",
       "      <td>10104927</td>\n",
       "      <td>PENNSYLVANIA</td>\n",
       "      <td>42.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>40.58</td>\n",
       "      <td>-75.70</td>\n",
       "      <td>40.65</td>\n",
       "      <td>-75.47</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>PUB</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>195007</td>\n",
       "      <td>5</td>\n",
       "      <td>1830</td>\n",
       "      <td>195007</td>\n",
       "      <td>5</td>\n",
       "      <td>1830</td>\n",
       "      <td>NaN</td>\n",
       "      <td>10104928</td>\n",
       "      <td>PENNSYLVANIA</td>\n",
       "      <td>42.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>40.60</td>\n",
       "      <td>-76.75</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>PUB</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>195007</td>\n",
       "      <td>24</td>\n",
       "      <td>1440</td>\n",
       "      <td>195007</td>\n",
       "      <td>24</td>\n",
       "      <td>1440</td>\n",
       "      <td>NaN</td>\n",
       "      <td>10104929</td>\n",
       "      <td>PENNSYLVANIA</td>\n",
       "      <td>42.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>41.63</td>\n",
       "      <td>-79.68</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>PUB</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 51 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   BEGIN_YEARMONTH  BEGIN_DAY  BEGIN_TIME  END_YEARMONTH  END_DAY  END_TIME  \\\n",
       "0           195004         28        1445         195004       28      1445   \n",
       "1           195004         29        1530         195004       29      1530   \n",
       "2           195007          5        1800         195007        5      1800   \n",
       "3           195007          5        1830         195007        5      1830   \n",
       "4           195007         24        1440         195007       24      1440   \n",
       "\n",
       "   EPISODE_ID  EVENT_ID         STATE  STATE_FIPS     ...      END_RANGE  \\\n",
       "0         NaN  10096222      OKLAHOMA        40.0     ...            0.0   \n",
       "1         NaN  10120412         TEXAS        48.0     ...            0.0   \n",
       "2         NaN  10104927  PENNSYLVANIA        42.0     ...            0.0   \n",
       "3         NaN  10104928  PENNSYLVANIA        42.0     ...            0.0   \n",
       "4         NaN  10104929  PENNSYLVANIA        42.0     ...            0.0   \n",
       "\n",
       "  END_AZIMUTH END_LOCATION BEGIN_LAT  BEGIN_LON END_LAT END_LON  \\\n",
       "0         NaN          NaN     35.12     -99.20   35.17  -99.20   \n",
       "1         NaN          NaN     31.90     -98.60   31.73  -98.60   \n",
       "2         NaN          NaN     40.58     -75.70   40.65  -75.47   \n",
       "3         NaN          NaN     40.60     -76.75     NaN     NaN   \n",
       "4         NaN          NaN     41.63     -79.68     NaN     NaN   \n",
       "\n",
       "  EPISODE_NARRATIVE EVENT_NARRATIVE DATA_SOURCE  \n",
       "0               NaN             NaN         PUB  \n",
       "1               NaN             NaN         PUB  \n",
       "2               NaN             NaN         PUB  \n",
       "3               NaN             NaN         PUB  \n",
       "4               NaN             NaN         PUB  \n",
       "\n",
       "[5 rows x 51 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "detail_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert the Time related strings to a DateTime column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 28.8 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Pad the \"DAY\" columns to length 2 so the date parser can do it's thing\n",
    "detail_data.BEGIN_DAY = detail_data.BEGIN_DAY.map(str).str.pad(2,fillchar='0')\n",
    "detail_data.END_DAY = detail_data.END_DAY.map(str).str.pad(2,fillchar='0')\n",
    "# ditto with the Time\n",
    "detail_data.BEGIN_TIME = detail_data.BEGIN_TIME.map(str).str.pad(4,fillchar='0')\n",
    "detail_data.END_TIME = detail_data.END_TIME.map(str).str.pad(4,fillchar='0')\n",
    "# create a new column by concating the three date/time related columns and convert the result to a datetime\n",
    "detail_data['BEGIN_DATE']=detail_data.BEGIN_YEARMONTH.map(str)+\" \"+detail_data.BEGIN_DAY.map(str)+\" \"+detail_data.BEGIN_TIME.map(str)\n",
    "detail_data.BEGIN_DATE=pd.to_datetime(detail_data.BEGIN_DATE, format='%Y%m %d %H%M', errors='coerce')\n",
    "\n",
    "detail_data['END_DATE']=detail_data.END_YEARMONTH.map(str)+\" \"+detail_data.END_DAY.map(str)+\" \"+detail_data.END_TIME.map(str)\n",
    "detail_data.END_DATE=pd.to_datetime(detail_data.END_DATE, format='%Y%m %d %H%M', errors='coerce')\n",
    "\n",
    "# drop the old columns\n",
    "detail_data = detail_data.drop(['BEGIN_YEARMONTH', \"END_YEARMONTH\",'BEGIN_DAY', \"END_DAY\", \"BEGIN_TIME\", \"END_TIME\"], axis=1)\n",
    "detail_data = detail_data.drop(['YEAR', \"MONTH_NAME\",'END_DATE_TIME', \"BEGIN_DATE_TIME\"], axis=1)\n",
    "\n",
    "# Columns I don't think I need at the momment\n",
    "detail_data = detail_data.drop(['WFO', 'SOURCE',\n",
    "       'MAGNITUDE', 'MAGNITUDE_TYPE', 'FLOOD_CAUSE', 'CATEGORY', 'TOR_F_SCALE',\n",
    "       'TOR_LENGTH', 'TOR_WIDTH', 'TOR_OTHER_WFO', 'TOR_OTHER_CZ_STATE',\n",
    "       'TOR_OTHER_CZ_FIPS', 'TOR_OTHER_CZ_NAME', 'BEGIN_RANGE',\n",
    "       'BEGIN_AZIMUTH', 'BEGIN_LOCATION', 'END_RANGE', 'END_AZIMUTH',\n",
    "       'END_LOCATION','EPISODE_NARRATIVE', 'EVENT_NARRATIVE', 'DATA_SOURCE'], axis=1)\n",
    "\n",
    "# columns I'm less sure I don't need, but it'll make things easier for the time being\n",
    "detail_data = detail_data.drop(['INJURIES_DIRECT',\n",
    "       'INJURIES_INDIRECT', 'DEATHS_DIRECT', 'DEATHS_INDIRECT',\n",
    "       'DAMAGE_PROPERTY', 'DAMAGE_CROPS'], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Correcting input errors\n",
    "There's some input errors with the Lat/Long coordinates - decimal place is just shifted, so dividing those by ten gives us the correct value (Affects 2080 rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['EPISODE_ID', 'EVENT_ID', 'STATE', 'STATE_FIPS', 'EVENT_TYPE',\n",
       "       'CZ_TYPE', 'CZ_FIPS', 'CZ_NAME', 'CZ_TIMEZONE', 'BEGIN_LAT',\n",
       "       'BEGIN_LON', 'END_LAT', 'END_LON', 'BEGIN_DATE', 'END_DATE'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "detail_data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "detail_data.loc[detail_data.BEGIN_LON < -180, 'BEGIN_LON']=detail_data.loc[detail_data.BEGIN_LON < -180, 'BEGIN_LON']/10\n",
    "detail_data.loc[detail_data.END_LON < -180, 'END_LON']=detail_data.loc[detail_data.END_LON < -180, 'END_LON']/10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fix inconsistencies in event type tagging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "detail_data.loc[detail_data.EVENT_TYPE.isin(['HAIL FLOODING', 'HAIL/ICY ROADS']), 'EVENT_TYPE'] = 'Hail'\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(['Heavy Wind']), 'EVENT_TYPE'] = 'High Wind'\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(['High Snow']), 'EVENT_TYPE'] = 'Heavy Snow'\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(['High Snow']), 'EVENT_TYPE'] = 'Heavy Snow'\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(\n",
    "    [\n",
    "        'TORNADOES, TSTM WIND, HAIL',\n",
    "        'THUNDERSTORM WINDS/FLOODING',\n",
    "        'THUNDERSTORM WINDS/FLASH FLOOD',\n",
    "        'THUNDERSTORM WINDS LIGHTNING',\n",
    "        'THUNDERSTORM WIND/ TREES',\n",
    "        'THUNDERSTORM WIND/ TREE',\n",
    "        'THUNDERSTORM WINDS FUNNEL CLOU',\n",
    "        'THUNDERSTORM WINDS/HEAVY RAIN',\n",
    "        'THUNDERSTORM WINDS HEAVY RAIN',\n",
    "        'THUNDERSTORM WINDS/ FLOOD'\n",
    "    ]\n",
    "), 'EVENT_TYPE'] = 'Thunderstorm Wind'\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(['TORNADO/WATERSPOUT']), 'EVENT_TYPE'] = 'Tornado'\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(['Landslide']), 'EVENT_TYPE'] = 'Debris Flow'\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(['Volcanic Ashfall']), 'EVENT_TYPE'] = 'Volcanic Ash'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop 14 rows tagged \"Northern Lights\" and a single row tagged 'OTHER'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "detail_data.drop(detail_data[detail_data.EVENT_TYPE.isin(['Northern Lights', 'OTHER'])].index, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trim the data to just the Continental US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "detail_data = detail_data[\n",
    "        ((detail_data.BEGIN_LON > -124.7844079) &\n",
    "        (detail_data.BEGIN_LON < -66.9513812) &\n",
    "        (detail_data.BEGIN_LAT > 24.7433195) &\n",
    "        (detail_data.BEGIN_LAT < 49.3457868)) |\n",
    "        (detail_data.BEGIN_LAT.isnull())\n",
    "    ].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create Categorical Index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create two categories - one based off of the event type and another based on my not very scientific classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "detail_data.EVENT_TYPE = detail_data.EVENT_TYPE.astype(\"category\")\n",
    "detail_data['EVENT_CODE'] = detail_data.EVENT_TYPE.cat.codes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "detail_data['META_TYPE'] = None\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(\n",
    "    [\n",
    "        'Avalanche',\n",
    "        'Blizzard',\n",
    "        'Cold/Wind Chill',\n",
    "        'Extreme Cold/Wind Chill',\n",
    "        'Freezing Fog',\n",
    "        'Frost/Freeze',\n",
    "        'Hail',\n",
    "        'Heavy Snow',\n",
    "        'Ice Storm',\n",
    "        'Lake-Effect Snow',\n",
    "        'Marine Hail',\n",
    "        'Sleet',\n",
    "        'Winter Storm',\n",
    "        'Winter Weather'\n",
    "    ]\n",
    "), 'META_TYPE'] = 'Cold'\n",
    "\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(\n",
    "    [\n",
    "        'Dense Smoke',\n",
    "        'Drought',\n",
    "        'Excessive Heat',\n",
    "        'Heat',\n",
    "        'Volcanic Ash',\n",
    "        'Wildfire'\n",
    "    ]\n",
    "), 'META_TYPE'] = 'Heat'\n",
    "\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(\n",
    "    [\n",
    "        'Debris Flow',\n",
    "        'Dense Fog',\n",
    "        'Lightning',\n",
    "        'Marine Dense Fog',\n",
    "        'Marine Lightning'\n",
    "    ]\n",
    "), 'META_TYPE'] = 'Other'\n",
    "\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(\n",
    "    [\n",
    "        'Astronomical Low Tide',\n",
    "        'Coastal Flood',\n",
    "        'Flash Flood',\n",
    "        'Flood',\n",
    "        'Heavy Rain',\n",
    "        'High Surf',\n",
    "        'Lakeshore Flood',\n",
    "        'Rip Current',\n",
    "        'Seiche',\n",
    "        'Sneakerwave',\n",
    "        'Storm Surge/Tide',\n",
    "        'Tsunami',\n",
    "        'Waterspout'\n",
    "    ]\n",
    "), 'META_TYPE'] = 'Water'\n",
    "\n",
    "detail_data.loc[detail_data.EVENT_TYPE.isin(\n",
    "    [\n",
    "        'Dust Devil',\n",
    "        'Dust Storm',\n",
    "        'Funnel Cloud',\n",
    "        'High Wind',\n",
    "        'Hurricane',\n",
    "        'Hurricane (Typhoon)',\n",
    "        'Marine High Wind',\n",
    "        'Marine Hurricane/Typhoon',\n",
    "        'Marine Strong Wind',\n",
    "        'Marine Thunderstorm Wind',\n",
    "        'Marine Tropical Depression',\n",
    "        'Marine Tropical Storm',\n",
    "        'Strong Wind',\n",
    "        'Thunderstorm Wind',\n",
    "        'Tornado',\n",
    "        'Tropical Depression',\n",
    "        'Tropical Storm'\n",
    "    ]\n",
    "), 'META_TYPE'] = 'Wind'\n",
    "detail_data.META_TYPE = detail_data.META_TYPE.astype(\"category\")\n",
    "detail_data['META_CODE'] = detail_data.META_TYPE.cat.codes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fill in missing data\n",
    "The US Census Gazetteer contains a file which provides Lat/Long for every county in the US. Combined with the FIPs codes, we can give estimated locations for events with no clear starting point (about 1/3 of the dataset - things like \"drought\" or \"Strong Wind\" that don't have a clearly defined geographic location).\n",
    "Dataset from https://www.census.gov/geo/maps-data/data/gazetteer2017.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>USPS</th>\n",
       "      <th>GEOID</th>\n",
       "      <th>ANSICODE</th>\n",
       "      <th>NAME</th>\n",
       "      <th>ALAND</th>\n",
       "      <th>AWATER</th>\n",
       "      <th>ALAND_SQMI</th>\n",
       "      <th>AWATER_SQMI</th>\n",
       "      <th>INTPTLAT</th>\n",
       "      <th>INTPTLONG</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FIPS</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>01001</th>\n",
       "      <td>AL</td>\n",
       "      <td>1001</td>\n",
       "      <td>161526</td>\n",
       "      <td>Autauga County</td>\n",
       "      <td>1539614693</td>\n",
       "      <td>25744269</td>\n",
       "      <td>594.449</td>\n",
       "      <td>9.940</td>\n",
       "      <td>32.532237</td>\n",
       "      <td>-86.646440</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>01003</th>\n",
       "      <td>AL</td>\n",
       "      <td>1003</td>\n",
       "      <td>161527</td>\n",
       "      <td>Baldwin County</td>\n",
       "      <td>4117605847</td>\n",
       "      <td>1133109409</td>\n",
       "      <td>1589.817</td>\n",
       "      <td>437.496</td>\n",
       "      <td>30.659218</td>\n",
       "      <td>-87.746067</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>01005</th>\n",
       "      <td>AL</td>\n",
       "      <td>1005</td>\n",
       "      <td>161528</td>\n",
       "      <td>Barbour County</td>\n",
       "      <td>2292144656</td>\n",
       "      <td>50538698</td>\n",
       "      <td>885.002</td>\n",
       "      <td>19.513</td>\n",
       "      <td>31.870253</td>\n",
       "      <td>-85.405104</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      USPS  GEOID  ANSICODE            NAME       ALAND      AWATER  \\\n",
       "FIPS                                                                  \n",
       "01001   AL   1001    161526  Autauga County  1539614693    25744269   \n",
       "01003   AL   1003    161527  Baldwin County  4117605847  1133109409   \n",
       "01005   AL   1005    161528  Barbour County  2292144656    50538698   \n",
       "\n",
       "       ALAND_SQMI  AWATER_SQMI   INTPTLAT  INTPTLONG  \n",
       "FIPS                                                  \n",
       "01001     594.449        9.940  32.532237 -86.646440  \n",
       "01003    1589.817      437.496  30.659218 -87.746067  \n",
       "01005     885.002       19.513  31.870253 -85.405104  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "county_df = pd.read_csv('data/2017_Gaz_counties_national.txt', sep='\\t', engine='python')\n",
    "county_df['FIPS'] = county_df['GEOID'].apply(lambda x: '{0:0>5}'.format(x))\n",
    "county_df.set_index(\"FIPS\", inplace=True)\n",
    "county_df.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fill in Lat/Lon for events based on county centroid\n",
    "Create a new column that combines the State and County FIPS so we can look them up (gets coordinates for ~250k entries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "detail_data['FIPS'] = detail_data['STATE_FIPS']*1000+detail_data['CZ_FIPS']\n",
    "detail_data['FIPS'].fillna(0, inplace=True)\n",
    "detail_data['FIPS'] = detail_data['FIPS'].apply(lambda x: '{0:0>5}'.format(int(x)))\n",
    "\n",
    "# drop the old columns, since we've condensed that data into one column\n",
    "detail_data = detail_data.drop(['EPISODE_ID', 'EVENT_ID', 'STATE_FIPS', 'CZ_TYPE', 'CZ_FIPS', 'CZ_TIMEZONE'], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try to match up FIPS ids for events with no location info, and fill in the begin lat/lon with that counties centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# temporaily assing FIPS as the index\n",
    "detail_data.set_index(\"FIPS\", inplace=True)\n",
    "\n",
    "# left join on the county dataframe to get a column of county centroids\n",
    "detail_data = detail_data.join(county_df[['INTPTLAT', 'INTPTLONG']], on=\"FIPS\")\n",
    "\n",
    "# reset the index to avoid duplicate issues in a sec \n",
    "detail_data.reset_index(inplace=True)\n",
    "\n",
    "# if the lat/lon are blank, and we were able to match it to a FIPS code, add the country centroid as the event lat/lon\n",
    "detail_data.loc[detail_data.BEGIN_LAT.isnull(), 'BEGIN_LAT'] = detail_data.INTPTLAT\n",
    "detail_data.loc[detail_data.BEGIN_LON.isnull(), 'BEGIN_LON'] = detail_data.INTPTLONG\n",
    "\n",
    "#drop the extra columns now that we don't need them\n",
    "detail_data.drop(['INTPTLAT', 'INTPTLONG'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fill in missing Date Entries\n",
    "Some entries have an end date but no starting date. We're just going to fill in the latter with the former (~110k entries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "detail_data.loc[detail_data.BEGIN_DATE.isnull(), 'BEGIN_DATE'] = detail_data.END_DATE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# redundant after all the date messing earlier, but these turn out to be real handy\n",
    "detail_data['MONTH'] = detail_data.BEGIN_DATE.dt.month\n",
    "detail_data['DAY'] = detail_data.BEGIN_DATE.dt.day"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO: look at filling in the remaining missing points with bp02oc18.dbx from https://www.weather.gov/gis/ZoneCounty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drop rows without enough data to do anything with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "detail_data = detail_data.drop(detail_data[\n",
    "    (detail_data.BEGIN_LON.isnull()) |\n",
    "    (detail_data.BEGIN_LAT.isnull()) |\n",
    "    (detail_data.BEGIN_DATE.isnull())\n",
    "].index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create the matricies and Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a matrix based on a given data and a +/- date span on eithe side\n",
    "def get_day_matrix(date,span=30,meta_type=None,event_type=None,coords=None,coord_range=None):\n",
    "    data_set = []\n",
    "    for year in range(1950,2017):\n",
    "        # set the year in the date\n",
    "        date = date.replace(year=year)\n",
    "\n",
    "        # set the start and end based on the \n",
    "        start = date - timedelta(days=span)\n",
    "        end = date + timedelta(days=span)\n",
    "        \n",
    "        \n",
    "    \n",
    "        year_data = detail_data[\n",
    "            (detail_data.BEGIN_DATE>pd.to_datetime(start)) &\n",
    "            (detail_data.BEGIN_DATE<pd.to_datetime(end))\n",
    "        ]\n",
    "        if meta_type is not None:\n",
    "            year_data = year_data[(detail_data.META_TYPE == meta_type)]\n",
    "        if event_type is not None:\n",
    "            year_data = year_data[(detail_data.EVENT_TYPE == event_type)]\n",
    "            \n",
    "        if coords is not None:\n",
    "            year_data = year_data[\n",
    "                (year_data.BEGIN_LAT > coords[0]-coord_range) &\n",
    "                (year_data.BEGIN_LAT < coords[0]+coord_range) &\n",
    "                (year_data.BEGIN_LON > coords[1]-coord_range) &\n",
    "                (year_data.BEGIN_LON < coords[1]+coord_range)\n",
    "            ]\n",
    "                \n",
    "        data_set.append(year_data)\n",
    "\n",
    "    data = pd.concat(data_set)\n",
    "\n",
    "    m = lil_matrix((5783,2760), dtype=np.float16)\n",
    "\n",
    "    for i, row in data.iterrows():\n",
    "        # calculate the discount applied to the matrix value based on how many days off we are\n",
    "        date = date.replace(year=row.BEGIN_DATE.year)\n",
    "        days_off = abs((row.BEGIN_DATE-date).days)\n",
    "        if days_off == 0:\n",
    "            mat_val = 1\n",
    "        else:\n",
    "            mat_val = 1/(1+math.log(days_off))\n",
    "\n",
    "        # is this a single point event or a line\n",
    "        if pd.isnull(row.END_LAT) or pd.isnull(row.END_LON):\n",
    "            # single point - plot a square\n",
    "            col_id = lat_to_index(row.BEGIN_LAT)\n",
    "            row_id = lon_to_index(row.BEGIN_LON)\n",
    "            for r in range(row_id-3, row_id+4):\n",
    "                for c in range(col_id-3, col_id+4):\n",
    "                    try:\n",
    "                        m[r,c] += mat_val\n",
    "                    except:\n",
    "                        print(\"bad row:\",row)\n",
    "        else:\n",
    "            # we have a start and an end, so it's a line\n",
    "            start = (lat_to_index(row.BEGIN_LAT),lon_to_index(row.BEGIN_LON))\n",
    "            end = (lat_to_index(row.END_LAT),lon_to_index(row.END_LON))\n",
    "            for pos in bresenham_line(start,end):\n",
    "                try:\n",
    "                    m[pos[1],pos[0]] += mat_val\n",
    "                except:\n",
    "                    print(\"bad row (line):\",pos)\n",
    "    return m\n",
    "\n",
    "# Examples:\n",
    "#\n",
    "# get_day_matrix(\n",
    "#     datetime.strptime('11/15/18', '%m/%d/%y'),\n",
    "#     event_type=\"Dense Fog\"\n",
    "# )\n",
    "\n",
    "# get_day_matrix(\n",
    "#     date, \n",
    "#     span=30, \n",
    "#     meta_type=\"Cold\", \n",
    "#     event_type=None,\n",
    "#     coords=coords,\n",
    "#     coord_range=2.5\n",
    "# )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creates a Geopandas Dataframe from the matrix output of get_day_matrix\n",
    "def create_geopandas_from_matrix(m):\n",
    "    # get the row and col arrays from the non zero elements of the matrix\n",
    "    row,col = m.nonzero()\n",
    "    # iterate through the arrays and convert them to their lat/lon approximations\n",
    "    data=[]\n",
    "    for x in range(len(row)):\n",
    "        data.append({\n",
    "            \"lat\":index_to_lat(col[x]),\n",
    "            \"lon\":index_to_lon(row[x]),\n",
    "            \"x\":m[row[x],col[x]]\n",
    "        })\n",
    "    # create a pandas dataframe from the data, then convert it to a geopandas dataframe\n",
    "    data = DataFrame(data)\n",
    "    if len(data) > 0:\n",
    "        data['Coordinates'] = list(zip(data.lon, data.lat))\n",
    "        data['Coordinates'] = data['Coordinates'].apply(Point)\n",
    "        return gpd.GeoDataFrame(data, geometry='Coordinates')\n",
    "    else:\n",
    "        return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "load in the us state map (from https://www.arcgis.com/home/item.html?id=b07a9393ecbd430795a6f6218443dccc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_state_shapemap():\n",
    "    us_map = gpd.read_file(\"data/states_shape/states.shp\")\n",
    "    return us_map[~us_map.STATE_NAME.isin(['Hawaii',\"Alaska\"])].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "load in the us county map (from https://www.census.gov/geo/maps-data/data/cbf/cbf_counties.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_county_shapemap():\n",
    "    us_county_map = gpd.read_file(\"data/shape/us_counties/cb_2017_us_county_500k.shp\")\n",
    "    us_county_map.STATEFP = us_county_map.STATEFP.astype(int)\n",
    "    # constrain the county map to the continental us\n",
    "    return us_county_map[\n",
    "        (us_county_map.STATEFP<60) &\n",
    "        (us_county_map.STATEFP!=2) &\n",
    "        (us_county_map.STATEFP!=15)\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the distance weights\n",
    "wq = pysal.weights.DistanceBand.from_dataframe(get_county_shapemap(), threshold=3)\n",
    "# weigh each county based on the total number of bordering counties\n",
    "wq.transform = 'r'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#DEMO PLOT - County plot by straight count\n",
    "date = datetime.strptime('04/15/18', '%m/%d/%y')\n",
    "m = get_day_matrix(date, span=10, event_type=\"Tornado\")\n",
    "data = create_geopandas_from_matrix(m)\n",
    "\n",
    "us_county_map = get_county_shapemap()\n",
    "\n",
    "data.crs = us_county_map.crs\n",
    "county_data = sjoin(data, us_county_map, how='left')\n",
    "us_county_map.set_index(\"GEOID\", inplace=True)\n",
    "us_county_map['counts'] = county_data.groupby(\"GEOID\")['x'].sum()\n",
    "us_county_map.counts.fillna(0, inplace=True)\n",
    "\n",
    "ax = us_county_map.plot(column=\"counts\", figsize=(5,5))\n",
    "\n",
    "# We can now plot our GeoDataFrame.\n",
    "plt.title(\"Points Condensed into Counties\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a gettis ord g plot provided a dataframe from 'create_geopandas_from_matrix(m)'\n",
    "def make_gettis_ord_g_plot(data, title=\"Gettis Ord G* Autocorrelation\", filename=None, color_map=\"GnBu\", ret_values=False):\n",
    "    # get the county shapemap\n",
    "    us_county_map = get_county_shapemap()\n",
    "\n",
    "    # join all of the data into the county dataframe\n",
    "    data.crs = us_county_map.crs\n",
    "    county_data = sjoin(data, us_county_map, how='left')\n",
    "    us_county_map.set_index(\"GEOID\", inplace=True)\n",
    "    us_county_map['counts'] = county_data.groupby(\"GEOID\")['x'].sum()\n",
    "    us_county_map.counts.fillna(0, inplace=True)\n",
    "    \n",
    "    # calculate the gettis ord gi of each county based on the data\n",
    "    us_county_map['counts_g_local'] = esda.getisord.G_Local(us_county_map['counts'], wq, permutations=5).Zs\n",
    "\n",
    "    # plot the result\n",
    "    f, ax = plt.subplots(1, figsize=(20,20))\n",
    "    us_county_map.plot(column='counts_g_local', cmap=color_map, linewidth=0.1, ax=ax, edgecolor='black')\n",
    "    ax.set_axis_off()\n",
    "    \n",
    "    # If there's a filename, just save it, if not, plot it to the notebook\n",
    "    plt.title(title)\n",
    "    if filename is not None:\n",
    "        plt.savefig(filename)\n",
    "        plt.close()\n",
    "    else:\n",
    "        plt.show()\n",
    "        \n",
    "    # in case I want to get the results back\n",
    "    if ret_values:\n",
    "        return us_county_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "make_gettis_ord_g_plot(data, title=\"Gettis Ord G* Autocorrelation - Tornado - 04/15\", color_map=\"Reds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get just the getti ord gi score for a given dataframe from 'create_geopandas_from_matrix(m)'\n",
    "def get_gettis_ord_g_scores(data):\n",
    "    # get the county shapemap\n",
    "    us_county_map = get_county_shapemap()\n",
    "\n",
    "    # join all of the data into the county dataframe\n",
    "    data.crs = us_county_map.crs\n",
    "    county_data = sjoin(data, us_county_map, how='left')\n",
    "    us_county_map.set_index(\"GEOID\", inplace=True)\n",
    "    us_county_map['counts'] = county_data.groupby(\"GEOID\")['x'].sum()\n",
    "    us_county_map.counts.fillna(0, inplace=True)\n",
    "    \n",
    "    return esda.getisord.G_Local(us_county_map['counts'], wq, permutations=5).Zs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "get_gettis_ord_g_scores(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create Gettis Ord Gi Dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_spatial_lag_quartiles(data, title=\"Spatial Lag Median Risk Quartiles\", filename=None, color_map=\"GnBu\", ret_values=False):\n",
    "    us_county_map = get_county_shapemap()\n",
    "\n",
    "    data.crs = us_county_map.crs\n",
    "    county_data = sjoin(data, us_county_map, how='left')\n",
    "    us_county_map.set_index(\"GEOID\", inplace=True)\n",
    "    us_county_map['counts'] = county_data.groupby(\"GEOID\")['x'].sum()\n",
    "    us_county_map.counts.fillna(0, inplace=True)\n",
    "\n",
    "    # create a list of counties and all the counties that county borders from the geometry of the shape file\n",
    "    # wq = pysal.weights.Queen.from_dataframe(us_county_map)\n",
    "    wq = pysal.weights.DistanceBand.from_dataframe(us_county_map, threshold=2, binary=True)\n",
    "    # weigh each county based on the total number of bordering counties\n",
    "    wq.transform = 'r'\n",
    "\n",
    "    # get the spatial lag for each county - based on all of the attrubute values we see in neighboring counties\n",
    "    y = us_county_map['counts']\n",
    "    ylag = pysal.weights.lag_spatial(wq,y)\n",
    "\n",
    "    f, ax = plt.subplots(1, figsize=(20,20))\n",
    "    us_county_map.assign(cl=ylag).plot(column='cl', categorical=True,\n",
    "            k=5, cmap=color_map, linewidth=0.1, ax=ax,\n",
    "            edgecolor='black')\n",
    "    ax.set_axis_off()\n",
    "    \n",
    "    plt.title(title)\n",
    "    if filename is not None:\n",
    "        plt.savefig(filename)\n",
    "        plt.close()\n",
    "    else:\n",
    "        plt.show()\n",
    "    \n",
    "    if ret_values:\n",
    "        return us_county_map.assign(cl=ylag)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# restricting predictions to these types at the moment\n",
    "event_types = ['Blizzard', 'Extreme Cold/Wind Chill', 'Hail', 'Heavy Snow', 'Ice Storm', 'Winter Storm','Drought', \n",
    "    'Excessive Heat', 'Wildfire', 'Flash Flood', 'Heavy Rain', 'High Wind', 'Hurricane', 'Thunderstorm Wind', 'Tornado']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# running overnight - tack on the others to run after the main ones\n",
    "# event_types = event_types + [e for e in detail_data.EVENT_TYPE.unique() if e not in event_types]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "for event_type in event_types:\n",
    "    date = datetime.strptime('01/01/18', '%m/%d/%y')\n",
    "    # set up the county data \n",
    "    county_data = get_county_shapemap()\n",
    "    county_data.set_index(\"GEOID\", inplace=True)\n",
    "    \n",
    "    while date.year<2019:\n",
    "        print(\"Plotting {} - {}\".format(event_type, str(date).split()[0]))\n",
    "\n",
    "        # make a matrix for this day, using events +/- 2 weeks\n",
    "        m = get_day_matrix(date, span=14, event_type=event_type)\n",
    "        # throw that data into a geopandas frame\n",
    "        data = create_geopandas_from_matrix(m)\n",
    "        if data is None:\n",
    "            print(\"No data found for\", date)\n",
    "            county_data[\"{}-{:%m-%d}_g_score\".format(re.sub(r'\\W+', '', event_type), date)] = 0\n",
    "        else:\n",
    "            vals = get_gettis_ord_g_scores(data)\n",
    "            county_data[\"{}-{:%m-%d}_g_score\".format(re.sub(r'\\W+', '', event_type), date)] = vals\n",
    "\n",
    "        date = date + timedelta(days=7)        \n",
    "    # save it\n",
    "    county_data.to_file(\"output/{}_weekly_g_score.geoJSON\".format(re.sub(r'\\W+', '', event_type)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make Predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# restricting predictions to these types at the moment\n",
    "event_types = ['Blizzard', 'Extreme Cold/Wind Chill', 'Hail', 'Heavy Snow', 'Ice Storm', 'Winter Storm','Drought', \n",
    "    'Excessive Heat', 'Wildfire', 'Flash Flood', 'Heavy Rain', 'High Wind', 'Hurricane', 'Thunderstorm Wind', 'Tornado']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting labeled events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# super simple - find an event at random that occured after 2016\n",
    "def get_positive_event():  \n",
    "    event = detail_data[detail_data.BEGIN_DATE > '2016-01-01'].sample()\n",
    "    return event.EVENT_TYPE.values[0], event.META_TYPE.values[0], (event.BEGIN_LAT.values[0], event.BEGIN_LON.values[0]), event.BEGIN_DATE.values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# trickier - we need to pick a location and date at random, then check to make sure nothing happened \n",
    "# if it did, we guess again until we come up with nothing\n",
    "def get_negative_event(max_range=1.5, day_range=10):\n",
    "    event_count = 99\n",
    "    event_type = random.choice(event_types)\n",
    "    while event_count != 0:\n",
    "        event_count = 0\n",
    "        # guess a random lat/long combo in the continental us\n",
    "        lat = random.uniform(24.7433195, 49.3457868)\n",
    "        lon = random.uniform(-124.7844079, -66.9513812)\n",
    "        \n",
    "        # guess a random date (days from 1-28 so I don't have to worry about Feb)\n",
    "        search_date = datetime(year=2016, month=random.randint(1,12), day=random.randint(1,28))\n",
    "        \n",
    "        # get all events from the dataset withing +/- max_range of that location\n",
    "        nearby_events = detail_data[\n",
    "            ((detail_data.BEGIN_LON > lon-max_range) &\n",
    "            (detail_data.BEGIN_LON < lon+max_range) &\n",
    "            (detail_data.BEGIN_LAT > lat-max_range) &\n",
    "            (detail_data.BEGIN_LAT < lat+max_range)) &\n",
    "            (detail_data.EVENT_TYPE==event_type)\n",
    "        ].copy()\n",
    "        \n",
    "        # check every year in the dataset and count how many hits we get\n",
    "        # TODO: There's got to be a simpler way of doing this\n",
    "        for year in range(1950, 2017):\n",
    "            search_date = search_date.replace(year=year)\n",
    "            event_count += len(nearby_events[\n",
    "                ((nearby_events.BEGIN_DATE >= (search_date - timedelta(days=day_range))) &\n",
    "                (nearby_events.BEGIN_DATE <= (search_date + timedelta(days=day_range))))\n",
    "            ])\n",
    "#         print(\"{} {} events at ({},{}) within {} days of {}\".format(event_count, event_type, lat, lon, day_range, search_date))\n",
    "    return event_type, \"\", (lat, lon), search_date"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prediction Approach #1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "rolling = 7\n",
    "threshold = 0.01\n",
    "\n",
    "def predict(event_type, coords, date):\n",
    "    # load in the geopandas file with the g scores\n",
    "    data_g_scores = gpd.read_file(\"output/{}_weekly_g_score.geoJSON\".format(re.sub(r'\\W+', '', event_type)))\n",
    "    \n",
    "    # get the data for the specific county these coordinates fall in \n",
    "    county_level_data = data_g_scores[data_g_scores.geometry.contains(Point(coords[1],coords[0]))]\n",
    "    \n",
    "    # not sure whey this is happening\n",
    "    if len(county_level_data) == 0:\n",
    "        return False, 0\n",
    "    \n",
    "    # screwy, but this works. Get an array of values by week\n",
    "    county_data = county_level_data.iloc[:,8:-1].transpose().iloc[:,0].values\n",
    "    \n",
    "    week_of_year = date.isocalendar()[1]\n",
    "    value = pd.Series(county_data).rolling(rolling).mean().values[week_of_year]\n",
    "    \n",
    "    return value > threshold, value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 17.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "tries = 10\n",
    "results = []\n",
    "\n",
    "# check positive events\n",
    "for x in range(tries):\n",
    "    event_type = \"\"\n",
    "    while event_type not in event_types:\n",
    "        event_type, meta_type, coords, date = get_positive_event()\n",
    "    # from https://stackoverflow.com/a/29753985\n",
    "    date = datetime.utcfromtimestamp(date.tolist()/1e9)\n",
    "    try:\n",
    "        prediction = predict(event_type, coords, date)\n",
    "        results.append({\n",
    "            \"truth\":True,\n",
    "            \"prediction\":prediction[0],\n",
    "            \"p_value\":prediction[1],\n",
    "            \"event_type\":event_type,\n",
    "            \"coords\":coords,\n",
    "            \"date\":date\n",
    "        })\n",
    "    except:\n",
    "        print(\"ERR (POS):\",event_type, coords, date)\n",
    "\n",
    "# check negative events\n",
    "for x in range(tries):\n",
    "    event_type = \"\"\n",
    "    while event_type not in event_types:\n",
    "        event_type, meta_type, coords, date = get_negative_event()\n",
    "    try:\n",
    "        prediction = predict(event_type, coords, date)\n",
    "        results.append({\n",
    "            \"truth\":False,\n",
    "            \"prediction\":prediction[0],\n",
    "            \"p_value\":prediction[1],\n",
    "            \"event_type\":event_type,\n",
    "            \"coords\":coords,\n",
    "            \"date\":date\n",
    "        })\n",
    "    except:\n",
    "        print(\"ERR (NEG):\",event_type, coords, date)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coords</th>\n",
       "      <th>date</th>\n",
       "      <th>event_type</th>\n",
       "      <th>p_value</th>\n",
       "      <th>prediction</th>\n",
       "      <th>truth</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>(44.288999999999994, -105.3636)</td>\n",
       "      <td>2016-08-07 18:55:00</td>\n",
       "      <td>Tornado</td>\n",
       "      <td>1.164462</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>(39.22, -94.94)</td>\n",
       "      <td>2017-06-16 23:18:00</td>\n",
       "      <td>Thunderstorm Wind</td>\n",
       "      <td>0.485024</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>(36.29, -91.57)</td>\n",
       "      <td>2017-03-01 03:05:00</td>\n",
       "      <td>Thunderstorm Wind</td>\n",
       "      <td>0.064905</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>(48.01, -98.0)</td>\n",
       "      <td>2016-06-10 02:00:00</td>\n",
       "      <td>Thunderstorm Wind</td>\n",
       "      <td>-0.451598</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>(40.26, -79.18)</td>\n",
       "      <td>2016-06-03 08:26:00</td>\n",
       "      <td>Flash Flood</td>\n",
       "      <td>0.239996</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>(37.939479999999996, -76.94187099999999)</td>\n",
       "      <td>2016-04-02 23:00:00</td>\n",
       "      <td>High Wind</td>\n",
       "      <td>0.036513</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>(38.14, -93.75)</td>\n",
       "      <td>2017-08-20 17:25:00</td>\n",
       "      <td>Thunderstorm Wind</td>\n",
       "      <td>-0.003817</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>(39.0476, -83.0712)</td>\n",
       "      <td>2017-05-27 18:15:00</td>\n",
       "      <td>Flash Flood</td>\n",
       "      <td>0.503592</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>(40.61, -102.47)</td>\n",
       "      <td>2016-06-27 16:30:00</td>\n",
       "      <td>Tornado</td>\n",
       "      <td>1.944154</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>(38.9, -94.37)</td>\n",
       "      <td>2016-05-26 14:32:00</td>\n",
       "      <td>Hail</td>\n",
       "      <td>0.423079</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>(49.16748863798082, -87.3481398519664)</td>\n",
       "      <td>2016-09-10 00:00:00</td>\n",
       "      <td>Winter Storm</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>(32.59063113488605, -68.3310720957892)</td>\n",
       "      <td>2016-10-13 00:00:00</td>\n",
       "      <td>Ice Storm</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>(24.903452881324505, -101.40765439450563)</td>\n",
       "      <td>2016-11-17 00:00:00</td>\n",
       "      <td>Tornado</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>(25.661652087001055, -90.19885645193351)</td>\n",
       "      <td>2016-09-11 00:00:00</td>\n",
       "      <td>Heavy Rain</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>(36.25214051698068, -87.82126043135372)</td>\n",
       "      <td>2016-11-08 00:00:00</td>\n",
       "      <td>Winter Storm</td>\n",
       "      <td>-0.206816</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>(48.36028341662422, -87.30865525580604)</td>\n",
       "      <td>2016-01-06 00:00:00</td>\n",
       "      <td>Excessive Heat</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>(30.538960204993828, -84.3035844545294)</td>\n",
       "      <td>2016-11-08 00:00:00</td>\n",
       "      <td>Ice Storm</td>\n",
       "      <td>-0.092268</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>(37.73191305104431, -81.41733360871362)</td>\n",
       "      <td>2016-04-20 00:00:00</td>\n",
       "      <td>Ice Storm</td>\n",
       "      <td>-0.011516</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>(47.31973834891245, -70.13225918619943)</td>\n",
       "      <td>2016-05-21 00:00:00</td>\n",
       "      <td>High Wind</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>(36.40887967930723, -69.77530881319103)</td>\n",
       "      <td>2016-03-09 00:00:00</td>\n",
       "      <td>Heavy Snow</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                       coords                date  \\\n",
       "0             (44.288999999999994, -105.3636) 2016-08-07 18:55:00   \n",
       "1                             (39.22, -94.94) 2017-06-16 23:18:00   \n",
       "2                             (36.29, -91.57) 2017-03-01 03:05:00   \n",
       "3                              (48.01, -98.0) 2016-06-10 02:00:00   \n",
       "4                             (40.26, -79.18) 2016-06-03 08:26:00   \n",
       "5    (37.939479999999996, -76.94187099999999) 2016-04-02 23:00:00   \n",
       "6                             (38.14, -93.75) 2017-08-20 17:25:00   \n",
       "7                         (39.0476, -83.0712) 2017-05-27 18:15:00   \n",
       "8                            (40.61, -102.47) 2016-06-27 16:30:00   \n",
       "9                              (38.9, -94.37) 2016-05-26 14:32:00   \n",
       "10     (49.16748863798082, -87.3481398519664) 2016-09-10 00:00:00   \n",
       "11     (32.59063113488605, -68.3310720957892) 2016-10-13 00:00:00   \n",
       "12  (24.903452881324505, -101.40765439450563) 2016-11-17 00:00:00   \n",
       "13   (25.661652087001055, -90.19885645193351) 2016-09-11 00:00:00   \n",
       "14    (36.25214051698068, -87.82126043135372) 2016-11-08 00:00:00   \n",
       "15    (48.36028341662422, -87.30865525580604) 2016-01-06 00:00:00   \n",
       "16    (30.538960204993828, -84.3035844545294) 2016-11-08 00:00:00   \n",
       "17    (37.73191305104431, -81.41733360871362) 2016-04-20 00:00:00   \n",
       "18    (47.31973834891245, -70.13225918619943) 2016-05-21 00:00:00   \n",
       "19    (36.40887967930723, -69.77530881319103) 2016-03-09 00:00:00   \n",
       "\n",
       "           event_type   p_value  prediction  truth  \n",
       "0             Tornado  1.164462        True   True  \n",
       "1   Thunderstorm Wind  0.485024        True   True  \n",
       "2   Thunderstorm Wind  0.064905        True   True  \n",
       "3   Thunderstorm Wind -0.451598       False   True  \n",
       "4         Flash Flood  0.239996        True   True  \n",
       "5           High Wind  0.036513        True   True  \n",
       "6   Thunderstorm Wind -0.003817       False   True  \n",
       "7         Flash Flood  0.503592        True   True  \n",
       "8             Tornado  1.944154        True   True  \n",
       "9                Hail  0.423079        True   True  \n",
       "10       Winter Storm  0.000000       False  False  \n",
       "11          Ice Storm  0.000000       False  False  \n",
       "12            Tornado  0.000000       False  False  \n",
       "13         Heavy Rain  0.000000       False  False  \n",
       "14       Winter Storm -0.206816       False  False  \n",
       "15     Excessive Heat  0.000000       False  False  \n",
       "16          Ice Storm -0.092268       False  False  \n",
       "17          Ice Storm -0.011516       False  False  \n",
       "18          High Wind  0.000000       False  False  \n",
       "19         Heavy Snow  0.000000       False  False  "
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DataFrame(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prediction Approach #2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "day_range=21\n",
    "max_range=1.5\n",
    "threshold = 0.75\n",
    "\n",
    "results = []\n",
    "\n",
    "conf={\"tt\":0, \"tf\":0, \"ft\":0, \"ff\":0}\n",
    "\n",
    "for x in range(3):\n",
    "    lat = random.uniform(24.7433195, 49.3457868)\n",
    "    lon = random.uniform(-124.7844079, -66.9513812)\n",
    "    correct = 0\n",
    "    print(x)\n",
    "\n",
    "    \n",
    "    for event_type in event_types:\n",
    "        try:\n",
    "            likely = True\n",
    "\n",
    "            data_g_scores = gpd.read_file(\"output/{}_weekly_g_score.geoJSON\".format(re.sub(r'\\W+', '', event_type)))\n",
    "\n",
    "            # get the data for the county from this data set\n",
    "            county_level_data = data_g_scores[data_g_scores.geometry.contains(Point(lon,lat))]\n",
    "\n",
    "            # screwy, but gets us a week number for our predictions\n",
    "            county_level_data = county_level_data.iloc[:,9:-1].transpose()\n",
    "            county_level_data.columns = ['value']\n",
    "            county_level_data.reset_index(inplace=True)\n",
    "            predicted_weeks = county_level_data.sort_values(\"value\", ascending=False).head(2).index.values\n",
    "\n",
    "            max_like = county_level_data.sort_values(\"value\", ascending=False).head(1).value.values[0]\n",
    "            if max_like < threshold:\n",
    "                likely=False\n",
    "\n",
    "            # derive a date from the sunday of that week we can use to search\n",
    "            search_date = datetime.strptime(\"2017-W{}\".format(predicted_weeks[0]) + '-0', \"%Y-W%W-%w\")\n",
    "\n",
    "            # filter the detail data down by distance and date to see if any of this kind of event are within our window\n",
    "            prediction_set = detail_data[\n",
    "                        ((detail_data.BEGIN_LON > lon-max_range) &\n",
    "                        (detail_data.BEGIN_LON < lon+max_range) &\n",
    "                        (detail_data.BEGIN_LAT > lat-max_range) &\n",
    "                        (detail_data.BEGIN_LAT < lat+max_range)) &\n",
    "                        (detail_data.EVENT_TYPE==event_type)\n",
    "                    ].copy()\n",
    "\n",
    "            prediction_set = prediction_set[\n",
    "                (prediction_set.BEGIN_DATE >= (search_date - timedelta(days=day_range))) &\n",
    "                (prediction_set.BEGIN_DATE <= (search_date + timedelta(days=day_range)))\n",
    "            ]\n",
    "\n",
    "            if (len(prediction_set) > 0 and likely) or (len(prediction_set) == 0 and not likely):\n",
    "                if len(prediction_set) > 0 and likely: conf['tt'] += 1\n",
    "                else: conf['ff'] += 1\n",
    "                    \n",
    "                correct += 1\n",
    "                print(\"Succcess - {} near {} on {} (Guess:{})\".format(event_type, (lat,lon), search_date, likely))\n",
    "            else:\n",
    "                if likely: conf['tf'] += 1\n",
    "                else: conf['ft'] += 1\n",
    "                    \n",
    "                print(\"Fail - {} near {} on {} (Guess:{})\".format(event_type, (lat,lon), search_date, likely))\n",
    "                pass\n",
    "        except:\n",
    "            continue\n",
    "\n",
    "    print(\"Total Correct:{}/{}\".format(correct,len(event_types)))\n",
    "    results.append({\n",
    "        \"lat\":lat,\n",
    "        \"lon\":lon,\n",
    "        \"correct\":correct\n",
    "    })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Old Code I don't want to throw away\n",
    "But I probably should"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get a matrix based on the start and end dates provided\n",
    "# can also specift a meta type or an event type\n",
    "def get_matrix(start, end, meta_type=None, event_type=None):\n",
    "    if meta_type is None and event_type is None:\n",
    "        # filter the detail data by those two dates and add it to the set\n",
    "        data = detail_data[\n",
    "            (detail_data.BEGIN_DATE>pd.to_datetime(start)) &\n",
    "            (detail_data.BEGIN_DATE<pd.to_datetime(end))\n",
    "        ]\n",
    "    else:\n",
    "        if meta_type is not None:\n",
    "            data = detail_data[\n",
    "                (detail_data.BEGIN_DATE>pd.to_datetime(start)) &\n",
    "                (detail_data.BEGIN_DATE<pd.to_datetime(end)) &\n",
    "                (detail_data.META_TYPE == meta_type)\n",
    "            ]\n",
    "        elif event_type is not None:\n",
    "            data = detail_data[\n",
    "                (detail_data.BEGIN_DATE>pd.to_datetime(start)) &\n",
    "                (detail_data.BEGIN_DATE<pd.to_datetime(end)) &\n",
    "                (detail_data.EVENT_TYPE == event_type)\n",
    "            ]\n",
    "\n",
    "    m = lil_matrix((5783,2760), dtype=np.float16)\n",
    "\n",
    "    for i, row in data.iterrows():\n",
    "        mat_val = 1\n",
    "\n",
    "\n",
    "        # is this a single point event or a line\n",
    "        if pd.isnull(row.END_LAT) or pd.isnull(row.END_LON):\n",
    "            # single point - plot a square\n",
    "            col_id = lat_to_index(row.BEGIN_LAT)\n",
    "            row_id = lon_to_index(row.BEGIN_LON)\n",
    "            for r in range(row_id-3, row_id+4):\n",
    "                for c in range(col_id-3, col_id+4):\n",
    "                    try:\n",
    "                        m[r,c] += mat_val\n",
    "                    except:\n",
    "                        print(\"bad row:\",row)\n",
    "        else:\n",
    "            # we have a start and an end, so it's a line\n",
    "            start = (lat_to_index(row.BEGIN_LAT),lon_to_index(row.BEGIN_LON))\n",
    "            end = (lat_to_index(row.END_LAT),lon_to_index(row.END_LON))\n",
    "            for pos in bresenham_line(start,end):\n",
    "                try:\n",
    "                    m[pos[1],pos[0]] += mat_val\n",
    "                except:\n",
    "                    print(\"bad row (line):\",pos)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_year_month_matrix(year, month):\n",
    "    start = \"{:0d}-{:02d}-01\".format(year, month)\n",
    "    if month == 12:\n",
    "        month = 1\n",
    "        year += 1\n",
    "    else:\n",
    "        month+=1\n",
    "    end = \"{:0d}-{:02d}-01\".format(year, month)\n",
    "    year_data = detail_data[(detail_data.BEGIN_DATE>start) & (detail_data.BEGIN_DATE<end)].dropna()\n",
    "    # trim it to just the continental US\n",
    "    year_data = year_data[\n",
    "        (year_data.BEGIN_LON > -124.7844079) &\n",
    "        (year_data.BEGIN_LON < -66.9513812) &\n",
    "        (year_data.BEGIN_LAT > 24.7433195) &\n",
    "        (year_data.BEGIN_LAT < 49.3457868)\n",
    "    ]\n",
    "    m = lil_matrix((5783,2760), dtype=np.int16)\n",
    "    for i, row in year_data.iterrows():\n",
    "        if row.END_LAT==np.nan or row.END_LON==np.nan:\n",
    "            col_id = lat_to_index(row.BEGIN_LAT)\n",
    "            row_id = lon_to_index(row.BEGIN_LON)\n",
    "            for r in range(row_id-3, row_id+4):\n",
    "                for c in range(col_id-3, col_id+4):\n",
    "                    try:\n",
    "                        m[r,c] = 1\n",
    "                    except:\n",
    "                        print(\"bad row:\",row)\n",
    "        else:\n",
    "            # we have a start and an end, so it's a line\n",
    "            start = (lat_to_index(row.BEGIN_LAT),lon_to_index(row.BEGIN_LON))\n",
    "            end = (lat_to_index(row.END_LAT),lon_to_index(row.END_LON))\n",
    "            for pos in bresenham_line(start,end):\n",
    "                try:\n",
    "                    m[pos[1],pos[0]] = 1\n",
    "                except:\n",
    "                    print(\"bad row:\",pos)\n",
    "    return m\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "county_data = get_county_shapemap()\n",
    "county_data.set_index('GEOID', inplace=True)\n",
    "\n",
    "start_date = datetime.strptime('01/01/1950', '%m/%d/%Y')\n",
    "end_date = datetime.strptime('01/01/2019', '%m/%d/%Y')\n",
    "\n",
    "for event_type in [x for x in detail_data.EVENT_TYPE.unique() if x not in [\"Avalanche\", \"Blizzard\", \"Flash Flood\", \"Hail\", \"Thunderstorm Wind\", \"Tornado\", \"Winter Storm\"]]:\n",
    "    print(event_type)\n",
    "    m = get_matrix(start_date, end_date, event_type=event_type)\n",
    "    data = create_geopandas_from_matrix(m)\n",
    "    print(\"Creating plots from {} entries\".format(len(data)))\n",
    "    val = get_spatial_lag_quartiles(\n",
    "        data, \n",
    "        title=\"Spatial Lag Median Risk Quartiles - {}\".format(event_type),\n",
    "        filename=\"event_output/spatial_lag_quartile_{0}.png\".format(event_type.replace('/',' ')),\n",
    "        ret_values=True)\n",
    "    county_data[event_type+\"_quartile\"] = val.cl\n",
    "    val = get_gettis_ord_g(\n",
    "        data, \n",
    "        title=\"Gettis Ord G* Autocorrelation - {}\".format(event_type), \n",
    "        filename=\"event_output/gettis_ord_g_{0}.png\".format(event_type.replace('/',' ')),\n",
    "        ret_values=True)\n",
    "    county_data[event_type+\"_g_local\"] = val.counts_g_local\n",
    "    county_data.to_csv(\"county_overall_output.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "meta_type_colors = {\n",
    "    'Cold':\"Blues\",\n",
    "    'All':\"Greys\",\n",
    "    'Wind':\"BuGn\",\n",
    "    'Water':\"PuBu\",\n",
    "    'Other':\"Greens\",\n",
    "    'Heat':\"YlOrRd\"\n",
    "}\n",
    "\n",
    "county_data = get_county_shapemap()\n",
    "county_data.set_index(\"GEOID\", inplace=True)\n",
    "\n",
    "for meta_type, color_map in meta_type_colors.items():\n",
    "    date = datetime.strptime('01/01/18', '%m/%d/%y')\n",
    "    while date.year<2019:\n",
    "        print(\"Plotting {} - {}\".format(meta_type, str(date).split()[0]))\n",
    "        if meta_type == \"All\":\n",
    "            m = get_day_matrix(date, span=14)\n",
    "        else:\n",
    "            m = get_day_matrix(date, span=14, meta_type=meta_type)\n",
    "        data = create_geopandas_from_matrix(m)\n",
    "        if data is None:\n",
    "            print(\"No data found for\", date)\n",
    "        else:\n",
    "            val = get_gettis_ord_g(\n",
    "                data,\n",
    "                title=\"Gettis Ord G* Autocorrelation - {} - ({})\".format(meta_type, str(date).split()[0]), \n",
    "                filename=\"temp.png\",\n",
    "                color_map=color_map,\n",
    "                ret_values=True\n",
    "            )\n",
    "            county_data[\"{}-{:%m-%d}_p_sim\".format(meta_type, date)] = val.counts_p_sim\n",
    "            county_data.to_csv(\"weekly_p_sim.csv\")\n",
    "        date = date + timedelta(days=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "county_data = get_county_shapemap()\n",
    "county_data.set_index(\"GEOID\", inplace=True)\n",
    "\n",
    "while date.year<2019:\n",
    "    print(\"Plotting {} - {}\".format(meta_type, str(date).split()[0]))\n",
    "    m = get_day_matrix(date, span=14, event_type='Tornado')\n",
    "    data = create_geopandas_from_matrix(m)\n",
    "    if data is None:\n",
    "        print(\"No data found for\", date)\n",
    "    else:\n",
    "        val = get_gettis_ord_g(\n",
    "            data,\n",
    "            title=\"Gettis Ord G* Autocorrelation - {} - ({})\".format(meta_type, str(date).split()[0]), \n",
    "            filename=\"temp.png\",\n",
    "            color_map=color_map,\n",
    "            ret_values=True\n",
    "        )\n",
    "        county_data[\"{}-{:%m-%d}_p_sim\".format(meta_type, date)] = val.counts_p_sim\n",
    "        county_data.to_csv(\"tornado_p_sim.csv\")\n",
    "    date = date + timedelta(days=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drought_p_sim = pd.read_csv(\"tornado_p_sim.csv\")\n",
    "drought_p_sim.set_index(\"GEOID\", inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drought_p_sim.loc[48221][9:].astype(float).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_id = 29095\n",
    "rolling = 4\n",
    "date = datetime.strptime('01/01/18', '%m/%d/%y')\n",
    "date_list = [\"{:%m-%d}\".format(date+timedelta(days=x)) for x in range(0, 365, 7)]\n",
    "\n",
    "f, ax = plt.subplots(1, figsize=(20,5))\n",
    "sns.lineplot(\n",
    "    y=running_mean(drought_p_sim.loc[geo_id][9:].astype(float).values, rolling), \n",
    "    x=date_list[:-rolling+1], palette=\"RdGy\")\n",
    "\n",
    "# ax.legend(['Drought Likelyhood'])\n",
    "\n",
    "plt.xticks(rotation=90)\n",
    "plt.title(\"Tornado Likelyhood for 39.0997° N, 94.5786° W (Kansas City, MO)\")\n",
    "plt.savefig(\"kc_plot.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weekly_p_sim = pd.read_csv(\"weekly_p_sim.csv\")\n",
    "weekly_p_sim.set_index(\"GEOID\", inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cold_p_sim = weekly_p_sim.iloc[:,9:62].copy()\n",
    "all_p_sim = weekly_p_sim.iloc[:,62:115].copy()\n",
    "wind_p_sim = weekly_p_sim.iloc[:,115:168].copy()\n",
    "water_p_sim = weekly_p_sim.iloc[:,168:221].copy()\n",
    "other_p_sim = weekly_p_sim.iloc[:,221:274].copy()\n",
    "heat_p_sim = weekly_p_sim.iloc[:,274:327].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set()\n",
    "geo_id = 26161\n",
    "rolling = 3\n",
    "date = datetime.strptime('01/01/18', '%m/%d/%y')\n",
    "date_list = [\"{:%m-%d}\".format(date+timedelta(days=x)) for x in range(0, 365, 7)]\n",
    "\n",
    "f, ax = plt.subplots(1, figsize=(20,5))\n",
    "sns.lineplot(y=running_mean(cold_p_sim.loc[geo_id].values, rolling), x=date_list[:-rolling+1])\n",
    "sns.lineplot(y=running_mean(wind_p_sim.loc[geo_id].values, rolling), x=date_list[:-rolling+1])\n",
    "sns.lineplot(y=running_mean(water_p_sim.loc[geo_id].values, rolling), x=date_list[:-rolling+1])\n",
    "sns.lineplot(y=running_mean(heat_p_sim.loc[geo_id].values, rolling), x=date_list[:-rolling+1])\n",
    "\n",
    "ax.legend(['Cold','Wind','Water','Heat'])\n",
    "plt.xticks(rotation=90)\n",
    "\n",
    "plt.savefig(\"aa_plot.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set()\n",
    "geo_id = 40109\n",
    "rolling = 3\n",
    "date = datetime.strptime('01/01/18', '%m/%d/%y')\n",
    "date_list = [\"{:%m-%d}\".format(date+timedelta(days=x)) for x in range(0, 365, 7)]\n",
    "\n",
    "f, ax = plt.subplots(1, figsize=(20,5))\n",
    "sns.lineplot(y=running_mean(cold_p_sim.loc[geo_id].values, rolling), x=date_list[:-rolling+1])\n",
    "sns.lineplot(y=running_mean(wind_p_sim.loc[geo_id].values, rolling), x=date_list[:-rolling+1])\n",
    "sns.lineplot(y=running_mean(water_p_sim.loc[geo_id].values, rolling), x=date_list[:-rolling+1])\n",
    "sns.lineplot(y=running_mean(heat_p_sim.loc[geo_id].values, rolling), x=date_list[:-rolling+1])\n",
    "\n",
    "ax.legend(['Cold','Wind','Water','Heat'])\n",
    "plt.xticks(rotation=90)\n",
    "plt.title(\"Annual Weather Type Risks for 35.4676° N, 97.5164° W (Oklahoma City, OK)\")\n",
    "plt.savefig(\"ok_plot.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rolling = 5\n",
    "date = datetime.strptime('01/01/18', '%m/%d/%y')\n",
    "date_list = [date+timedelta(days=x) for x in range(0, 365, 7)][:-rolling+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(geo_id, date):\n",
    "    date = date.replace(year=2018)\n",
    "    p_sims = DataFrame({\n",
    "        \"Cold\":running_mean(cold_p_sim.loc[geo_id].values, rolling), \n",
    "        \"Wind\":running_mean(wind_p_sim.loc[geo_id].values, rolling),\n",
    "        \"Water\":running_mean(water_p_sim.loc[geo_id].values, rolling),\n",
    "        \"Heat\":running_mean(heat_p_sim.loc[geo_id].values, rolling),\n",
    "        \"date\":date_list})\n",
    "\n",
    "    start = date - timedelta(days=4)\n",
    "    end = date + timedelta(days=4)\n",
    "    predictions = p_sims[\n",
    "        (p_sims.date >= start) &\n",
    "        (p_sims.date <= end)\n",
    "    ].mean().sort_values(ascending=False)\n",
    "    \n",
    "    return predictions\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "date = datetime.strptime('05/24/18', '%m/%d/%y')\n",
    "geo_id = 36111\n",
    "\n",
    "predictions = predict(geo_id, date)"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}