smoke_gets_in_your_eyes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Smoke Gets in Your Eyes\n",
    "\n",
    "This notebook contains the essential final Python code for the results of the _Smoke Gets in Your Eyes_ miniproject."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sqlite3\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn import tree, preprocessing\n",
    "import sklearn.ensemble as ske\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LinearRegression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocess wildfire data for severity prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   FIRE_YEAR DISCOVERY_TIME STAT_CAUSE_DESCR  CONT_DATE CONT_TIME   LATITUDE  \\\n",
      "0       2005           1300    Miscellaneous  2453403.5      1730  40.036944   \n",
      "1       2004           0845        Lightning  2453137.5      1530  38.933056   \n",
      "2       2004           1921   Debris Burning  2453156.5      2024  38.984167   \n",
      "3       2004           1600        Lightning  2453189.5      1400  38.559167   \n",
      "4       2004           1600        Lightning  2453189.5      1200  38.559167   \n",
      "\n",
      "    LONGITUDE STATE  DISCOVERY_DATE  FIRE_SIZE FIRE_SIZE_CLASS  \n",
      "0 -121.005833    CA       2453403.5       0.10               A  \n",
      "1 -120.404444    CA       2453137.5       0.25               A  \n",
      "2 -120.735556    CA       2453156.5       0.10               A  \n",
      "3 -119.913333    CA       2453184.5       0.10               A  \n",
      "4 -119.933056    CA       2453184.5       0.10               A  \n"
     ]
    }
   ],
   "source": [
    "# Read wildfire data from SQLITE database\n",
    "cnx = sqlite3.connect('FPA_FOD_20170508.sqlite')\n",
    "df = pd.read_sql_query(\"SELECT FIRE_YEAR,DISCOVERY_TIME,STAT_CAUSE_DESCR,CONT_DATE,CONT_TIME,LATITUDE,LONGITUDE,STATE,DISCOVERY_DATE,FIRE_SIZE,FIRE_SIZE_CLASS FROM 'Fires'\", cnx)\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   FIRE_YEAR DISCOVERY_TIME STAT_CAUSE_DESCR  CONT_DATE CONT_TIME   LATITUDE  \\\n",
      "0       2005           1300    Miscellaneous 2005-02-02      1730  40.036944   \n",
      "1       2004           0845        Lightning 2004-05-12      1530  38.933056   \n",
      "2       2004           1921   Debris Burning 2004-05-31      2024  38.984167   \n",
      "3       2004           1600        Lightning 2004-07-03      1400  38.559167   \n",
      "4       2004           1600        Lightning 2004-07-03      1200  38.559167   \n",
      "\n",
      "    LONGITUDE STATE DISCOVERY_DATE  FIRE_SIZE FIRE_SIZE_CLASS  \n",
      "0 -121.005833    CA     2005-02-02       0.10               A  \n",
      "1 -120.404444    CA     2004-05-12       0.25               A  \n",
      "2 -120.735556    CA     2004-05-31       0.10               A  \n",
      "3 -119.913333    CA     2004-06-28       0.10               A  \n",
      "4 -119.933056    CA     2004-06-28       0.10               A  \n"
     ]
    }
   ],
   "source": [
    "# Convert fire discovery and containment dates to yyyy-mm-dd\n",
    "df['DISCOVERY_DATE'] = pd.to_datetime(df['DISCOVERY_DATE'] - pd.Timestamp(0).to_julian_date(), unit='D')\n",
    "df['CONT_DATE'] = pd.to_datetime(df['CONT_DATE'] - pd.Timestamp(0).to_julian_date(), unit='D')\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   FIRE_YEAR DISCOVERY_TIME STAT_CAUSE_DESCR  CONT_DATE CONT_TIME   LATITUDE  \\\n",
      "0       2005           1300    Miscellaneous 2005-02-02      1730  40.036944   \n",
      "1       2004           0845        Lightning 2004-05-12      1530  38.933056   \n",
      "2       2004           1921   Debris Burning 2004-05-31      2024  38.984167   \n",
      "3       2004           1600        Lightning 2004-07-03      1400  38.559167   \n",
      "4       2004           1600        Lightning 2004-07-03      1200  38.559167   \n",
      "\n",
      "    LONGITUDE STATE DISCOVERY_DATE  FIRE_SIZE FIRE_SIZE_CLASS  \\\n",
      "0 -121.005833    CA     2005-02-02       0.10               A   \n",
      "1 -120.404444    CA     2004-05-12       0.25               A   \n",
      "2 -120.735556    CA     2004-05-31       0.10               A   \n",
      "3 -119.913333    CA     2004-06-28       0.10               A   \n",
      "4 -119.933056    CA     2004-06-28       0.10               A   \n",
      "\n",
      "   DISCOVERY_MONTH  DISCOVERY_DAY DISCOVERY_DAY_OF_WEEK  CONT_MONTH  CONT_DAY  \\\n",
      "0                2              2             Wednesday         2.0       2.0   \n",
      "1                5             12             Wednesday         5.0      12.0   \n",
      "2                5             31                Monday         5.0      31.0   \n",
      "3                6             28                Monday         7.0       3.0   \n",
      "4                6             28                Monday         7.0       3.0   \n",
      "\n",
      "  CONT_DAY_OF_WEEK  \n",
      "0        Wednesday  \n",
      "1        Wednesday  \n",
      "2           Monday  \n",
      "3         Saturday  \n",
      "4         Saturday  \n"
     ]
    }
   ],
   "source": [
    "# Add new columns for month, date, and weekday for discovery and containment dates\n",
    "\n",
    "df['DISCOVERY_MONTH'] = pd.DatetimeIndex(df['DISCOVERY_DATE']).month\n",
    "df['DISCOVERY_DAY'] = pd.DatetimeIndex(df['DISCOVERY_DATE']).day\n",
    "df['DISCOVERY_DAY_OF_WEEK'] = df['DISCOVERY_DATE'].dt.weekday_name\n",
    "\n",
    "#df['CONT_DATE'].fillna(0)\n",
    "df['CONT_MONTH'] = pd.DatetimeIndex(df['CONT_DATE']).month\n",
    "df['CONT_DAY'] = pd.DatetimeIndex(df['CONT_DATE']).day\n",
    "df['CONT_DAY_OF_WEEK'] = df['CONT_DATE'].dt.weekday_name\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   FIRE_YEAR DISCOVERY_TIME  STAT_CAUSE_DESCR  CONT_DATE CONT_TIME   LATITUDE  \\\n",
      "0       2005           1300                 7 2005-02-02      1730  40.036944   \n",
      "1       2004           0845                 6 2004-05-12      1530  38.933056   \n",
      "2       2004           1921                 3 2004-05-31      2024  38.984167   \n",
      "3       2004           1600                 6 2004-07-03      1400  38.559167   \n",
      "4       2004           1600                 6 2004-07-03      1200  38.559167   \n",
      "\n",
      "    LONGITUDE  STATE DISCOVERY_DATE  FIRE_SIZE FIRE_SIZE_CLASS  \\\n",
      "0 -121.005833      4     2005-02-02       0.10               A   \n",
      "1 -120.404444      4     2004-05-12       0.25               A   \n",
      "2 -120.735556      4     2004-05-31       0.10               A   \n",
      "3 -119.913333      4     2004-06-28       0.10               A   \n",
      "4 -119.933056      4     2004-06-28       0.10               A   \n",
      "\n",
      "   DISCOVERY_MONTH  DISCOVERY_DAY  DISCOVERY_DAY_OF_WEEK  CONT_MONTH  \\\n",
      "0                2              2                      6         2.0   \n",
      "1                5             12                      6         5.0   \n",
      "2                5             31                      1         5.0   \n",
      "3                6             28                      1         7.0   \n",
      "4                6             28                      1         7.0   \n",
      "\n",
      "   CONT_DAY CONT_DAY_OF_WEEK  \n",
      "0       2.0        Wednesday  \n",
      "1      12.0        Wednesday  \n",
      "2      31.0           Monday  \n",
      "3       3.0         Saturday  \n",
      "4       3.0         Saturday  \n"
     ]
    }
   ],
   "source": [
    "le = preprocessing.LabelEncoder()\n",
    "df['STAT_CAUSE_DESCR'] = le.fit_transform(df['STAT_CAUSE_DESCR'])\n",
    "df['STATE'] = le.fit_transform(df['STATE'])\n",
    "df['DISCOVERY_DAY_OF_WEEK'] = le.fit_transform(df['DISCOVERY_DAY_OF_WEEK'])\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   FIRE_YEAR DISCOVERY_TIME  STAT_CAUSE_DESCR  CONT_DATE CONT_TIME   LATITUDE  \\\n",
      "0       2005           1300                 7 2005-02-02      1730  40.036944   \n",
      "1       2004           0845                 6 2004-05-12      1530  38.933056   \n",
      "2       2004           1921                 3 2004-05-31      2024  38.984167   \n",
      "3       2004           1600                 6 2004-07-03      1400  38.559167   \n",
      "4       2004           1600                 6 2004-07-03      1200  38.559167   \n",
      "\n",
      "    LONGITUDE  STATE DISCOVERY_DATE  FIRE_SIZE  FIRE_SIZE_CLASS  \\\n",
      "0 -121.005833      4     2005-02-02       0.10                0   \n",
      "1 -120.404444      4     2004-05-12       0.25                0   \n",
      "2 -120.735556      4     2004-05-31       0.10                0   \n",
      "3 -119.913333      4     2004-06-28       0.10                0   \n",
      "4 -119.933056      4     2004-06-28       0.10                0   \n",
      "\n",
      "   DISCOVERY_MONTH  DISCOVERY_DAY  DISCOVERY_DAY_OF_WEEK CONT_MONTH CONT_DAY  \\\n",
      "0                2              2                      6          2        2   \n",
      "1                5             12                      6          5       12   \n",
      "2                5             31                      1          5       31   \n",
      "3                6             28                      1          7        3   \n",
      "4                6             28                      1          7        3   \n",
      "\n",
      "   CONT_DAY_OF_WEEK  \n",
      "0                 7  \n",
      "1                 7  \n",
      "2                 1  \n",
      "3                 2  \n",
      "4                 2  \n"
     ]
    }
   ],
   "source": [
    "df['CONT_DAY_OF_WEEK']=df['CONT_DAY_OF_WEEK'].fillna(\"Unknown\")\n",
    "df['CONT_DAY']=df['CONT_DAY'].fillna(\"0\")\n",
    "df['CONT_MONTH']=df['CONT_MONTH'].fillna(\"0\")\n",
    "df['CONT_TIME']=df['CONT_TIME'].fillna(\"0\")\n",
    "df['DISCOVERY_TIME']=df['DISCOVERY_TIME'].fillna(\"0\")\n",
    "\n",
    "df['CONT_DAY_OF_WEEK'] = le.fit_transform(df['CONT_DAY_OF_WEEK'])\n",
    "df['FIRE_SIZE_CLASS'] = le.fit_transform(df['FIRE_SIZE_CLASS'])\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['CONT_MONTH']=df['CONT_MONTH'].astype('Float64')\n",
    "df['CONT_DAY']=df['CONT_DAY'].astype('Float64')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ML for predicting wildfire severity (size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow\n",
    "import keras\n",
    "from keras.models import Sequential, load_model\n",
    "from keras.layers import Dense, Activation, BatchNormalization\n",
    "from keras.optimizers import SGD\n",
    "from keras.activations import relu\n",
    "from keras.utils import plot_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0\n",
       "1    0\n",
       "2    0\n",
       "3    0\n",
       "4    0\n",
       "Name: FIRE_SIZE_CLASS, dtype: int64"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "labels=df['FIRE_SIZE_CLASS']\n",
    "labels.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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>FIRE_YEAR</th>\n",
       "      <th>STAT_CAUSE_DESCR</th>\n",
       "      <th>LATITUDE</th>\n",
       "      <th>LONGITUDE</th>\n",
       "      <th>DISCOVERY_MONTH</th>\n",
       "      <th>DISCOVERY_DAY</th>\n",
       "      <th>DISCOVERY_DAY_OF_WEEK</th>\n",
       "      <th>CONT_MONTH</th>\n",
       "      <th>CONT_DAY</th>\n",
       "      <th>CONT_DAY_OF_WEEK</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2005</td>\n",
       "      <td>7</td>\n",
       "      <td>40.036944</td>\n",
       "      <td>-121.005833</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>6</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2004</td>\n",
       "      <td>6</td>\n",
       "      <td>38.933056</td>\n",
       "      <td>-120.404444</td>\n",
       "      <td>5</td>\n",
       "      <td>12</td>\n",
       "      <td>6</td>\n",
       "      <td>5.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2004</td>\n",
       "      <td>3</td>\n",
       "      <td>38.984167</td>\n",
       "      <td>-120.735556</td>\n",
       "      <td>5</td>\n",
       "      <td>31</td>\n",
       "      <td>1</td>\n",
       "      <td>5.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2004</td>\n",
       "      <td>6</td>\n",
       "      <td>38.559167</td>\n",
       "      <td>-119.913333</td>\n",
       "      <td>6</td>\n",
       "      <td>28</td>\n",
       "      <td>1</td>\n",
       "      <td>7.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2004</td>\n",
       "      <td>6</td>\n",
       "      <td>38.559167</td>\n",
       "      <td>-119.933056</td>\n",
       "      <td>6</td>\n",
       "      <td>28</td>\n",
       "      <td>1</td>\n",
       "      <td>7.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   FIRE_YEAR  STAT_CAUSE_DESCR   LATITUDE   LONGITUDE  DISCOVERY_MONTH  \\\n",
       "0       2005                 7  40.036944 -121.005833                2   \n",
       "1       2004                 6  38.933056 -120.404444                5   \n",
       "2       2004                 3  38.984167 -120.735556                5   \n",
       "3       2004                 6  38.559167 -119.913333                6   \n",
       "4       2004                 6  38.559167 -119.933056                6   \n",
       "\n",
       "   DISCOVERY_DAY  DISCOVERY_DAY_OF_WEEK  CONT_MONTH  CONT_DAY  \\\n",
       "0              2                      6         2.0       2.0   \n",
       "1             12                      6         5.0      12.0   \n",
       "2             31                      1         5.0      31.0   \n",
       "3             28                      1         7.0       3.0   \n",
       "4             28                      1         7.0       3.0   \n",
       "\n",
       "   CONT_DAY_OF_WEEK  \n",
       "0                 7  \n",
       "1                 7  \n",
       "2                 1  \n",
       "3                 2  \n",
       "4                 2  "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logits=df.drop(['FIRE_SIZE','FIRE_SIZE_CLASS','DISCOVERY_DATE','CONT_DATE','STATE','CONT_TIME','DISCOVERY_TIME'],axis=1)\n",
    "logits.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels_cat=keras.utils.to_categorical(labels,7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train, x_test, y_train, y_test = train_test_split(logits,labels_cat,test_size=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels_cat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model=Sequential()\n",
    "model.add(BatchNormalization(input_shape=[10]))\n",
    "model.add(Dense(500,kernel_initializer='truncated_normal'))\n",
    "model.add(Activation('sigmoid'))\n",
    "model.add(Dense(300,kernel_initializer='truncated_normal'))\n",
    "model.add(Activation('sigmoid'))\n",
    "model.add(Dense(7,kernel_initializer='truncated_normal'))\n",
    "model.add(Activation('softmax'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.compile(optimizer=SGD(lr=0.01, momentum=0.9,decay=1e-6),\n",
    "              loss='categorical_crossentropy',\n",
    "              metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.fit(x_train,y_train,epochs=10,validation_data=(np.array(x_test),np.array(y_test)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "model=load_model('64%model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.evaluate(x_test,y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=model.predict(x_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_model(model, to_file='network.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "batch_normalization_67 (Batc (None, 10)                40        \n",
      "_________________________________________________________________\n",
      "dense_82 (Dense)             (None, 1000)              11000     \n",
      "_________________________________________________________________\n",
      "activation_81 (Activation)   (None, 1000)              0         \n",
      "_________________________________________________________________\n",
      "batch_normalization_68 (Batc (None, 1000)              4000      \n",
      "_________________________________________________________________\n",
      "dropout_33 (Dropout)         (None, 1000)              0         \n",
      "_________________________________________________________________\n",
      "dense_83 (Dense)             (None, 1200)              1201200   \n",
      "_________________________________________________________________\n",
      "activation_82 (Activation)   (None, 1200)              0         \n",
      "_________________________________________________________________\n",
      "batch_normalization_69 (Batc (None, 1200)              4800      \n",
      "_________________________________________________________________\n",
      "dropout_34 (Dropout)         (None, 1200)              0         \n",
      "_________________________________________________________________\n",
      "dense_84 (Dense)             (None, 1200)              1441200   \n",
      "_________________________________________________________________\n",
      "activation_83 (Activation)   (None, 1200)              0         \n",
      "_________________________________________________________________\n",
      "batch_normalization_70 (Batc (None, 1200)              4800      \n",
      "_________________________________________________________________\n",
      "dropout_35 (Dropout)         (None, 1200)              0         \n",
      "_________________________________________________________________\n",
      "dense_85 (Dense)             (None, 500)               600500    \n",
      "_________________________________________________________________\n",
      "activation_84 (Activation)   (None, 500)               0         \n",
      "_________________________________________________________________\n",
      "batch_normalization_71 (Batc (None, 500)               2000      \n",
      "_________________________________________________________________\n",
      "dropout_36 (Dropout)         (None, 500)               0         \n",
      "_________________________________________________________________\n",
      "dense_86 (Dense)             (None, 256)               128256    \n",
      "_________________________________________________________________\n",
      "activation_85 (Activation)   (None, 256)               0         \n",
      "_________________________________________________________________\n",
      "batch_normalization_72 (Batc (None, 256)               1024      \n",
      "_________________________________________________________________\n",
      "dense_87 (Dense)             (None, 7)                 1799      \n",
      "_________________________________________________________________\n",
      "activation_86 (Activation)   (None, 7)                 0         \n",
      "=================================================================\n",
      "Total params: 3,400,619\n",
      "Trainable params: 3,392,287\n",
      "Non-trainable params: 8,332\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating Air quality data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def preprocess_site_list(site_file):\n",
    "    \"\"\"Get the location of air quality measurement sites.\n",
    "        \n",
    "    Return the latitude and longitude coordinates in two dictionaries. \n",
    "    The keys are site codes as strings and the values are coordinates as floats.\n",
    "    \"\"\"\n",
    "\n",
    "    # Make sure that state, county, and site codes are read as strings\n",
    "    # instead of numbers\n",
    "    col_types = {'State Code': str, 'County Code': str, 'Site Number': str}\n",
    "\n",
    "    sites = pd.read_csv(site_file, dtype=col_types)\n",
    "\n",
    "    # Columns that can be dropped from the site listing\n",
    "    cols_to_drop = [\n",
    "        'Land Use', 'Location Setting', 'Met Site State Code',\n",
    "        'Met Site County Code', 'Met Site Site Number', 'Met Site Type',\n",
    "        'Met Site Distance', 'Met Site Direction', 'GMT Offset',\n",
    "        'Owning Agency', 'Local Site Name', \"Address\", \"Zip Code\",\n",
    "        \"State Name\", \"County Name\", \"City Name\", \"CBSA Name\", \"Tribe Name\",\n",
    "        \"Extraction Date\"\n",
    "    ]\n",
    "\n",
    "    sites = sites.drop(columns=cols_to_drop)\n",
    "\n",
    "    # Convert site closing dates to datetime64 objects (does not work\n",
    "    # with read_csv())\n",
    "    sites['Site Closed Date'] = pd.to_datetime(sites['Site Closed Date'])\n",
    "\n",
    "    # Get sites that have been closed down before 1992 and drop them\n",
    "    # from the list\n",
    "    outdated = sites.loc[\n",
    "        sites['Site Closed Date'] < pd.to_datetime('01-01-1992')]\n",
    "    sites = sites.drop(outdated.index)\n",
    "\n",
    "    # Get the complete site code, add it as a column, and finally\n",
    "    # use the site code as the index\n",
    "    sites['site-code'] = sites['State Code'] + '-' + sites[\n",
    "        'County Code'] + '-' + sites['Site Number']\n",
    "    sites = sites.set_index('site-code')\n",
    "\n",
    "    # Divide the data into two Series: one with latitude coordinates and one\n",
    "    # with longitude coordinates\n",
    "    sites_lat = sites['Latitude']\n",
    "    sites_lon = sites['Longitude']\n",
    "\n",
    "    # Transform the site listings into a dictionary where keys are site\n",
    "    # codes and values are coordinates}\n",
    "    sites_lat = sites_lat.to_dict()\n",
    "    sites_lon = sites_lon.to_dict()\n",
    "\n",
    "    # Return the two dictionaries with the coordinates\n",
    "    return sites_lat, sites_lon\n",
    "\n",
    "\n",
    "def preprocess_aqi(sites_lat, sites_lon, aqi_file):\n",
    "    \"\"\"Convert a CSV file with AQI data into a pandas DataFrame.\n",
    "\n",
    "    Get the following columns from an AQI data file: Date, AQI, Category,\n",
    "    Defining Parameter, and Defining Site. Ensure that dates are parsed as \n",
    "    datetime objects.\n",
    "\n",
    "    This function also adds coordinate data for the defining site to the data \n",
    "    frame.\n",
    "    \"\"\"\n",
    "\n",
    "    # Specify used column names and types to make CSV parsing more efficient\n",
    "    col_types = {\n",
    "        'AQI': int,\n",
    "        'Category': str,\n",
    "        'Defining Parameter': str,\n",
    "        'Defining Site': str\n",
    "    }\n",
    "\n",
    "    col_names = list(col_types.keys()) + ['Date']\n",
    "\n",
    "    aqi = pd.read_csv(\n",
    "        aqi_file, usecols=col_names, dtype=col_types, parse_dates=['Date'])\n",
    "\n",
    "    # Add column containing the site information (including location\n",
    "    # coordinates) to each AQI measurement\n",
    "    aqi['Latitude'] = aqi['Defining Site'].map(sites_lat)\n",
    "    aqi['Longitude'] = aqi['Defining Site'].map(sites_lon)\n",
    "\n",
    "    return aqi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get dictionaries containing latitude and longitude coordinates for AQI measurement sites\n",
    "sites_lat, sites_lon = preprocess_site_list('data/aqi/aqs_sites.csv')\n",
    "\n",
    "# Generate a list of years 1992–2015\n",
    "years = [str(y) for y in range(1992, 2016)]\n",
    "\n",
    "# Generate a dictionary where keys are years (as strings) and values are DataFrames\n",
    "# containing the AQI data with coordinates added\n",
    "aqi = {\n",
    "    y: preprocess_aqi(sites_lat, sites_lon,\n",
    "                      'data/aqi/daily_aqi_by_county_{}.csv'.format(y))\n",
    "    for y in years\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from pandas.tseries.offsets import DateOffset\n",
    "\n",
    "\n",
    "LON_AT_EQUATOR = 111\n",
    "LAT_IN_KM = 111\n",
    "\n",
    "\n",
    "def get_distance(lat1, lon1, lat2, lon2):\n",
    "    dist = math.sqrt((lat1 - lat2)**2 + (lon1 - lon2)**2)\n",
    "    return round(dist, 6)\n",
    "\n",
    "\n",
    "def dist_as_lon_degree(distance, latitude):\n",
    "    \"\"\"A simplified conversion of a specified distance into degrees longitude\n",
    "    at a specified latitude\"\"\"\n",
    "    lon_deg_as_km = LON_AT_EQUATOR * math.cos(math.radians(latitude))\n",
    "    return distance / lon_deg_as_km\n",
    "\n",
    "\n",
    "def get_aqis_during_fire(year, start_date, end_date):\n",
    "    # Narrow down search by getting measurements from the time of the fire.\n",
    "    # Note that here the date range is defined as DISCOVERY_DATE + 1 to CONT_DATE + 1:\n",
    "    # it seems likely that measurements from the discovery date might have been made before\n",
    "    # the fire was discovered.\n",
    "    # If either of the given dates is a null value, return an empty dataframe\n",
    "    if pd.isnull(start_date) or pd.isnull(end_date):\n",
    "        return pd.DataFrame()\n",
    "    aqi_df = aqi[year]\n",
    "    date_range = pd.date_range(start_date + DateOffset(days=1),\n",
    "                               end_date + DateOffset(days=1))\n",
    "    return aqi_df.loc[aqi_df['Date'].isin(date_range)]\n",
    "\n",
    "\n",
    "def get_nearest_measurement(fire_lat, fire_lon, aqi_df):\n",
    "    # Narrow down the search space by only including measurements that are +/- 0.5 degrees\n",
    "    # from the fire coordinates (this translates to a square area with the fire in the\n",
    "    # center: the square is approx. 110km x 110km\n",
    "    # If the dataframe is empty, return an empty Series\n",
    "    if aqi_df.empty:\n",
    "        return pd.Series()\n",
    "    lon_degrees = dist_as_lon_degree(LAT_IN_KM, fire_lat)\n",
    "    candidates = aqi_df.loc[(\n",
    "        (abs(aqi_df['Latitude'] - fire_lat) <= 0.5)\n",
    "        & (abs(aqi_df['Longitude'] - fire_lon) <= lon_degrees / 2))]\n",
    "    # If no measurements are within defined range, return an empty Series\n",
    "    if candidates.empty:\n",
    "        return pd.Series()\n",
    "    # Loop through the candidate measurements and get the nearest location\n",
    "    nearest = np.argmin([\n",
    "        get_distance(fire_lat, fire_lon, row.Latitude, row.Longitude)\n",
    "        for row in candidates.itertuples()\n",
    "    ])\n",
    "    return candidates.iloc[nearest]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "import sys\n",
    "\n",
    "# Define the number of rows, the interval at which to save partial results,\n",
    "# and create arrays for the data\n",
    "rows = df.shape[0]\n",
    "cutoff = int(rows / 10)\n",
    "percent = int(rows / 100)\n",
    "nearest_measurements = np.array(np.zeros(rows), dtype=object)\n",
    "aqi_readings = np.array(np.zeros(rows), dtype=np.float32)\n",
    "\n",
    "# Used to measure time elapsed in running the algorithm\n",
    "t1 = datetime.now()\n",
    "\n",
    "# Counter variable used for filenames for partial results\n",
    "counter = 1\n",
    "\n",
    "for row in df.itertuples():\n",
    "    # Write out progress\n",
    "    if row.Index % percent == 0:\n",
    "        sys.stdout.write(\"Progress: {}%  \\r\".format(\n",
    "            round((row.Index / rows) * 100, 2)))\n",
    "        sys.stdout.flush()\n",
    "    # Get the AQI measurement closest to the fire\n",
    "    aqi_candidates = get_aqis_during_fire(\n",
    "        str(row.FIRE_YEAR), row.DISCOVERY_DATE, row.CONT_DATE)\n",
    "    nearest = get_nearest_measurement(row.LATITUDE, row.LONGITUDE,\n",
    "                                      aqi_candidates)\n",
    "    nearest_measurements[row.Index] = nearest.to_dict()\n",
    "    if nearest.empty:\n",
    "        aqi_readings[row.Index] = np.nan\n",
    "    else:\n",
    "        aqi_readings[row.Index] = nearest['AQI']\n",
    "\n",
    "    # Save partial results in case script fails to run to the end\n",
    "    if row.Index > 0 and row.Index % cutoff == 0:\n",
    "        start = row.Index - cutoff\n",
    "        stop = row.Index\n",
    "        subframe = df.iloc[start:stop, :]\n",
    "        subframe['AQI'] = aqi_readings[start:stop]\n",
    "        subframe['AIR_QUALITY'] = nearest_measurements[start:stop]\n",
    "        subframe.to_csv(\n",
    "            'data/wildfire/wildfires_with_aqi_v2_{}.csv'.format(\n",
    "                str(counter).zfill(2)))\n",
    "        counter += 1\n",
    "\n",
    "# Add the air quality data to the wildfires dataframe, drop rows that do \n",
    "# not have an AQI reading (because no nearby measurements were found),\n",
    "# and save the results into a CSV file\n",
    "df['AQI'] = aqi_readings\n",
    "df['AIR_QUALITY'] = nearest_measurements\n",
    "df.dropna(subset=['AQI'], inplace=True)\n",
    "df.to_csv('data/wildfire/wildfires_with_aqi_all.csv')\n",
    "\n",
    "t2 = datetime.now()\n",
    "delta = t2 - t1\n",
    "\n",
    "print()\n",
    "print(\"Processing {} entries took {} minutes and {} seconds\".format(\n",
    "    rows, delta.seconds // 60, delta.seconds % 60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load in the new created data and delete measurement mistakes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "df=pd.read_csv('data/wildfire/wildfires_with_aqi_1.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "df=df[df['AQI']<500]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot of fire size and AQI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aqi=df['AQI'].as_matrix().tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fire=df['FIRE_SIZE'].as_matrix().tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bokeh\n",
    "from bokeh.plotting import figure\n",
    "from bokeh.models import ColumnDataSource\n",
    "from bokeh.plotting import figure\n",
    "from bokeh.io import show, output_notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = figure(plot_width = 600, plot_height = 600, \n",
    "           title = 'Example Glyphs',\n",
    "           x_axis_label = 'firesize', y_axis_label = 'aqi',\n",
    "          x_axis_type='log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "p.square(x=fire, y=aqi, size = 1, color = 'black')\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# linear regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import linear_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm = linear_model.LinearRegression()\n",
    "model = lm.fit(np.array(fire).reshape(-1,1),aqi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00020529618176290398"
      ]
     },
     "execution_count": 162,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lm.score(np.array(fire).reshape(-1,1),aqi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th>  <td>   0.002</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th>  <td>   0.002</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th>  <td>   509.7</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Sat, 27 Oct 2018</td> <th>  Prob (F-statistic):</th>  <td>9.17e-113</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>13:56:01</td>     <th>  Log-Likelihood:    </th> <td>-1.6156e+06</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>289728</td>      <th>  AIC:               </th>  <td>3.231e+06</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>289727</td>      <th>  BIC:               </th>  <td>3.231e+06</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>      <td> </td>     \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "   <td></td>     <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th> <td>    0.0010</td> <td> 4.57e-05</td> <td>   22.577</td> <td> 0.000</td> <td>    0.001</td> <td>    0.001</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>102089.827</td> <th>  Durbin-Watson:     </th>  <td>   0.300</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>   <td> 0.000</td>   <th>  Jarque-Bera (JB):  </th> <td>406148.491</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>            <td> 1.732</td>   <th>  Prob(JB):          </th>  <td>    0.00</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>        <td> 7.652</td>   <th>  Cond. No.          </th>  <td>    1.00</td> \n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.002\n",
       "Model:                            OLS   Adj. R-squared:                  0.002\n",
       "Method:                 Least Squares   F-statistic:                     509.7\n",
       "Date:                Sat, 27 Oct 2018   Prob (F-statistic):          9.17e-113\n",
       "Time:                        13:56:01   Log-Likelihood:            -1.6156e+06\n",
       "No. Observations:              289728   AIC:                         3.231e+06\n",
       "Df Residuals:                  289727   BIC:                         3.231e+06\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "x1             0.0010   4.57e-05     22.577      0.000       0.001       0.001\n",
       "==============================================================================\n",
       "Omnibus:                   102089.827   Durbin-Watson:                   0.300\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):           406148.491\n",
       "Skew:                           1.732   Prob(JB):                         0.00\n",
       "Kurtosis:                       7.652   Cond. No.                         1.00\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 163,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "\n",
    "X = fire\n",
    "#X = X[\"FIRE_SIZE\"]\n",
    "y = aqi\n",
    "\n",
    "# Note the difference in argument order\n",
    "model = sm.OLS(y, X).fit()\n",
    "predictions = model.predict(X) # make the predictions by the model\n",
    "\n",
    "# Print out the statistics\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['STAT_CAUSE_DESCR'].value_counts().plot(kind='barh',color='coral')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_arson = df[df['STAT_CAUSE_DESCR']=='Arson']\n",
    "df_arson['DAY_OF_WEEK'].value_counts().plot(kind='barh',color='coral')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_corr(df,size=10):\n",
    "    corr = df.corr()  #the default method is pearson\n",
    "    fig, ax = plt.subplots(figsize=(size, size))\n",
    "    ax.matshow(corr,cmap=plt.cm.Oranges)\n",
    "    plt.xticks(range(len(corr.columns)), corr.columns)\n",
    "    plt.yticks(range(len(corr.columns)), corr.columns)\n",
    "    for tick in ax.get_xticklabels():\n",
    "        tick.set_rotation(45)    \n",
    "    plt.show()        \n",
    "plot_corr(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfplot=df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfplot=dfplot.sort_index(by=['FIRE_SIZE'],axis=0, ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dflol=dfplot[:500]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dflol.plot(kind='scatter',x='LONGITUDE',y='LATITUDE',color='coral',alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_lightning = df#[df['STAT_CAUSE_DESCR']=='Lightning']\n",
    "df_lightning['DISCOVERY_DAY_OF_WEEK'].value_counts().plot(kind='barh',color='coral')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['DISCOVERY_DAY_OF_WEEK'].value_counts().plot(kind='barh',color='coral')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['STATE'].value_counts().head(n=10).plot(kind='barh',color='coral')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_CA = df[df['STATE']=='CA']\n",
    "df_CA['STAT_CAUSE_DESCR'].value_counts().plot(kind='barh',color='coral',title='causes of fires for CA')\n",
    "plt.show()"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}