{"cells":[{"metadata":{"id":"x0b46VULTw3T","colab_type":"text"},"cell_type":"markdown","source":"# DS4G : Environmental Insights Explorer\n# Exploring Alternatives for Emission Factor Calculations\n# (Domain - Environment) (Analytics Prediction)"},{"metadata":{"id":"rbePaaeoZVuf","colab_type":"text"},"cell_type":"markdown","source":"## The entire solution is divided into following notebooks\n\n####1.   [DS4G Environmental Insights Explorer- Introduction.(Part 1)](https://www.kaggle.com/meenakshiramaswamy/ds4g-eie-no2emission-earthenginereading)\n\n\n####2.   [Exploration on data points and detailed analysis.(Part 2)](https://www.kaggle.com/meenakshiramaswamy/ds4g-eie-no2emission-eda)\n\n####3.   [Emission Factor Calculation and Predictive model building.(Part 3)](https://www.kaggle.com/meenakshiramaswamy/ds4g-eie-calcef-buildmodel/) \n\n"},{"metadata":{"id":"_HcMVBq8bKtg","colab_type":"text"},"cell_type":"markdown","source":"#### [Part 3 -  Model Building and Emission Factor Calculation] subdivided into following sections\n\n\n*   Emission factor Calculation \n*   Model Building\n*   Time Series - SARIMAX\n*   Time Series - Prophat\n*   KMeans Clustering\n*   Supervised Algorithm\n*   End Note\n*   Acknowledgement\n\n\n"},{"metadata":{"trusted":true,"id":"EUzPzb76gNEy","colab_type":"code","colab":{}},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport math as m\nfrom pandas import read_csv\n\nfrom datetime import datetime\n\nimport os\nimport glob\n\nimport rasterio as rio\nimport folium \n\nimport geopandas\nimport tifffile as tiff\nfrom folium import plugins\nfrom shapely.geometry import Point\nimport rasterstats\nfrom rasterstats import zonal_stats, point_query","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"id":"HJU3uTye5SdE","colab_type":"code","colab":{},"_kg_hide-input":true},"cell_type":"code","source":"pd.set_option('max_columns', 500)\npd.set_option('max_rows', 500)\nimport warnings\nwarnings.filterwarnings('ignore')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Code for displaying plotly express plot\ndef configure_plotly_browser_state():\n  import IPython\n  display(IPython.core.display.HTML('''\n        <script src=\"/static/components/requirejs/require.js\"></script>\n        <script>\n          requirejs.config({\n            paths: {\n              base: '/static/base',\n              plotly: 'https://cdn.plot.ly/plotly-latest.min.js?noext',\n            },\n          });\n        </script>\n        '''))\n  ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from plotly import tools, subplots\nimport plotly.offline as py\npy.init_notebook_mode(connected = True)\nimport plotly.graph_objs as go\nimport plotly.express as px\nconfigure_plotly_browser_state()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"id":"LIDs4Md05SdL","colab_type":"code","colab":{},"_kg_hide-input":true},"cell_type":"code","source":"from IPython.display import display\nfrom IPython.core.interactiveshell import InteractiveShell\nInteractiveShell.ast_node_interactivity = \"all\"","execution_count":null,"outputs":[]},{"metadata":{"id":"O19u1TkN5Sef","colab_type":"text"},"cell_type":"markdown","source":"## Read all the data"},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path = \"/kaggle/input/ds4g-environmental-insights-explorer\"\n\npower_plants_df = pd.read_csv(\"/kaggle/input/ds4g-eie-no2emission-eda/power_plants.csv\")\npower_plants_fossil_df = pd.read_csv(\"/kaggle/input/ds4g-eie-no2emission-eda/fossil.csv\")\n\nno2_dfs_gldas_df = pd.read_csv(\"/kaggle/input/ds4g-eie-no2emission-eda/no2_dfs_gldas_df.csv\")\npower_plants_renew_df = pd.read_csv(\"/kaggle/input/ds4g-eie-no2emission-eda/renew.csv\")\n#pp_gdf = pd.read_csv(\"/kaggle/input/kaggle-ds4g-part2/pp_gdf.csv\")","execution_count":null,"outputs":[]},{"metadata":{"id":"UPW_AQCUWo4F","colab_type":"text"},"cell_type":"markdown","source":"### Using Geopandas we convert the coordinates into geometry and plot them in the map"},{"metadata":{"id":"5Ki1OdFMWm49","colab_type":"code","outputId":"61fb9c70-7d29-4375-9e58-5fb57456ce21","executionInfo":{"status":"ok","timestamp":1584718561524,"user_tz":-330,"elapsed":2708,"user":{"displayName":"Meenakshi Ramaswamy","photoUrl":"https://lh4.googleusercontent.com/-3sXdlOXrEr8/AAAAAAAAAAI/AAAAAAAACC8/ngMx9EiU3T4/s64/photo.jpg","userId":"09422112087559126593"}},"colab":{"base_uri":"https://localhost:8080/","height":225},"trusted":true},"cell_type":"code","source":"pp_gdf = geopandas.GeoDataFrame(\n    power_plants_df, geometry=geopandas.points_from_xy(power_plants_df.longitude, power_plants_df.latitude))\n\nworld = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))\n\n# We restrict to Puerto Rica.\nax = world[world.continent == 'Puerto Rica'].plot(\n    color='white', edgecolor='black')\n\n# We can now plot our ``GeoDataFrame``.\npp_gdf.plot(ax=ax, color='magenta')\n\nplt.show()\npp_gdf['name'] = pp_gdf['name'] + ':' + pp_gdf['primary_fuel']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"id":"r_ulvW4w5Sgl","colab_type":"code","cellView":"form","colab":{},"_kg_hide-input":true},"cell_type":"code","source":"#@title\n'''\nimport ee\nfrom kaggle_secrets import UserSecretsClient\nfrom google.oauth2.credentials import Credentials\nuser_secrets = UserSecretsClient()\nsecret_value_0 = user_secrets.get_secret(\"akshi\")\n'''","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true,"id":"2zz6uHVy5Sgp","colab_type":"code","cellView":"form","colab":{},"_kg_hide-input":true},"cell_type":"code","source":"#@title\n'''\nee.Authenticate()\n!cat ~/.config/earthengine/credentials\n'''","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"id":"iSOKsm7q5Sgr","colab_type":"code","cellView":"form","colab":{},"_kg_hide-input":true},"cell_type":"code","source":"#@title\n'''\n#user_secret = \"\" # Your user secret, defined in the add-on menu of the notebook editor\n#refresh_token = UserSecretsClient().get_secret(user_secret)\ncredentials = Credentials(\n        None,\n        refresh_token=secret_value_0,\n        token_uri=ee.oauth.TOKEN_URI,\n        client_id=ee.oauth.CLIENT_ID,\n        client_secret=ee.oauth.CLIENT_SECRET,\n        scopes=ee.oauth.SCOPES)\nee.Initialize()\n#ee.Initialize(credentials=credentials)\n'''","execution_count":null,"outputs":[]},{"metadata":{"id":"5t1jOsVTJtdW","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"def initMap(df, lat, lon):\n    location = [lat, lon]\n    Map = folium.Map(location=location, zoom_start=8)\n    return Map\n","execution_count":null,"outputs":[]},{"metadata":{"id":"8QDjcTedhtIP","colab_type":"text"},"cell_type":"markdown","source":"### from https://www.kaggle.com/maxlenormand/simplified-emission-for-each-plant-wip"},{"metadata":{"id":"13Qb3VrMhFPf","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"pp_geo_df_ = pp_gdf.copy()\npp_geo_df_['geometry'] = pp_gdf.geometry.buffer(0.05)","execution_count":null,"outputs":[]},{"metadata":{"id":"BfR9icPUhpgV","colab_type":"code","outputId":"83b451a2-7efc-4bc9-e8de-7ba6178d16da","executionInfo":{"status":"ok","timestamp":1584720693376,"user_tz":-330,"elapsed":2010015,"user":{"displayName":"Meenakshi Ramaswamy","photoUrl":"https://lh4.googleusercontent.com/-3sXdlOXrEr8/AAAAAAAAAAI/AAAAAAAACC8/ngMx9EiU3T4/s64/photo.jpg","userId":"09422112087559126593"}},"colab":{"base_uri":"https://localhost:8080/","height":1000},"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"lat=18.140178; lon=-66.664513 #puerto rico\nMap = initMap(power_plants_df, lat, lon)   \n#plot_points_on_map(power_plants,lat,lon,10)\nfor i in range(pp_geo_df_.shape[0]):\n    folium.GeoJson(pp_geo_df_.geometry[i]).add_to(Map)\n\nMap","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### What is Emissions Measurement ?\nEmissions measurement is the process of measuring the amount of pollutants, in a gaseous or particulate form, being emitted to the air from a specific source, such as an industrial process.\nMeasurements of emissions can be used to understand the relative importance of a given source compared to other sources and in developing emissions inventories.\n### What is Emission Factor ?\nRepresentative value that attempts to relate the quantity of a pollutant released to the atmosphere with an activity associated with the release of that pollutant.\nUsually expressed as the weight of pollutant divided by a unit weight, volume, distance, or duration of the activity emitting the pollutant (e.g., lbs of NOx emitted per gallon fuel combusted)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def img_stat(pp_geometry, name, primary_fuel, j):\n  global res , pp_monthwise\n\n  stats = []\n  dates = []\n  #s5p_no2_20180701T161259_20180707T175356\n  s5p_no2_path = glob.glob(data_path+'/eie_data/s5p_no2/*')\n  for (i, image) in enumerate(s5p_no2_path):\n      date = datetime.strptime(image[76:84], '%Y%m%d')\n      stat = rasterstats.zonal_stats(pp_geometry,\n                                    image,\n                                    band=2, #2: tropospheric_NO2_column_number_density\n                                    stats=['mean'])\n      stat = stat[0] # get location of pp\n      stat = stat['mean'] # retrieve stat\n      dates.append(date)\n      stats.append(stat)\n      \n  results = pd.DataFrame(index=dates, data=stats, columns=[primary_fuel])\n  results.plot()\n  plt.title('tropospheric NO2_column_number_density over time period\\n')   \n  # chking -------\n  pp_all_month = pd.DataFrame(list(zip(dates, stats)), columns= ['Date','no2-mean'], index=None)\n  pp_all_month['pp'] = name\n  if j==0 :\n    res = pp_all_month\n  else:\n    res = res.append(pp_all_month, ignore_index= True)  \n     \n  # perform GroupBy operation over monthly frequency\n  res.Date = pd.to_datetime(res.Date,format='%Y-%m-%d')\n  a = res.set_index('Date').groupby(pd.Grouper(freq='M'))['no2-mean'].sum().reset_index()\n  a['pp'] = name\n  if j==0 :\n    pp_monthwise = a\n  else:\n    pp_monthwise = pp_monthwise.append(a, ignore_index = True )\n\n  # chking -------\n  x = results.sum()\n  #print (x)\n  new_row = {\"plant_name\": name, \"no2_emission\": x}\n  emission_df.loc[j] = new_row\n\n  return pp_monthwise   ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"emission_df = pd.DataFrame()\n\n#pp_all_month = pd.DataFrame(columns= ['Date','no2-mean','pp'])\n\nemission_df['plant_name'] =\"\"\nemission_df['no2_emission'] = 0.0\nfor i in range(0,pp_gdf.shape[0]):\n  img_stat(pp_gdf.geometry[i], pp_gdf.name[i], pp_gdf.primary_fuel[i], i )\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"emission_df.info()\nemission_df['plant_name'] = emission_df['plant_name'].astype(str)\nemission_df['no2_emission'] = emission_df['no2_emission'].astype(float)\n\npp_monthwise.info()\npp_monthwise['Date'] = pd.to_datetime(pp_monthwise['Date'])\npp_monthwise['no2-mean'] = pp_monthwise['no2-mean'].astype(float)\npp_monthwise['pp'] = pp_monthwise['pp'].astype(str)\n\npp_monthwise.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"emission_df.to_csv('emission_df.csv')\npp_monthwise.to_csv('pp_monthwise.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"'''\nplt.figure(figsize=(40,20))  \nfig, ax = plt.subplots()    \nwidth = 0.5 # the width of the bars \nind = np.arange(len(emission_df['no2_emission']))  # the x locations for the groups\nax.barh(ind, emission_df['no2_emission'], width, color=\"blue\")\nax.set_yticks(ind+width/2)\nax.set_yticklabels(emission_df['plant_name'], minor=False)\nplt.title('NO2 Emission')\nplt.xlabel('no2 emission')\nplt.ylabel('plant name')      \n#plt.show()\nplt.savefig(os.path.join('test.png'), dpi=300, format='png', bbox_inches='tight') # use format='svg' or 'pdf' for vectorial pictur\n'''","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from matplotlib import style\nstyle.use('ggplot')\nplt.figure(figsize=(12,8))                                 # Setting the figure size\nsns.distplot(emission_df['no2_emission'])                     # Creating the histogram\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### All PR Power plants NO2 annual emission "},{"metadata":{"trusted":true},"cell_type":"code","source":"import plotly.graph_objects as go\nconfigure_plotly_browser_state()\nfig = go.Figure(data=[go.Bar(\n            x=emission_df['plant_name'], y=emission_df['no2_emission'],\n            text=\"\",\n            textposition='auto'\n        )])\n\nfig.show()\n#print (\"All PR Power plants NO2 emission\\n\")","execution_count":null,"outputs":[]},{"metadata":{"id":"h_C1KMmp_eYw","colab_type":"text"},"cell_type":"markdown","source":"### Step 6: Simple Emission Factor Calculation"},{"metadata":{"id":"pmmTnDra_OKY","colab_type":"text"},"cell_type":"markdown","source":"E = A x EF x (1-ER/100)\nwhere:\n\nE = emissions; A = activity rate; EF = emission factor, and ER =overall emission reduction efficiency, %\n\ntherefore\n\nEF = E / [A x (1-ER/100)]\nTo simplify things a bit, I'll reduce that equation to: EF = E / A\n\nSimplified Emissions Factor = Emissions / Activity Rate\nWhich again can be similified to the following:\n\nSimplified Emissions Factor = (Measure of NO2 emissions) / (Quanity of electricity generated)"},{"metadata":{"id":"DWCLMNWu_QA7","colab_type":"code","outputId":"49a95ebf-5858-45c7-adf7-a977ca7f00aa","executionInfo":{"status":"ok","timestamp":1584720693377,"user_tz":-330,"elapsed":2009997,"user":{"displayName":"Meenakshi Ramaswamy","photoUrl":"https://lh4.googleusercontent.com/-3sXdlOXrEr8/AAAAAAAAAAI/AAAAAAAACC8/ngMx9EiU3T4/s64/photo.jpg","userId":"09422112087559126593"}},"colab":{"base_uri":"https://localhost:8080/","height":98},"trusted":true},"cell_type":"code","source":"power_plants_df = power_plants_df.sort_values('capacity_mw',ascending=False).reset_index()\npower_plants_df[['name','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh']][29:30]\nquantity_of_electricity_generated = power_plants_df['estimated_generation_gwh'][29:30].values\nprint('Quanity of Electricity Generated: ', quantity_of_electricity_generated)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### The ground truth assumption on no2 emission , we take from our csv"},{"metadata":{"id":"uqRJCk7F-0rV","colab_type":"code","outputId":"19ab60fa-a5f6-4db0-ec1e-cd6c8d29b24a","executionInfo":{"status":"ok","timestamp":1584720693378,"user_tz":-330,"elapsed":2009958,"user":{"displayName":"Meenakshi Ramaswamy","photoUrl":"https://lh4.googleusercontent.com/-3sXdlOXrEr8/AAAAAAAAAAI/AAAAAAAACC8/ngMx9EiU3T4/s64/photo.jpg","userId":"09422112087559126593"}},"colab":{"base_uri":"https://localhost:8080/","height":136},"trusted":true},"cell_type":"code","source":"power_plants_fossil_df[power_plants_fossil_df.primary_fuel=='Coal']\nannual_EF_coal = power_plants_fossil_df.loc[power_plants_fossil_df['primary_fuel'] == 'Coal', 'no2_emission_in_mols'].values\nprint (\"Emission factor for Coal PP A.E.S. Corp {:.8f} mol * h / m^2 * gw\".format(float(annual_EF_coal/((365*24)))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Coal plant EF is not accurate due to data issues."},{"metadata":{"id":"-_lIpGlAUQaq","colab_type":"text"},"cell_type":"markdown","source":"### Lets see oil powerplant Emission Factor\n\n"},{"metadata":{"id":"PTI-CZrQUU_p","colab_type":"code","outputId":"03c1f670-1816-4c12-e07b-b14ab2230570","executionInfo":{"status":"ok","timestamp":1584720693381,"user_tz":-330,"elapsed":2009839,"user":{"displayName":"Meenakshi Ramaswamy","photoUrl":"https://lh4.googleusercontent.com/-3sXdlOXrEr8/AAAAAAAAAAI/AAAAAAAACC8/ngMx9EiU3T4/s64/photo.jpg","userId":"09422112087559126593"}},"colab":{"base_uri":"https://localhost:8080/","height":465},"trusted":true},"cell_type":"code","source":"power_plants_fossil_df[power_plants_fossil_df.primary_fuel=='Oil']\nannual_EF_oil = power_plants_fossil_df.loc[power_plants_fossil_df['primary_fuel'] == 'Oil', 'no2_emission_in_mols'].values\nprint (\"Emission factor for Vieques EPP {:.8f} mol * h / m^2 * gw\".format(float(annual_EF_oil[3]/(365*24))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Emission factor would be no2 emission in mols divided by the generated electricity (gwh)"},{"metadata":{},"cell_type":"markdown","source":"### The hourly emission of No2 for fossil power plants could be calculated as below"},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_fossil_df['no2_emission(mol / hr)'] = power_plants_fossil_df['no2_emission_in_mols'] / (365*24)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_fossil_df[['capacity_mw','estimated_generation_gwh','name','primary_fuel','input_energy_Btu',\n                        #'capacity_utilization',\n                        'Efficiency_percent',\n                        'no2_Emission_for_input_Btu_in_gms','no2_emission_in_tons', 'no2_emission_in_mols', 'no2_emission(mol / hr)']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Lets calculate the EF for hourly , daily and yearly"},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_fossil_df['EF(mol * h / m^2 * gw)'] = power_plants_fossil_df['no2_emission(mol / hr)'] / power_plants_fossil_df['estimated_generation_gwh'] \npower_plants_fossil_df['EF_daily'] = power_plants_fossil_df['no2_emission(mol / hr)'] / power_plants_fossil_df['estimated_generation_gwh'] * 24\npower_plants_fossil_df['EF_yearly'] = power_plants_fossil_df['no2_emission(mol / hr)'] / power_plants_fossil_df['estimated_generation_gwh'] * (24 *365)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_fossil_df[['capacity_mw','estimated_generation_gwh','name','primary_fuel','input_energy_Btu',#'capacity_utilization',\n                        'Efficiency_percent','no2_Emission_for_input_Btu_in_gms','no2_emission_in_tons', 'no2_emission_in_mols', 'no2_emission(mol / hr)','EF(mol * h / m^2 * gw)','EF_daily','EF_yearly']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_fossil_df[['capacity_mw','estimated_generation_gwh','name','primary_fuel','input_energy_Btu',#'capacity_utilization',\n                        'Efficiency_percent','no2_Emission_for_input_Btu_in_gms','no2_emission_in_tons', 'no2_emission_in_mols', 'no2_emission(mol / hr)','EF(mol * h / m^2 * gw)',\n                        'EF_daily','EF_yearly']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Marginal Emission Factor \n#### MEF would be the ratio of total no2 emission from all power plants corresponding to their efficiency\n**** MEF of plant = NO2 Emission of plant / Efficiency of plant"},{"metadata":{"trusted":true},"cell_type":"code","source":"MEF_df = power_plants_fossil_df\npower_plants_fossil_df['plant_name'] = power_plants_fossil_df['name'] + ':' + power_plants_fossil_df['primary_fuel']\npower_plants_fossil_df = pd.merge(power_plants_fossil_df, emission_df, left_on=('plant_name'), right_on='plant_name')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_fossil_df['MEF'] = power_plants_fossil_df['no2_emission'] / (power_plants_fossil_df['Efficiency_percent'] / 100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#power_plants_fossil_df.head()\nprint ('Marginal Emission Factor : ', power_plants_fossil_df['MEF'].sum())\nprint ('Average No2 Emission : ', power_plants_fossil_df['no2_emission'].sum())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_fossil_df[['capacity_mw','estimated_generation_gwh','name','primary_fuel','input_energy_Btu','capacity_utilization',\n                        'Efficiency_percent','no2_Emission_for_input_Btu_in_gms','no2_emission_in_tons', 'no2_emission_in_mols', 'no2_emission(mol / hr)','EF(mol * h / m^2 * gw)',\n                        'EF_daily','EF_yearly','MEF']]","execution_count":null,"outputs":[]},{"metadata":{"id":"FwIyuwa-2w6o","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"power_plants_fossil_df.to_csv('power_plants_fossil_df.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### We are not concentrating on green power plants as they are not contributing towards no2 emissions"},{"metadata":{},"cell_type":"markdown","source":"###  Models"},{"metadata":{"id":"fpaRCWz9rVKj","colab_type":"text"},"cell_type":"markdown","source":"#### WE will try some time series and (unsupervised )clustering models.\n\n#### To Start with try, *time series*\n\n"},{"metadata":{},"cell_type":"markdown","source":"#### Lets try Autoregression Models for Time Series Forecasting With Python\n\n\nA regression model, such as linear regression, models an output value based on a linear combination of input values.\n\nAn autoregression model makes an assumption that the observations at previous time steps are useful to predict the value at the next time step.\n\nIf both variables change in the same direction (e.g. go up together or down together), this is called a positive correlation. If the variables move in opposite directions as values change (e.g. one goes up and one goes down), then this is called negative correlation\n"},{"metadata":{},"cell_type":"markdown","source":"### ARIMA Model"},{"metadata":{},"cell_type":"markdown","source":"Referred \nhttps://www.digitalocean.com/community/tutorials/a-guide-to-time-series-forecasting-with-arima-in-python-3"},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_dfs_gldas_df.date_from = pd.to_datetime(no2_dfs_gldas_df.date_from)\nno2_weather_values_stats = no2_dfs_gldas_df.groupby('date_from').mean().reset_index()\nprint('We have data for {} days'.format(no2_dfs_gldas_df['date_from'].nunique()))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train= pd.DataFrame()\ntrain['ds'] = pd.to_datetime(no2_weather_values_stats[\"date_from\"])\ntrain['y']=no2_weather_values_stats[\"tropospheric_NO2\"]\nindexedData = train.set_index('ds')\nindexedData.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"moving_average = indexedData.rolling(window=12).mean()\nplt.plot(indexedData, color='blue', label='Original')\nplt.plot(moving_average, color='red', label='Rolling Mean')\nplt.legend(loc='best')\nplt.title('Moving Average-Annual')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# The 'MS' string groups the data in buckets by start of the month\n#train.set_index('y', inplace=True)\ny = indexedData['y'].resample('MS').mean()\n\n# The term bfill means that we use the value before filling in missing values\ny = y.fillna(y.bfill())\n\n#print(y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import itertools\nimport statsmodels.api as sm\npd.plotting.register_matplotlib_converters() # Add this \nplt.style.use('fivethirtyeight')\nf, ax = plt.subplots(figsize=(14,8))\nax.set_xlabel('date_from', fontsize=15)\nax.set_ylabel('N02 values', fontsize=15);\nindexedData['y'].plot(figsize=(15, 6))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"# Define the p, d and q parameters to take any value between 0 and 2\np = d = q = range(0, 2)\n\n# Generate all different combinations of p, q and q triplets\npdq = list(itertools.product(p, d, q))\n\n# Generate all different combinations of seasonal p, q and q triplets\nseasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]\n\nprint('Examples of parameter combinations for Seasonal ARIMA...')\nprint('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))\nprint('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))\n#print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))\n#print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"warnings.filterwarnings(\"ignore\") # specify to ignore warning messages\n\nfor param in pdq:\n    for param_seasonal in seasonal_pdq:\n        try:\n            mod = sm.tsa.statespace.SARIMAX(indexedData['y'],\n                                            order=param,\n                                            seasonal_order=param_seasonal,\n                                            enforce_stationarity=False,\n                                            enforce_invertibility=False)\n\n            results = mod.fit()\n\n            print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))\n        except:\n            continue","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mod = sm.tsa.statespace.SARIMAX(indexedData['y'],\n                                order=(4, 2, 1),\n                                seasonal_order=(1, 1, 1, 12),\n                                enforce_stationarity=False,\n                                enforce_invertibility=False)\n\nresults = mod.fit()\n\nprint(results.summary())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The coef column shows the weight (i.e. importance) of each feature and how each one impacts the time series. The P>|z| column informs us of the significance of each feature weight. Here, each weight has a p-value close to 0, so it is reasonable to include the features in our model."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"results.plot_diagnostics(figsize=(10, 8))\n#plt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred = results.get_prediction(start=270,end =330, dynamic=False)\npred_ci = pred.conf_int()\npred_ci.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ax = y['2018':].plot(label='observed')\npred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)\n\nax.fill_between(pred_ci.index,\n                pred_ci.iloc[:, 0],\n                pred_ci.iloc[:, 1], color='r', alpha=.2)\n\nax.set_xlabel('Date')\nax.set_ylabel('NO2 Levels')\nplt.legend()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"mte_forecast = pred.predicted_mean\n#print (indexedData['2019-03-05':])\nmte_truth = indexedData['2019-03-05':]\n\n\n# Compute the mean square error\nmse = ((mte_forecast - mte_truth) ** 2).mean()\nprint('The Mean Squared Error (MSE) of the forecast is {}'.format(round(mse, 2)))\n#print('The Root Mean Square Error (RMSE) of the forcast: {:.4f}'.format(np.sqrt(sum((mte_forecast-mte_truth)**2)/len(mte_forecast))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"### Prophet"},{"metadata":{},"cell_type":"markdown","source":"#### Lets try Prophet model to find the trend and pattern"},{"metadata":{"trusted":true},"cell_type":"code","source":"from fbprophet import Prophet","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Reference kernel\nhttps://www.kaggle.com/armamut/predicting-transactions-fb-prophet-tutorial"},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"m = Prophet()\nm.fit(train)\nfuture = m.make_future_dataframe(periods=100)\nforecast = m.predict(future)\nforecast","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"m.plot_components(forecast)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### We could see the trend, monday and thursday are *peak No2* days. "},{"metadata":{"trusted":true},"cell_type":"code","source":"import plotly.offline as py\nimport plotly.graph_objs as go\npy.init_notebook_mode()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"py.iplot([\n    go.Scatter(x=train['ds'], y=train['y'], name='y'),\n    go.Scatter(x=forecast['ds'], y=forecast['yhat'], name='yhat'),\n    go.Scatter(x=forecast['ds'], y=forecast['yhat_upper'], fill='tonexty', mode='none', name='upper'),\n    go.Scatter(x=forecast['ds'], y=forecast['yhat_lower'], fill='tonexty', mode='none', name='lower'),\n    go.Scatter(x=forecast['ds'], y=forecast['trend'], name='Trend')\n])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"m = Prophet(changepoint_prior_scale=2.5)\nm.fit(train)\nfuture = m.make_future_dataframe(periods=180)\nforecast = m.predict(future)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Calculate root mean squared error.\nprint('RMSE: %f' % np.sqrt(np.mean((forecast.loc[:100, 'yhat']-train['y'])**2)) )\npy.iplot([\n    go.Scatter(x=train['ds'], y=train['y'], name='y'),\n    go.Scatter(x=forecast['ds'], y=forecast['yhat'], name='yhat'),\n    go.Scatter(x=forecast['ds'], y=forecast['yhat_upper'], fill='tonexty', mode='none', name='upper'),\n    go.Scatter(x=forecast['ds'], y=forecast['yhat_lower'], fill='tonexty', mode='none', name='lower'),\n    go.Scatter(x=forecast['ds'], y=forecast['trend'], name='Trend')\n])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### We can see there is sudden surge in no2 density during April 2019 and Sep 2018. "},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"### Lets try with the clustering alogorithm k-means"},{"metadata":{"trusted":true},"cell_type":"code","source":"def overlay_image_on_puerto_rico_df(df, img, zoom):\n    lat_map=df.iloc[[0]].loc[:,[\"latitude\"]].iat[0,0]\n    lon_map=df.iloc[[0]].loc[:,[\"longitude\"]].iat[0,0]\n    m = folium.Map([lat_map, lon_map], zoom_start=zoom)\n    color={ 'Hydro' : 'lightblue', 'Solar' : 'orange', 'Oil' : 'darkblue', 'Coal' : 'black', 'Gas' : 'lightgray', 'Wind' : 'green' }\n    folium.raster_layers.ImageOverlay(\n        image=img,\n        bounds = [[18.56,-67.32,],[17.90,-65.194]],\n        colormap=lambda x: (1, 0, 0, x),\n    ).add_to(m)\n    \n    for i in range(0,len(df)):\n        popup = folium.Popup(str(df.primary_fuel[i:i+1]))\n        folium.Marker([df[\"latitude\"].iloc[i],df[\"longitude\"].iloc[i]],\n                     icon=folium.Icon(icon_color='red',icon ='bolt',prefix='fa',color=color[df.primary_fuel.iloc[i]])).add_to(m)\n        \n    return m","execution_count":null,"outputs":[]},{"metadata":{"id":"k6kAR8HLlYIQ","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"image = data_path + '/eie_data/s5p_no2/s5p_no2_20190510T164338_20190516T183223.tif'\n#overlay_image_on_puerto_rico(file_name,band_layer,lat,lon,zoom):\noverlay_image_on_puerto_rico_df(power_plants_df,image,zoom=8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.cluster import KMeans\nfrom sklearn.preprocessing import MinMaxScaler\nscaler = MinMaxScaler() ","execution_count":null,"outputs":[]},{"metadata":{"id":"BIDRRKEEl0lG","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"lon = []\nlat = []\nNO2 = []\nimp_img = np.flip(tiff.imread(image))\nfor i in range(imp_img[:,:,1].shape[0]):\n    for j in range(imp_img[:,:,1].shape[1]):\n        #print(imp_img[:,:,1][i,j])\n        NO2.append(imp_img[:,:,1][i,j])\n        lon.append(i)\n        lat.append(j)\n        \nNO2 = np.array(NO2)\nlon = np.array(lon)\nlat = np.array(lat)","execution_count":null,"outputs":[]},{"metadata":{"id":"-IDL7vOIl9NV","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"results = pd.DataFrame(columns=['NO2', 'lat', 'lon'])\nresults = pd.DataFrame({'NO2': NO2/max(NO2),\n                    'lat': lat/max(lat),\n                    'lon': lon/max(lon)})","execution_count":null,"outputs":[]},{"metadata":{"id":"nHtc1Y_TmCg2","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"data_scaled = scaler.fit_transform(results)\npred = KMeans(n_clusters=15).fit_predict(results)","execution_count":null,"outputs":[]},{"metadata":{"id":"NjAOww0VmEAo","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"plt.figure()\npred_image = pred.reshape(148, 475)\nsns.heatmap(pred_image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"Error =[]\nfor i in range(1, 15):\n    kmeans = KMeans(n_clusters = i).fit(results)\n    kmeans.fit(results)\n    Error.append(kmeans.inertia_)\n\nplt.plot(range(1, 15), Error)\nplt.title('Elbow method')\nplt.xlabel('No of clusters')\nplt.ylabel('Error')\nplt.show()\nkmeans.cluster_centers_","execution_count":null,"outputs":[]},{"metadata":{"id":"zr_o43tkmHZo","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"overlay_image_on_puerto_rico_df(power_plants_df,pred_image,zoom=8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### we applied the idea from [this great kernel](https://www.kaggle.com/nayuts/can-we-attribute-emissions-to-power-plants)\n### changing the dates for our analysis and adding extra parameters temperature spechumidity relhumidity  totalprecipitationsurface"},{"metadata":{},"cell_type":"markdown","source":"### We take data from Apr-12-2019 till Apr-18-2019,for a week ."},{"metadata":{"trusted":true},"cell_type":"code","source":"gldas_files = glob.glob('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gldas/*')\ngldas_files = sorted(gldas_files)\ngfs_files = glob.glob('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/*')\ngfs_files = sorted(gfs_files)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#(gldas_files[2280:2329])\n#gfs_files[1140:1165]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gldas_files_daily = []\nfor i in range(0,len(gldas_files[2280:2329]),8):\n    #print(gldas_files[i:i+8])\n    gldas_files_daily.append(gldas_files[i:i+8])\n\ngfs_files_daily = []\nfor i in range(0,len(gfs_files[1140:1165]),4):\n    #print(gfs_files[i:i+4])\n    gfs_files_daily.append(gfs_files[i:i+4])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"temperature spechumidity relhumidity uwind vwind totalprecipitationsurface"},{"metadata":{"trusted":true},"cell_type":"code","source":"image_reglession_u = []\nimage_reglession_v = []\nimage_reglession_t = []\nimage_reglession_sh = []\nimage_reglession_rh = []\nimage_reglession_tp = []\nimage_reglession_rn = []\n\nfor i in range(len(gfs_files_daily)):\n    gfs_tmp = gfs_files_daily[i]\n    gldas_tmp = gldas_files_daily[i]\n    array_wind_u = []\n    array_wind_v = []\n    array_t = []\n    array_sh = []\n    array_rh = []\n    array_tp = []\n    array_rn = []\n    for j in range(len(gfs_tmp)):\n        gfs_image_t = tiff.imread(gfs_tmp[j])[:,:,0]\n        gfs_image_sh = tiff.imread(gfs_tmp[j])[:,:,1]\n        gfs_image_rh = tiff.imread(gfs_tmp[j])[:,:,2]\n        gfs_image_u = tiff.imread(gfs_tmp[j])[:,:,3]\n        gfs_image_v = tiff.imread(gfs_tmp[j])[:,:,4]\n        gfs_image_tp = tiff.imread(gfs_tmp[j])[:,:,5]        \n        gldas_image_rn = tiff.imread(gldas_tmp[2*j])[:,:,7]\n        gldas_image1 = tiff.imread(gldas_tmp[2*j])[:,:,11]\n        gldas_image2 = tiff.imread(gldas_tmp[2*j + 1])[:,:,11]\n\n        #fill na by mean\n        gfs_image_t = np.nan_to_num(gfs_image_t, nan=np.nanmean(gfs_image_t))\n        gfs_image_sh = np.nan_to_num(gfs_image_sh, nan=np.nanmean(gfs_image_sh))\n        gfs_image_rh= np.nan_to_num(gfs_image_rh, nan=np.nanmean(gfs_image_rh))\n        gfs_image_tp = np.nan_to_num(gfs_image_tp, nan=np.nanmean(gfs_image_tp))\n        \n        gfs_image_u = np.nan_to_num(gfs_image_u, nan=np.nanmean(gfs_image_u))\n        gfs_image_v = np.nan_to_num(gfs_image_v, nan=np.nanmean(gfs_image_v))\n        \n        gldas_image_rn = np.nan_to_num(gldas_image_rn, nan=np.nanmean(gldas_image_rn))\n        gldas_image1 = np.nan_to_num(gldas_image1, nan=np.nanmean(gldas_image1))\n        gldas_image2 = np.nan_to_num(gldas_image2, nan=np.nanmean(gldas_image2))\n       \n        gldas_image = (gldas_image1 + gldas_image2)/2\n        wind_u = gfs_image_u * gldas_image\n        wind_v = gfs_image_v * gldas_image\n       \n        \n        array_wind_u.append(wind_u)\n        array_wind_v.append(wind_v)\n        \n        array_t.append(gfs_image_t)\n        array_sh.append(gfs_image_sh)\n        array_rh.append(gfs_image_rh)\n        array_tp.append(gfs_image_tp)\n\n        array_rn.append(gldas_image_rn)\n        \n        image_reglession_u.append(np.nanmean(np.array(array_wind_u), axis=0))\n        image_reglession_v.append(np.nanmean(np.array(array_wind_v), axis=0))\n        image_reglession_rn.append(np.nanmean(np.array(array_rn), axis=0))\n        \n        image_reglession_t.append(np.nanmean(np.array(array_t), axis=0))\n        image_reglession_sh.append(np.nanmean(np.array(array_sh), axis=0))\n        image_reglession_tp.append(np.nanmean(np.array(array_tp), axis=0))\n        image_reglession_rh.append(np.nanmean(np.array(array_rh), axis=0))\n\n\n        \nimage_reglession_u = np.nanmean(np.array(image_reglession_u), axis=0)\nimage_reglession_v = np.nanmean(np.array(image_reglession_v), axis=0)\nimage_reglession_rn = np.nanmean(np.array(array_rn), axis=0)\n\nimage_reglession_t = np.nanmean(np.array(array_t), axis=0)\nimage_reglession_sh = np.nanmean(np.array(array_sh), axis=0)\nimage_reglession_tp = np.nanmean(np.array(array_tp), axis=0)\nimage_reglession_rh = np.nanmean(np.array(array_rh), axis=0)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_u.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_v.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_t.reshape((148, 475)))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_sh.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_rh.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = tiff.imread('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20190412T170633_20190418T190433.tif')\nlon = []\nlat = []\nNO2 = []\nwind_u = []\nwind_v = []\n\ntemp = []\nspechum = []\nrelhum = []\ntotper = []\nrain = []\n\n\nfor i in range(image[:,:,0].shape[0]):\n    for j in range(image[:,:,0].shape[1]):\n        #print(image[:,:,0][i,j])\n        NO2.append(image[:,:,0][i,j])\n        lon.append(i)\n        lat.append(j)\n        wind_u.append(image_reglession_u.reshape((148, 475))[i,j])\n        wind_v.append(image_reglession_v.reshape((148, 475))[i,j])\n        temp.append(image_reglession_u.reshape((148, 475))[i,j])\n        spechum.append(image_reglession_v.reshape((148, 475))[i,j])\n        relhum.append(image_reglession_u.reshape((148, 475))[i,j])\n        totper.append(image_reglession_v.reshape((148, 475))[i,j])\n        rain.append(image_reglession_u.reshape((148, 475))[i,j])\n\n\n        \nNO2 = np.array(NO2)\nlon = np.array(lon)\nlat = np.array(lat)\nwind_u = np.array(wind_u)\nwind_v = np.array(wind_v)\ntemp = np.array(temp)\nspechum = np.array(spechum)\nrelhum = np.array(relhum)\ntotpert = np.array(totper)\nrain = np.array(rain)\n        \nresults_wind = pd.DataFrame(columns=['NO2', 'lat', 'lon', 'wind_u', 'wind_v','temp','spechum' ,'relhum','totper','rain'])\nresults_wind = pd.DataFrame({\n                    'NO2': NO2/max(NO2),\n                    'lat': lat/max(lat),\n                    'lon': lon/max(lon),\n                    'wind_u' : wind_u/(- min(wind_u)),\n                    'wind_v' : wind_v/(- min(wind_v)),\n                    'temp': temp/max(temp),\n                    'spechum': spechum/max(spechum),\n                    'relhum': relhum/max(relhum),\n                    'totper': totper/max(totper) #,'rain': rain/max(rain)\n                    })","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results_wind[\"NO2\"].fillna(0, inplace = True) \nresults_wind[\"temp\"].fillna(0, inplace = True) \nresults_wind[\"spechum\"].fillna(0, inplace = True) \nresults_wind[\"relhum\"].fillna(0, inplace = True) \nresults_wind[\"totper\"].fillna(0, inplace = True) \n#results_wind[\"rain\"].fillna(0, inplace = True) \ndata_scaled = scaler.fit_transform(results_wind)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_all_factors = KMeans(n_clusters=15).fit_predict(results_wind)\nplt.figure()\nsns.heatmap(pred_all_factors.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"overlay_image_on_puerto_rico_df(power_plants_df, pred_all_factors.reshape((148, 475)), 12)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Northern and southern coast are having high no2 as we have coal plant in the south coast and cluster of oil plants in the north coast"},{"metadata":{},"cell_type":"markdown","source":"#### Lets try with some supervised models"},{"metadata":{},"cell_type":"markdown","source":"#### Simple Linear Regression"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nfrom sklearn import metrics\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import r2_score","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x= no2_weather_values_stats[['temp_2m_above_ground_mean','sp_humidity_2m_above_ground_mean',\n                     'rh_2m_above_ground_mean','wind_velocity',#'date_from','tropospheric_NO2','no2_mean',\n                     'tot_percip_surface_mean', 'aai_mean','tropopause_pressure_mean','cloud_fraction_mean', 'Qair_f_inst_mean',\n                     'Tair_f_inst_mean','Wind_f_inst_mean']]\ny = no2_weather_values_stats['no2_mean']\nx_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.25, random_state=2020)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"lm = LinearRegression()\nlm_fit = lm.fit(x_train, y_train)\nlm_train_fit = lm_fit.predict(x_train)\ny_pred_lm = lm_fit.predict(x_test)\n#print(\"r2 score on Test:\", r2_score(y_test,y_pred_lm))\nlm.score(x_test, y_test)\nprint('RMSE :', np.sqrt(metrics.mean_squared_error(y_test, y_pred_lm)))\nprint('MSE :', metrics.mean_squared_error(y_test, y_pred_lm))\n\ncoeff_df = pd.DataFrame(lm.coef_, x.columns, columns=['Coefficient'])\ncoeff_df\ndf = pd.DataFrame({'Actual': y_test, 'Predicted': y_pred_lm, 'difference': y_test - y_pred_lm})\ndf.head()\n#lm.intercept_","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Import Eli5 package\nimport eli5\nfrom eli5.sklearn import PermutationImportance\n\n# Find the importance of columns for prediction\nperm = PermutationImportance(lm_fit, random_state=2020).fit(x_test,y_pred_lm)\neli5.show_weights(perm, feature_names = x_test.columns.tolist())\n#Understanding how each feature influences the prediction\neli5.show_prediction(lm_fit, doc=x_test.iloc[[5]], feature_names=list(x_test.columns))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#Import SHAP package\nimport shap\n\n#Create explainer for linear model\nexplainer = shap.LinearExplainer(lm_fit,data=x_test.values)\nshap_values = explainer.shap_values(x_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#Understanding how each feature influences the prediction\n\nshap.initjs()\nind = 55\n\n\nshap.force_plot(\n    explainer.expected_value, shap_values[ind,:], x_test.iloc[ind,:],\n    feature_names=x_test.columns.tolist()\n)\n\nshap.summary_plot(shap_values,x_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Definitiely temperature and air plays major role in no2 density"},{"metadata":{},"cell_type":"markdown","source":"#### Lets try with xgb"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from xgboost import XGBRegressor\nfrom sklearn.model_selection import GridSearchCV\ncv_params = {'max_depth': [3,5,7], 'min_child_weight': [1,3,5]}\nind_params = {'eta': 0.1, \n              'eval_metric': 'rmse',\n              'n_estimators': 1000,\n              'seed':0,\n              'subsample': 0.8,\n              'colsample_bytree': 0.8, \n             'objective': 'reg:linear'}\noptimized_GBM = GridSearchCV(XGBRegressor(**ind_params), \n                            cv_params, cv = 5, n_jobs = -1) \ncv_res = optimized_GBM.fit(x_train,y_train)\noptimized_GBM.cv_results_","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Random Forest"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from sklearn.ensemble import RandomForestRegressor\n\nforest = RandomForestRegressor(max_depth = 10, n_estimators = 1000, random_state = 2020)\nmodel_rf = forest.fit(x_train, y_train)\nprint(model_rf.score(x_train, y_train))\n# Making predictions\ny_pred_rf = model_rf.predict(x_test)\nlist(zip(x_train.columns,model_rf.feature_importances_))\nmodel_rf.score(x_test,y_test)\nprint('RMSE :',np.sqrt(metrics.mean_squared_error(y_test,y_pred_rf)))\nprint('MSE :', metrics.mean_squared_error(y_test, y_pred_rf))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## End Note"},{"metadata":{},"cell_type":"markdown","source":"*** emission_df - Provided the Power plantwise annual emission\n\n*** pp_monthwise - Provides the monthly emission on each power plant\n\n*** Marginal emission factor - Included on the power_plant_fossil_df\n\n*** power_plant_fossil_df - Provided the emission factor of fossil powerplants"},{"metadata":{},"cell_type":"markdown","source":"### Thanks to kaggle organizers and community for providing such a different environmental project, it was really challenging project and great opportunity to learn."},{"metadata":{},"cell_type":"markdown","source":"### Acknowledgement\n\n\n"},{"metadata":{"id":"0VVg4Exb4b0I","colab_type":"text"},"cell_type":"markdown","source":"### Starter kernels\n\nhttps://www.kaggle.com/paultimothymooney/how-to-get-started-with-the-earth-engine-data/comments\n\nhttps://www.kaggle.com/paultimothymooney/overview-of-the-eie-analytics-challenge\n\nhttps://www.kaggle.com/paultimothymooney/explore-image-metadata-s5p-gfs-gldas"},{"metadata":{"id":"lj-9v_eponCs","colab_type":"text"},"cell_type":"markdown","source":"### Referred kernels & Discussion Links\n"},{"metadata":{},"cell_type":"markdown","source":"https://www.kaggle.com/gpoulain/eda-ef-with-n2o-time-series-earth-engine\n\nhttps://www.kaggle.com/caesarlupum/ds4g-anomaly-analysis\n\nhttps://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson\n\nhttps://www.kaggle.com/ragnar123/exploratory-data-analysis-and-factor-model-idea\n\nhttps://www.kaggle.com/nayuts/can-we-attribute-emissions-to-power-plants\n\nhttps://www.kaggle.com/parulpandey/understanding-the-data-wip\n\n\nhttps://www.kaggle.com/c/ds4g-environmental-insights-explorer/discussion/129991\n\nhttps://www.kaggle.com/c/ds4g-environmental-insights-explorer/discussion/134727\n\nhttps://www.kaggle.com/c/ds4g-environmental-insights-explorer/discussion/130221\n\nhttps://www.kaggle.com/c/ds4g-environmental-insights-explorer/discussion/130055\n\nhttps://www.kaggle.com/c/ds4g-environmental-insights-explorer/discussion/130537"},{"metadata":{"id":"AdekdG1i5SlV","colab_type":"text"},"cell_type":"markdown","source":"#### Thank you for reading !!! Suggestions and comments are welcome."}],"metadata":{"colab":{"name":"DS4G-Part-3.ipynb","provenance":[{"file_id":"1e_y8GV4f_i2cdBmHrln2LdJRucSccTfR","timestamp":1584533637402},{"file_id":"1uSRpNrAyduJ8RuWYd9o0RDVfsfLvBS19","timestamp":1584279176959}],"collapsed_sections":[],"machine_shape":"hm"},"kernelspec":{"name":"python3","display_name":"Python 3"},"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.6"},"accelerator":"GPU"},"nbformat":4,"nbformat_minor":4}