{"cells":[{"metadata":{},"cell_type":"markdown","source":"## 0. Introduction\nIn this challenge, we used LightGBM (tree-based algorithm) to predict monthly emission factor for each power plant in Puerto Rico with a mean absolute percentage error of 21%. We think this is due to a lack of high quality data.\n\n### Response: emission factor\nEmission factor of a power plant is defined as the amount of NO2 around the plant per fuel combusted. \n\n`emission_factor` = `ColumnAmountNO2` / `fuel_combusted`\n\nUsing OMI data from NASA, we extracted the `ColumnAmountNO2` for each power plant and +/- 15 pixels from its surrounding (to be explained later). Since `ColumnAmountNO2` is recorded at different intervals throughout a month, we took the average reading for every month.\n\nWe also managed to collect monthly electricity generated `generated` for 22 power plants from [EIA](https://www.eia.gov/electricity/data/eia923/) in year 2018. Since data for 2019 is not available, we decided to use data from 2017 as a proxy. To calculate `fuel_combusted`, we multiplied `generated` with average the `heat_rate` for different energy source provided by [EIA](https://www.eia.gov/electricity/annual/html/epa_08_01.html). `heat_rate` is the amount of fuel combusted per Kilowatt of electricity generated.\n\n### Predictors: aggregate functions of band 1 to band 7 from Sentinel-5P OFFL NO2\nThe input to LightGBM is mainly data from Sentinel-5P OFFL NO2. Using `rasterstats`, we are able to generate the reading of every band based on the respective power plant. We then take the `mean`, `max`, `min`, and `std` of each plant based on the month."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport glob \nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom datetime import datetime\n\nimport rasterio as rio\nimport folium\nimport tifffile\n\n# Geospatial libraries that we will be using for this\nimport rasterio\nimport rasterstats\nimport folium\nimport geopandas as gpd\n\nfrom tqdm.notebook import trange, tqdm\n\nimport warnings\nwarnings.simplefilter(action='ignore')\n\npd.options.display.max_columns = None\n\nDATA_DIR = '../input/ds4g-environmental-insights-explorer/eie_data/'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 1. Dataset\nWe used an external [dataset](https://disc.gsfc.nasa.gov/datasets/OMNO2_003/summary?keywords=OMI) `omi` from the NASA. According to several research papers, such as [de Foy et. al](https://www.sciencedirect.com/science/article/abs/pii/S1352231015301291) and [Kharol et. al](https://www.sciencedirect.com/science/article/abs/pii/S1352231015302569), OMI (Ozone Monitoring Instrument) is used to capture the actual NO2 level. \n\nThe variable `ColumnAmountNO2` consists of [estimates](https://aura.gesdisc.eosdis.nasa.gov/data/Aura_OMI_Level2/OMNO2.003/doc/README.OMNO2.pdf) of NO2 in the total vertical column density and its uncertainty. We think this is a better measure of NO2 compared to `NO2_column_number_density` given in Sentinel-5P OFFL NO2 because of its popularity among research papers. \n\nWe also performed a statistical test to show that NO2 levels around power plants with coal, oil, and gas as primary fuels are higher than other areas. Please refer to 1. [here](https://www.kaggle.com/wjshenggggg/team-sunway-ds4g-appendix)."},{"metadata":{"trusted":true},"cell_type":"code","source":"omi = pd.read_csv(\"../input/omiv1/OMI-Aura_L2-OMNO2.csv\")\n\n# drop na\nomi = omi[~omi['ColumnAmountNO2'].isna()]\n\nomi.reset_index(inplace=True, drop=True)\n\nomi.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`power_plants` contains information on each plant and is provided."},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants = pd.read_csv(os.path.join(DATA_DIR, 'gppd', 'gppd_120_pr.csv'))\n\npower_plants.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Using `rasterstats`, we are able to generate the reading of every band from Sentinel-5P OFFL NO2 based on the respective power plant. `df_ras_gen` is the outcome after we take the `mean`, `max`, `min`, and `std` of each plant based on the month. Please refer to 2. [here](https://www.kaggle.com/wjshenggggg/team-sunway-ds4g-appendix) for more information on how `df_ras_gen` is constructed. Thanks to [Maxime](https://www.kaggle.com/maxlenormand/simplified-emission-for-each-plant-update-1) for sharing this approach.\n\nAdditionally, there is a `generated` column. It records the total amount of electricity generated by the respective plant in a given month. This [monthly data](https://www.eia.gov/electricity/data/eia923/) is made available by EIA. Unfortunately, we only have data for 22 out of 35 plants for the year 2018. We are less worried because the other 13 plants are all renewable energy plants. \n\nTo extrapolate the amount of electricity generated in 2019, we used the data in 2017 as a proxy since the electricity generated in the first half of 2018 was low due to a hurricane."},{"metadata":{"trusted":true},"cell_type":"code","source":"df_ras_gen = pd.read_csv('../input/df-ras-gen/df_ras_gen.csv')\n\ndf_ras_gen.head()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"COL = 'generated'\n\ng = df_ras_gen.groupby(['system:index'])[COL].mean()\n\nmerged_df = pd.merge(g, power_plants, on='system:index')\n\nnon_renew_energy = ['Oil', 'Coal', 'Gas']\n\nmerged_df['non_renew'] = 0\nmerged_df.loc[merged_df['primary_fuel'].isin(non_renew_energy), 'non_renew'] = 1\n\nkeep = ['system:index', COL, 'non_renew']\n\nno2_df = merged_df[keep]\n\nno2_df['system:index'] = no2_df['system:index'].str[-5:]\n\n# no2_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The plot below shows the average `generated` for each power plant. Red represents renewable energy while blue represents non-renewable energy power plants."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('Average `{}` for each power plant'.format(COL))\nbarlist = plt.bar(no2_df['system:index'], no2_df[COL])\n\nfor idx, row in no2_df.iterrows():\n    if row['non_renew'] == 0:\n        barlist[idx].set_color('r')\n        \nplt.xticks(rotation=90)        \nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"unique_plants = df_ras_gen['system:index'].unique()\n\n# remove the 13 plants that are not present in EIA data\npower_plants = power_plants[power_plants['system:index'].isin(unique_plants)]\n\npower_plants.reset_index(drop=True, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Extract `latitude` and `longtitude` from `power_plants`."},{"metadata":{"trusted":true},"cell_type":"code","source":"def split_column_into_new_columns(dataframe,\n                                  column_to_split,\n                                  new_column_one,\n                                  begin_column_one,\n                                  end_column_one):\n    '''Function provided by Paul,\n    https://www.kaggle.com/paultimothymooney/overview-of-the-eie-analytics-challenge\n    '''\n    for i in range(0, len(dataframe)):\n        dataframe.loc[i, new_column_one] = dataframe.loc[i, column_to_split][begin_column_one:end_column_one]\n    return dataframe\n\npower_plants = split_column_into_new_columns(power_plants,'.geo','latitude',50,66)\n\npower_plants = split_column_into_new_columns(power_plants,'.geo','longitude',31,48)\n\npower_plants['latitude'] = power_plants['latitude'].astype(float)\n\npower_plants['longitude'] = power_plants['longitude'].astype(float)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We scale the `latitude` and `longitude` to pixels of 148 and 475 to caputure the `ColumnAmountNO2` data around each power plant. Since OMI captures `ColumnAmountNO2` reading at different location, we decide to attribute the `ColumnAmountNO2` reading +/- 15 pixels from a power plant to itself. \n\nWe chose the `DEFINED_BOUNDS` to be 15 as it shows the highest level of NO2 in power plant areas compared to areas with no power plant."},{"metadata":{"trusted":true},"cell_type":"code","source":"def scale_lat_to_img_idx(lat):\n    lat = float(lat)\n    lat_img_idx = (18.6 - lat) * 148 / (18.6 - 17.9)\n    return int(lat_img_idx)\n\ndef scale_lon_to_img_idx(lon):\n    lon = float(lon)\n    lon_img_idx = (67.3 + lon) * 475 / (67.3 - 65.2)\n    return int(lon_img_idx)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants['img_idx_lt'] = power_plants['latitude'].apply(scale_lat_to_img_idx)\n\npower_plants['img_idx_lg'] = power_plants['longitude'].apply(scale_lon_to_img_idx)\n\nomi['img_idx_lt'] = omi['Latitude'].apply(scale_lat_to_img_idx)\n\nomi['img_idx_lg'] = omi['Longitude'].apply(scale_lon_to_img_idx)\n\n# filter out of bounds\nomi = omi[(omi['img_idx_lt'] >= 0) & (omi['img_idx_lt'] <= 148)]\nomi = omi[(omi['img_idx_lg'] >= 0) & (omi['img_idx_lt'] <= 475)]\n\nomi.reset_index(inplace=True, drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DEFINED_BOUNDS = 15\n\ndic = {}\n\nfor idx, row in power_plants.iterrows():\n    dic[row['system:index']] = []\n    \n    lat, lon = row['img_idx_lt'], row['img_idx_lg']\n    \n    for _, omi_row in omi.iterrows():\n        \n        omi_lat, omi_lon = omi_row['img_idx_lt'], omi_row['img_idx_lg']\n        \n        lat_in_bounds = omi_lat - DEFINED_BOUNDS <= lat <= omi_lat + DEFINED_BOUNDS\n        lon_in_bounds = omi_lon - DEFINED_BOUNDS <= lon <= omi_lon + DEFINED_BOUNDS\n        \n        # if the OMI reading of a given location is\n        # +/- 15 pixels from a power plant,\n        # assign the reading to the power plant.\n        if lat_in_bounds and lon_in_bounds:\n            d = {}\n            d['time'] = omi_row['Time']\n            d['ColumnAmountNO2'] = omi_row['ColumnAmountNO2']\n            \n            dic[row['system:index']].append(d)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We then take the mean of `ColumnAmountNO2` based on power plant and month. Unfortunately, some plants do not have any `ColumnAmountNO2` recorded for the given month."},{"metadata":{"trusted":true},"cell_type":"code","source":"final = pd.DataFrame()\n\n# construct a df for all the plants\n# with the recorded ColumnAmountNO2 and time\nfor k in dic.keys():\n    tmp = pd.DataFrame()\n    tmp['system:index'] = [k] * len(dic[k])\n    tmp = pd.concat([tmp, pd.DataFrame(dic[k])], axis=1)\n    final = pd.concat([final, tmp], axis=0, ignore_index=True)\n\nfinal['month'] = pd.to_datetime(final['time']).dt.month\n\ng = final.groupby(['system:index', 'month'])['ColumnAmountNO2'].mean()\n\ng = g.reset_index()\n\ng.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Merge `ColumnAmountNO2` with `df_ras_gen`! We remove any entry that has `generated` of 0. We suspect this is due to insufficient information."},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.merge(df_ras_gen, g, on=['system:index', 'month'])\n\ndf = df.drop_duplicates()\n\ndf['generated'] = df['generated'].astype(np.float32)\n\n# remove `generated` of 0\ndf = df[df['generated'] > 0.01]\n\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"monthly_elec = df.groupby('month')['generated'].sum()\n\nplt.title('Total electricity generated')\nplt.xlabel('Months')\nplt.plot(monthly_elec)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We believe the dip is due to insufficient information by EIA."},{"metadata":{},"cell_type":"markdown","source":"## 2. Compute emission factor (EF)\nSo far, we have the amount of electricity generated by each power plant for each month. But as outlined in the challenge, we need the estimate the emission factor.\n\nAccording to a [post](https://www.enveraconsulting.com/how-to-calculate-air-emissions/) by Envera Consulting, we can think of emission factor as the pounds of NOx per million cubic feet of fuel combusted. \n\nUsing [data](https://www.eia.gov/electricity/annual/html/epa_08_01.html) provided by EIA, we can estimate fuel combusted based on the source. For instance, for 1 Kilowatthour of electricity produced, 10,481 Btu of coal has to be combusted. The conversion rate table is as follows. "},{"metadata":{"trusted":true},"cell_type":"code","source":"heat_rate_dict = {\n    'Coal': 10_481,\n    'Oil': 11_095,\n    'Gas': 7_821\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The fact that [EIA considers](https://www.eia.gov/tools/faqs/faq.php?id=74&t=11) electricity generation from biomass, hydro, solar, and wind to be carbon neutral, we consider those plants as having 0 emission factor. Hence, they are removed from our dataframe."},{"metadata":{"trusted":true},"cell_type":"code","source":"plant_fuel = power_plants[['system:index', 'primary_fuel']]\n\nplant_fuel_dict = {}\n\nfor _, row in plant_fuel.iterrows():\n    plant_fuel_dict[row['system:index']] = row['primary_fuel'] \n\ndf['primary_fuel'] = df['system:index'].map(plant_fuel_dict)\n\ndf = df[df['primary_fuel'].isin(non_renew_energy)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Assign the `heat_rate` to power plants. As defined above, emission factor is now equal to `ColumnAmountNO2` / `fuel_combusted`\n\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"df['heat_rate'] = df['primary_fuel'].map(heat_rate_dict)\n\ndf['fuel_combusted'] = df['heat_rate'] * df['generated']\n\ndf['emission_factor'] = df['ColumnAmountNO2'] / (df['fuel_combusted'])\n\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We visualize the distribution of emission factors and the emission factors based on each month below."},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(df['emission_factor'])\nplt.title('Distribution of emission factors')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We observe that the emission factor is right skew. We decide to clip the emission factor greater than 90th quantile and lesser than 10th quantile."},{"metadata":{"trusted":true},"cell_type":"code","source":"top_quantile = df['emission_factor'].quantile(.9)\nbottom_quantile = df['emission_factor'].quantile(.1)\n\ndf.loc[df['emission_factor'] > top_quantile, 'emission_factor'] = top_quantile\ndf.loc[df['emission_factor'] < bottom_quantile, 'emission_factor'] = bottom_quantile","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Based on our methodology, the monthly emission factor is as follows."},{"metadata":{"trusted":true},"cell_type":"code","source":"monthly_ef = df.groupby('month')['emission_factor'].sum()\n\nplt.title('Total emission factors')\nplt.xlabel('Months')\nplt.plot(monthly_ef)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Model\nWe will use `lightgbm` to predict emission factors using data from Sentinel-5P OFFL NO2 after taking the monthly aggregate functions for each plant."},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgb\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"exclude = ['generated', 'primary_fuel', 'ColumnAmountNO2', 'heat_rate',\n          'year', 'month',\n           'fuel_combusted', 'emission_factor']\n\nle = LabelEncoder()\ndf['system:index'] = le.fit_transform(df['system:index'])\n\ncols = [c for c in df.columns if c not in exclude]\n\nX_tr = df[cols]\ny_tr = df['emission_factor']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since emission factors vary a lot, we decided to use mean absolute percentage error (MAPE) to gauge the peformance of our model. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Setting hyperparameters\nlgb_params = {\n    'learning_rate': 0.05,\n    'metric': 'mape',\n    'objective': 'regression',\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create train and validation set\nX_train, X_test, y_train, y_test = train_test_split(X_tr, y_tr, test_size=0.2, random_state=42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"trn_data = lgb.Dataset(X_train.values, label=y_train)\nval_data = lgb.Dataset(X_test.values, label=y_test)\n\nlgb_clf = lgb.train(lgb_params,\n                    trn_data,\n                    1000,\n                    valid_sets=[trn_data, val_data],\n                    early_stopping_rounds=200,\n                    verbose_eval=200)\n\ny_pred = lgb_clf.predict(X_test, \n                         num_iteration=lgb_clf.best_iteration)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('Prediction vs Actual')\nplt.plot(y_pred, 'r--', label='predict')\nplt.plot(y_test.values, color='b', label='actual')\nplt.legend(loc='upper right')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Our model has a MAPE of 21%, likely due to insufficient data. We expect our model to perform better with more data provided."},{"metadata":{},"cell_type":"markdown","source":"## 4. Analysis: feature importance\nAs seen below, other than `system:index`, statistics of `troposheric_NO2_column_density`, `cloud_fraction`, and `NO2_column_number_density` play important roles in predicting the `emission_factor`."},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_score = lgb_clf.feature_importance()\n\nfeatures_df = pd.DataFrame({'features': cols,\n                            'importance': feature_score})\n\nfeatures_df = features_df.sort_values('importance', ascending=False)\n\nplt.figure(figsize=(8,6))\nsns.barplot(x=\"importance\", y=\"features\", \n            data=features_df)\nplt.title('LightGBM Features')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 5. Recommendation\n### Scalability and accuracy\nThis `lightgbm` is certainly scalable. It does not fare well on accuracy due to limited data.\n\n### Pros and cons of the model\n* Pros: Cost-effective as it only requires satelite data. \n* Cons: The model assumes `ColumnAmountNO2` of +/- 15 pixels from a power plant belongs to its \"emission\"; it only accounts for power plants with primary fuel of oil, gas, and coal and assumes renewable energy plant does not contribute to emission factor; construction of the response variable `emission_factor` relies heavily on accurate data on electricity generated by a power plant and heat rate.\n\n## 6. Conclusion\nThanks to Kaggle for organizing such an interesting competition! Although we did not have sufficient high quality data, we relished in using our creativity to explore various methods to model and predict the emission factor."}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}