Environmental_Features_Creation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <HR>Creation of Environmental Features for Wildfires Incidents Using NCEP data<HR>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import the Required Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import json\n",
    "import glob\n",
    "import time\n",
    "import pickle\n",
    "import warnings\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import concurrent.futures\n",
    "import dask.dataframe as dd\n",
    "from tqdm.notebook import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_PATH = \"data/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 'EnvironmentalFeaturesGenerator' class to extract the environmental features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class EnvironmentalFeaturesGenerator:\n",
    "    \"\"\"FeaturesGenerator to generate environmental features.\n",
    "    \n",
    "    It generates the features like Air_Temperature('air'), Geopotential_Height('hgt'), \n",
    "    Relative_Humidity('rhum'), East_Wind/Zonal_Wind('uwnd'), North_Wind/Meridional_Wind('vwnd'),  \n",
    "    Tropopause_Pressure('trpp'), Tropopause_Temperature('trpt'), Surface_Potential_Temperature('srfpt') at different pressure level using NCEP data. \n",
    "    \n",
    "    NCEP data link:\n",
    "    https://psl.noaa.gov/data/gridded/data.ncep.html\n",
    "    https://www.ncei.noaa.gov/erddap/convert/oceanicAtmosphericVariableNames.htmlTable\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, wildfires_data, ncep_data, ncep_sfc_data):\n",
    "        self.wildfires_data = wildfires_data\n",
    "        self.ncep_data = ncep_data\n",
    "        self.levels = [100, 150, 200, 500, 1000]\n",
    "        self.ncep_data_vars = set(ncep_data[list(ncep_data)[0]]) - set(['head'])\n",
    "        self.ncep_sfc_data = ncep_sfc_data\n",
    "        self.ncep_sfc_data_vars = set(ncep_sfc_data[list(ncep_sfc_data)[0]]) - set(['head'])\n",
    "        \n",
    "    def _extract_features(self, row):\n",
    "        \"\"\"extract the environmental features for particular wildfire incident\"\"\"\n",
    "        ncep_data = self.ncep_data\n",
    "        ncep_sfc_data = self.ncep_sfc_data\n",
    "        date = row['date']\n",
    "        features = dict(row)\n",
    "        #reduce the dimensions of ncep_data(xarray dataset) by fixing coordinates(lon,lat)\n",
    "        #and then convert it to dataframe\n",
    "        ncep_data = ncep_data[date.year] \\\n",
    "                            .sel(lon=row['longitude'], lat=row['latitude'], method='nearest') \\\n",
    "                            .to_dask_dataframe() \\\n",
    "                            .compute() \\\n",
    "                            .set_index(['level','time'])\n",
    "        #reduce the dimensions of ncep_sfc_data(xarray dataset) by fixing coordinates(lon,lat)\n",
    "        #and then convert it to dataframe\n",
    "        ncep_sfc_data = ncep_sfc_data[date.year] \\\n",
    "                            .sel(lon=row['longitude'], lat=row['latitude'], method='nearest') \\\n",
    "                            .to_dask_dataframe() \\\n",
    "                            .compute() \\\n",
    "                            .set_index(['time'])\n",
    "\n",
    "        for level in self.levels:\n",
    "            #features at different pressure level\n",
    "            point = ncep_data.loc[level]\n",
    "            p1w = point.rolling(7).mean()  # 1 Week mean\n",
    "            p2w = point.rolling(14).mean() # 2 Week mean\n",
    "            p3w = point.rolling(21).mean() # 3 Week mean\n",
    "            # \n",
    "            v0w = point.loc[date]\n",
    "            v1w = p1w.loc[date]\n",
    "            v2w = p2w.loc[date]\n",
    "            v3w = p3w.loc[date]\n",
    "            #\n",
    "            for data_var in self.ncep_data_vars:\n",
    "                features[\"{0}_0w_lvl_{1}\".format(data_var,level)] = v0w[data_var]\n",
    "                features[\"{0}_1w_lvl_{1}\".format(data_var,level)] = v1w[data_var]\n",
    "                features[\"{0}_2w_lvl_{1}\".format(data_var,level)] = v2w[data_var]\n",
    "                features[\"{0}_3w_lvl_{1}\".format(data_var,level)] = v3w[data_var]\n",
    "        #features at surface level\n",
    "        point = ncep_sfc_data\n",
    "        p1w = point.rolling(7).mean()  # 1 Week mean\n",
    "        p2w = point.rolling(14).mean() # 2 Week mean\n",
    "        p3w = point.rolling(21).mean() # 3 Week mean\n",
    "        # \n",
    "        v0w = point.loc[date]\n",
    "        v1w = p1w.loc[date]\n",
    "        v2w = p2w.loc[date]\n",
    "        v3w = p3w.loc[date]\n",
    "        #\n",
    "        for data_var in self.ncep_sfc_data_vars:\n",
    "            features[\"{0}_0w\".format(data_var)] = v0w[data_var]\n",
    "            features[\"{0}_1w\".format(data_var)] = v1w[data_var]\n",
    "            features[\"{0}_2w\".format(data_var)] = v2w[data_var]\n",
    "            features[\"{0}_3w\".format(data_var)] = v3w[data_var] \n",
    "\n",
    "        return features\n",
    "    \n",
    "    def _get_day_info(self, df_features):\n",
    "        \"\"\" \"\"\"\n",
    "        df_features['day'] = df_features.date.dt.day\n",
    "        df_features['month'] = df_features.date.dt.month\n",
    "        df_features['year']= df_features.date.dt.year\n",
    "        df_features['week_day'] = df_features.date.dt.weekday\n",
    "        df_features['weekofyear'] = df_features.date.dt.weekofyear\n",
    "        df_features['is_winter'] = df_features.apply(lambda x: int(x['month'] in [12,1,2]), axis=1)\n",
    "        df_features['is_autumn'] = df_features.apply(lambda x: int(x['month'] in [9,10,11]), axis=1)\n",
    "        df_features['is_summer'] = df_features.apply(lambda x: int(x['month'] in [6,7,8]), axis=1)\n",
    "        df_features['is_spring'] = df_features.apply(lambda x: int(x['month'] in [3,4,5]), axis=1)\n",
    "        return df_features\n",
    "        \n",
    "    def _thread_func(self, year):\n",
    "        wildfires_data = self.wildfires_data\n",
    "        file_path = DATA_PATH + 'features/' + str(year) + '.pickle'\n",
    "        df_features = {}\n",
    "        try:\n",
    "            if os.path.isfile(file_path):\n",
    "                with open(file_path, 'rb') as f:\n",
    "                    df_features = pickle.load(f)\n",
    "                    print('Loaded: ',len(df_features)) \n",
    "            start_year, end_year = year, year+1\n",
    "            df_subsample = wildfires_data.query('(date >= @start_year) & (date < @end_year)')\n",
    "            for row_id, row in tqdm(df_subsample.iterrows(), total=df_subsample.shape[0], desc=\"Year \"+str(year)):\n",
    "                if row_id not in df_features:\n",
    "                    df_features[row_id] = self._extract_features(row)\n",
    "                    with open(file_path, 'wb') as f:\n",
    "                        pickle.dump(df_features, f)\n",
    "            df_features = self._get_day_info(pd.DataFrame(df_features.values()))           \n",
    "            df_features.set_index('fire_id').to_csv(DATA_PATH + 'features/' + str(year) + '.csv')\n",
    "        except Exception as e:\n",
    "            print(e)\n",
    "            \n",
    "    def start(self, years):\n",
    "        \"\"\"start extracting the features using multiple threads, one thread per year \"\"\"\n",
    "        print(\"=\"*5,\"Features Extraction has started\",\"=\"*5)\n",
    "        time.sleep(0.5)\n",
    "        start = time.time()\n",
    "        with concurrent.futures.ThreadPoolExecutor(max_workers=8) as executor:\n",
    "            for year in years:\n",
    "                executor.submit(self._thread_func, year)\n",
    "        end = time.time()\n",
    "        print(\"\\n\",\"=\"*5,\"Features Extraction has finished\",\"=\"*5)\n",
    "        print(\"Total time taken\", end-start, \" sec\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 'load_Dataset' method to load the data into the memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_Dataset(years):\n",
    "    \"\"\"load the wildfire data and NCEP data into the memory\"\"\"\n",
    "    ncep_data, ncep_sfc_data = dict(), dict()\n",
    "    #Mutli-dimensional ncep data at different pressure level and surface level\n",
    "    for year in years:\n",
    "        ncep_data[year] = xr.open_mfdataset('data/ncep/*.'+str(year)+'.nc', combine='by_coords', parallel=True)\n",
    "        ncep_sfc_data[year] = xr.open_mfdataset('data/ncep/ncep_sfc/*.'+str(year)+'.nc', combine='by_coords', parallel=True)\n",
    "    #dataset for wildfire incidents\n",
    "    wildfires_train = pd.read_csv('data/wildfires_train.csv', parse_dates=['date'])\n",
    "    wildfires_train['fire_type_name_en'] = \\\n",
    "        wildfires_train.fire_type_name.map(\n",
    "        json.load(open('data/fire_type_name_en.json','r',encoding='utf-8')),\n",
    "        na_action='ignore')\n",
    "    return wildfires_train, ncep_data, ncep_sfc_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extract Environmental Features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Features Extraction has started =====\n",
      "Loaded:  9593\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "56d2c2489f1f47cf875bf28128c4f67a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Year 2012', max=9593.0, style=ProgressStyle(description_w…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Loaded:  16418\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "643fa5f0dbc549f68930b6bf0a995a94",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Year 2013', max=16418.0, style=ProgressStyle(description_…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Loaded:  28179\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "af9fc19104f94082b706cb0e1e3e4e4a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Year 2015', max=28179.0, style=ProgressStyle(description_…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Loaded:  35530\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "37c2a5a02a1a48978c9919cce297e41b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Year 2014', max=35530.0, style=ProgressStyle(description_…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      " ===== Features Extraction has finished =====\n",
      "Total time taken 98.7444679737091  sec\n"
     ]
    }
   ],
   "source": [
    "years = range(2012,2016)\n",
    "generator = EnvironmentalFeaturesGenerator(*load_Dataset(years))\n",
    "generator.start(years)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Features Extraction has started =====\n",
      "Loaded:  16011\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "57caa5912998447ea0c48fbc9d053f85",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Year 2019', max=16011.0, style=ProgressStyle(description_…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Loaded:  23070\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6d08164b4605460f8f28a7c2b1248696",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Year 2016', max=23070.0, style=ProgressStyle(description_…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded:  22820\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "adc0437df0f349b6884669af0fe6f020",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Year 2018', max=22820.0, style=ProgressStyle(description_…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded:  23250\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "21b99781473240748192e217c2956f4e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='Year 2017', max=23250.0, style=ProgressStyle(description_…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "\n",
      "\n",
      " ===== Features Extraction has finished =====\n",
      "Total time taken 91.26349234580994  sec\n"
     ]
    }
   ],
   "source": [
    "years = range(2016,2020)\n",
    "generator = EnvironmentalFeaturesGenerator(*load_Dataset(years))\n",
    "generator.start(years)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Merge all years features data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pickle_env_features = {}\n",
    "# for file in glob.glob(DATA_PATH + 'features/'+'*.pickle'):\n",
    "#     with open(file, 'rb') as f:\n",
    "#         pickle_env_features.update(pickle.load(f))  \n",
    "# with open(DATA_PATH + 'env_features.pickle', 'wb') as f:\n",
    "#     pickle.dump(pickle_env_featuress, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_env_features = dd.read_csv(DATA_PATH + 'features/*.csv', parse_dates=['date']).compute()\n",
    "df_env_features.sort_values(\"fire_id\", inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(174871, 188)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_env_features.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Missing Values in features data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": true
   },
   "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>Count</th>\n",
       "      <th>Pert</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>srfpt_3w_lvl_150</th>\n",
       "      <td>59391</td>\n",
       "      <td>33.962750</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>srfpt_3w_lvl_200</th>\n",
       "      <td>59391</td>\n",
       "      <td>33.962750</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>srfpt_3w_lvl_100</th>\n",
       "      <td>59391</td>\n",
       "      <td>33.962750</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>srfpt_2w_lvl_200</th>\n",
       "      <td>59219</td>\n",
       "      <td>33.864391</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>srfpt_2w_lvl_100</th>\n",
       "      <td>59219</td>\n",
       "      <td>33.864391</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>srfpt_1w_lvl_1000</th>\n",
       "      <td>231</td>\n",
       "      <td>0.132097</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rhum_1w</th>\n",
       "      <td>231</td>\n",
       "      <td>0.132097</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rhum_1w_lvl_1000</th>\n",
       "      <td>231</td>\n",
       "      <td>0.132097</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>air_1w_lvl_200</th>\n",
       "      <td>231</td>\n",
       "      <td>0.132097</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>air_1w_lvl_150</th>\n",
       "      <td>231</td>\n",
       "      <td>0.132097</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>144 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                   Count       Pert\n",
       "srfpt_3w_lvl_150   59391  33.962750\n",
       "srfpt_3w_lvl_200   59391  33.962750\n",
       "srfpt_3w_lvl_100   59391  33.962750\n",
       "srfpt_2w_lvl_200   59219  33.864391\n",
       "srfpt_2w_lvl_100   59219  33.864391\n",
       "...                  ...        ...\n",
       "srfpt_1w_lvl_1000    231   0.132097\n",
       "rhum_1w              231   0.132097\n",
       "rhum_1w_lvl_1000     231   0.132097\n",
       "air_1w_lvl_200       231   0.132097\n",
       "air_1w_lvl_150       231   0.132097\n",
       "\n",
       "[144 rows x 2 columns]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "missing_value = pd.DataFrame(df_env_features.isnull().sum(axis=0), columns=['Count'])\n",
    "missing_value['Pert'] = missing_value.Count*100/df_env_features.shape[0]\n",
    "missing_value.sort_values(by='Count', ascending=False, inplace=True)\n",
    "missing_value[missing_value.Count>0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fill the missing values\n",
    "df_env_features.fillna(0, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save environmental features into csv file\n",
    "df_env_features.to_csv(DATA_PATH + 'env_features.csv', index=False)"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}