{"cells":[{"metadata":{},"cell_type":"markdown","source":"## The Engie Factor Team NOTEBOOK"},{"metadata":{},"cell_type":"markdown","source":"\n"},{"metadata":{},"cell_type":"markdown","source":"## Authors\nName, Kaggle ID, email
\nSergio Garcia, SergioGarciadea, QK5945@engie.com
\nLouis Giron, LGiron, Louis.Giron@engie.com
\nAbderrahim Benseguane,AbderrahimBenseguane, abderrahim.benseguane@engie.com
\n## CREDITS:\n### The following example Jupyter Notebooks helped get us get started:\n#### https://www.kaggle.com/paultimothymooney/overview-of-the-eie-analytics-challenge\n#### https://www.kaggle.com/parulpandey/understanding-the-data-wip\n#### https://www.kaggle.com/paultimothymooney/how-to-get-started-with-the-earth-engine-data\n#### https://www.kaggle.com/ragnar123/exploratory-data-analysis-and-factor-model-idea\n\n## Useful links:\n#### https://www.atmos-chem-phys.net/16/5283/2016/acp-16-5283-2016.pdf'"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"from scipy import integrate\nfrom IPython.display import Math, Latex, display, HTML\n\nfrom matplotlib import animation\nfrom matplotlib import pyplot as plt\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nimport datetime as dt\nfrom datetime import datetime \n\ntry:\n import rasterio as rio\n import folium\nexcept Exception:\n pass\nimport tifffile as tiff \nimport random \n\nimport os\nfrom tqdm import tqdm\nimport ast\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\n%matplotlib inline\nimport seaborn as sns\n\ndef plot_points_on_map(dataframe,begin_index,end_index,latitude_column,latitude_value,longitude_column,longitude_value,zoom):\n df = dataframe[begin_index:end_index]\n location = [latitude_value,longitude_value]\n plot = folium.Map(location=location,zoom_start=zoom)\n for i in range(0,len(df)):\n popup = folium.Popup(str(df.primary_fuel[i:i+1]))\n folium.Marker([df[latitude_column].iloc[i],df[longitude_column].iloc[i]],popup=popup).add_to(plot)\n return(plot)\n\ndef overlay_puerto_rico(file_name,band_layer,lat,lon,zoom, bnds=[[18.6,-67.3,],[17.9,-65.7]]):\n band = rio.open(file_name).read(band_layer)\n band_below = rio.open(file_name).read(band_layer-1)\n m = folium.Map([lat, lon], zoom_start=zoom)\n folium.raster_layers.ImageOverlay(\n image=band,\n bounds = bnds,\n colormap=lambda x: (1, 0, 0, x),\n ).add_to(m)\n\n return m\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"WEATHER_images = []\n\ndata_dir='/kaggle/input/ds4g-environmental-insights-explorer/eie_data/'\n\nimport os\nfor dirname, _, filenames in os.walk(data_dir + 'gfs/'):\n for filename in filenames:\n WEATHER_images.append(os.path.join(dirname, filename))\n\nprint(len(WEATHER_images))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NO2_images = []\n\nimport os\nfor dirname, _, filenames in os.walk(data_dir + 's5p_no2/'):\n for filename in filenames:\n NO2_images.append(os.path.join(dirname, filename))\n\nprint(len(NO2_images))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example = tiff.imread(NO2_images[89])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"latitude=18.1429005246921; longitude=-65.4440010699994\noverlay_puerto_rico('/kaggle/input/my-backgrounds/division.png',band_layer=4,lat=latitude,lon=longitude,zoom=8,\n )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"WIND_data = pd.DataFrame()\nfor i, num in tqdm(enumerate(range(len(WEATHER_images))), total=1460):\n img = tiff.imread(WEATHER_images[i])\n WIND_data[pd.to_datetime(WEATHER_images[i][-14:-6])+dt.timedelta(hours=int(WEATHER_images[i][-6:-4]))] = pd.Series(\n {\n 'path':WEATHER_images[i],'WIND_X_M_SEC':np.average(np.nan_to_num(img[:, :, 3])), \n 'WIND_Y_M_SEC':np.average(np.nan_to_num(img[:, :, 4]))\n }\n )\n #print(dt.timedelta(hours=int(WEATHER_images[i][-6:-4])))\nWIND_data = WIND_data.T\nWIND_data = WIND_data.sort_index()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### We are going to focus on bands 3 and 4, so that we centered around 600 nm\n### Source: Wikipedia, https://en.m.wikipedia.org/wiki/File:Mid-infrared_absorption_spectra_of_Gases.png\n## ---------\n![](https://upload.wikimedia.org/wikipedia/commons/thumb/3/34/Mid-infrared_absorption_spectra_of_Gases.png/440px-Mid-infrared_absorption_spectra_of_Gases.png)\n## ---------"},{"metadata":{"trusted":true},"cell_type":"code","source":"# NOTE: We are using nansum, given that sometimes the images had some NaN values in these band layers\n\nNO2_data = pd.DataFrame()\nfor i, num in tqdm(enumerate(range(len(NO2_images))), total=387):\n img = tiff.imread(NO2_images[i])\n NO2_data[pd.to_datetime(NO2_images[i][-19:-11]) + dt.timedelta(hours=int(NO2_images[i][-10:-8]))] = pd.Series(\n {'path':NO2_images[i],'NO2':np.nansum(img[:, :, 2 : 4]),\n 'NW_NO2':np.nansum(img[:int(example.shape[0]/2), :int(example.shape[1]/3), 2:4]),\n 'SW_NO2':np.nansum(img[int(example.shape[0]/2):, :int(example.shape[1]/3), 2:4]),\n 'N_NO2':np.nansum(img[:int(example.shape[0]/2), int(example.shape[1]/3):int(example.shape[1]*2/3), 2:4]),\n 'S_NO2':np.nansum(img[int(example.shape[0]/2):, int(example.shape[1]/3):int(example.shape[1]*2/3), 2:4]),\n 'NE_NO2':np.nansum(img[:int(example.shape[0]/2), int(example.shape[1]*2/3):, 2:4]),\n 'SE_NO2':np.nansum(img[int(example.shape[0]/2):, int(example.shape[1]*2/3):, 2:4])}\n )\n\nNO2_data = NO2_data.T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set_style(\"darkgrid\") #Available styles: darkgrid, whitegrid, dark, white, and ticks","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(14,11))\nax = NO2_data[['NO2']].plot(figsize=(14,11))\nax.set(xlabel='Date', ylabel='NO2 level')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NO2_data = NO2_data.sort_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Lets compare a couple of regions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(14,11))\nax = NO2_data[['NE_NO2', 'SW_NO2']].plot(figsize=(14,11))\nax.set(xlabel='Date', ylabel='NO2 level')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NO2_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"graphs = ['NW_NO2','SW_NO2','N_NO2','S_NO2','NE_NO2','SE_NO2']\nfig, ax = plt.subplots(2, 3, sharex='col', figsize=(17,9))\n\nfor i in range(2):\n for j in range(3):\n ax_tmp = plt.subplot(ax[i, j])\n ax_tmp.set(xlabel='Date', ylabel='NO2 level')\n ax_tmp.plot(NO2_data[graphs[i+j]])\n ax_tmp.text(0.05,0.95,str(graphs[i+j]),fontsize=11, color='grey',transform=ax_tmp.transAxes)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2, 3, sharex='col', figsize=(17,9))\n\nfor i in range(2):\n for j in range(3):\n ax_tmp = plt.subplot(ax[i, j])\n ax_tmp.set(xlabel='Date', ylabel='NO2 daily variation')\n ax_tmp.plot(NO2_data[graphs[i+j]].diff())\n ax_tmp.text(0.05,0.95,str(graphs[i+j])+' Delta',fontsize=11, color='grey',transform=ax_tmp.transAxes)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Heatmap of the different region's NO2 level correlations\n\nVar_Corr = NO2_data.fillna(0).corr()\nsns.heatmap(Var_Corr, xticklabels=Var_Corr.columns, yticklabels=Var_Corr.columns, annot=True, cmap = 'coolwarm')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for g in graphs:\n NO2_data[g+'_diff'] = NO2_data[g].diff()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Heatmap of the different region's NO2 level variations correlations\nVar_Corr = NO2_data[['NW_NO2_diff','SW_NO2_diff','N_NO2_diff','S_NO2_diff','NE_NO2_diff','SE_NO2_diff']].fillna(0).corr()\nsns.heatmap(Var_Corr, xticklabels=Var_Corr.columns, yticklabels=Var_Corr.columns, annot=True, cmap = 'coolwarm')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NO2_data.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NO2_data.median()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NO2_data.var()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Loading the Power Plant data frame\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"#Power plants in Puerto Rico\n\npower_df = pd.read_csv(data_dir+'gppd/gppd_120_pr.csv')\n#power_df.head(3).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_df.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"# Removing non-polluting Hydro, Solar, and Wind"},{"metadata":{"trusted":true},"cell_type":"code","source":"dirty_power_df = power_df[(power_df.primary_fuel != 'Hydro')&(power_df.primary_fuel != 'Solar') & (power_df.primary_fuel !='Wind')]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dirty_power_df.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Total Installed Capacity - Polluting Capacity\n### Removing renewable plans that do not contribute to emissions"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Total capacity of all the plants\ndirty_capacity_mw = dirty_power_df['capacity_mw'].sum()\nprint('Total Installed Dirty Capacity: '+'{:.2f}'.format(dirty_capacity_mw) + ' MW')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dirty_capacity = (dirty_power_df.groupby(['primary_fuel'])['capacity_mw'].sum()).to_frame()\ndirty_capacity = dirty_capacity.sort_values('capacity_mw',ascending=False)\ndirty_capacity['percentage_of_total'] = (dirty_capacity['capacity_mw']/dirty_capacity_mw)*100\ndirty_capacity","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Estimated generation"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Total generation of all the plants\ntotal_gen_mw = dirty_power_df['estimated_generation_gwh'].sum()\nprint('Total Generatation: '+'{:.2f}'.format(total_gen_mw) + ' GW')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"generation = (dirty_power_df.groupby(['primary_fuel'])['estimated_generation_gwh'].sum()).to_frame()\ngeneration = generation.sort_values('estimated_generation_gwh',ascending=False)\ngeneration['percentage_of_total'] = (generation['estimated_generation_gwh']/total_gen_mw)*100\ngeneration","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Wind Information\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"WIND_data.head(5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Applying an exponential moving average function to the various wind speed measurements to get a better idea of the Average Wind Speed"},{"metadata":{"trusted":true},"cell_type":"code","source":"WIND_data['WIND_X_M_SEC'] = WIND_data['WIND_X_M_SEC'].ewm(com=0.67).mean()\nWIND_data['WIND_Y_M_SEC'] = WIND_data['WIND_Y_M_SEC'].ewm(com=0.67).mean()\n#df['Wind_Speed_X'].rolling(4).apply(integrate.trapz)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def wind_to_X_pixels(wind_speed):\n kmd = (3600 / 1000.0 * 6.0)/2.0 # converting to average km / period displacement\n pto_rico_x = 200 # the island is approx. 200 km wide\n img_x = 475\n response = ((wind_speed * kmd)/pto_rico_x)*img_x\n return response\n \ndef wind_to_Y_pixels(wind_speed):\n kmd = (3600 / 1000.0 * 6.0)/2.0 # converting to average km / period displacement\n pto_rico_y = 75 # the island's height is approx. 75 km\n img_y = 148\n response = ((wind_speed * kmd)/pto_rico_y)*img_y \n return response","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"WIND_data['DAILY_DISP_X_PX'] = WIND_data['WIND_X_M_SEC'].apply(wind_to_X_pixels)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"WIND_data['DAILY_DISP_Y_PX'] = WIND_data['WIND_Y_M_SEC'].apply(wind_to_Y_pixels)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"WIND_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"WIND_data.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### An interesting thing to note is that wind in the East to West direction is very constant, but in the vertical direction component it sometimes goes North to South and sometimes South to North."},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data = pd.concat([NO2_data, WIND_data], axis=1, join='inner')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Building the simulation DataFrame"},{"metadata":{},"cell_type":"markdown","source":"## Simple Generative Model\n
\n“What I cannot create, I do not understand”\n\n
\n"},{"metadata":{},"cell_type":"markdown","source":"### We are going to simulate particles being emitted everyday in the morning proportional to the daily energy generated by that plant times the Emission Factor for that type of energy (we are only working with dirty types of energy). We will then calculate the correlation between estimated pollution in the simulation and measured NO2 in each sector. Finally we will repeat the process varying the emission factors for each of the dirty fuels, seeking to minize the following formula:"},{"metadata":{"trusted":true},"cell_type":"code","source":"display(Math(r'\\min F(C^F,O^F,G^F) = 1- \\sqrt{\\frac{corr_{NW}^2 + corr_{NE}^2 + corr_{N}^2 …}{6}}'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### The emission factors that minimizes the above are our proposed emission factor (EF) estimates "},{"metadata":{"trusted":true},"cell_type":"code","source":"def convert_geo_px(long, lat):\n\n\n image = data_dir + 's5p_no2/s5p_no2_20180701T161259_20180707T175356.tif'\n sat_data = rio.open(image)\n\n width_in_projected_units = sat_data.bounds.right - sat_data.bounds.left\n height_in_projected_units = sat_data.bounds.top - sat_data.bounds.bottom\n\n # Upper left pixel\n row_min = 0\n col_min = 0\n # Lower right pixel. Rows and columns are zero indexing.\n row_max = sat_data.height - 1\n col_max = sat_data.width - 1\n # Transform coordinates with the dataset's affine transformation.\n topleft = sat_data.transform * (row_min, col_min)\n botright = sat_data.transform * (row_max, col_max)\n\n return sat_data.width*((long-topleft[0])/width_in_projected_units),\\\n sat_data.height*((topleft[1]-lat)/height_in_projected_units)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sim_df = pd.DataFrame(index=all_data.index, columns = ['NW_NO2', 'SW_NO2','N_NO2','S_NO2',\n 'NE_NO2','SE_NO2'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sim_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emission_sources = dirty_power_df[['primary_fuel', 'estimated_generation_gwh' ,'.geo']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emission_sources['estimated_generation_gwh']/=365","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emission_sources = emission_sources.rename(columns={\"estimated_generation_gwh\": \"daily_gwh\"})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emission_sources.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emission_sources['.geo'] = emission_sources['.geo'].apply(lambda x: ast.literal_eval(x))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emission_sources['.geo'] = emission_sources['.geo'].apply(lambda x: x.get('coordinates'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emission_sources['source_px'] = emission_sources['.geo'].apply(lambda x: convert_geo_px(x[0],x[1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"emission_sources","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def in_interval(x, interval):\n if x < interval[1] and x >= interval[0]:\n return True\n else:\n return False","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def assign_to_zone(pdot_x, pdot_y, pdot_size):\n x_size = 475; y_size=148\n zone_NW = [[0, int(y_size/2)], [0, int(x_size/3)]]\n zone_SW = [[int(y_size/2), y_size], [0, int(x_size/3)]]\n zone_N = [[0, int(y_size/2)], [int(x_size/3), int(x_size*2/3)]]\n zone_S = [[int(y_size/2), y_size], [int(x_size/3), int(x_size*2/3)]]\n zone_NE = [[0, int(y_size/2)], [int(x_size*2/3), x_size]]\n zone_SE = [[int(y_size/2), y_size], [int(x_size*2/3), x_size]]\n \n r = np.sqrt(pdot_size/np.math.pi)\n \n # pdot_size is considered as a surface\n\n # first assign to the main zone\n zones_and_percent = None\n \n # zone NW\n if in_interval(pdot_y, zone_NW[0]) and in_interval(pdot_x, zone_NW[1]):\n # y axis\n if in_interval(pdot_y + r, zone_NW[0]) and in_interval(pdot_y - r, zone_NW[0]):\n # x axis\n # Case 1 : the entire circle is in the zone\n if in_interval(pdot_x + r, zone_NW[1]) and in_interval(pdot_x - r, zone_NW[1]):\n zones_and_percent = [['NW_NO2', pdot_size]]\n # Case 2 : the circle is out of the zone - West side\n elif in_interval(pdot_x + r, zone_NW[1]) and not in_interval(pdot_x - r, zone_NW[1]):\n surface_diff = np.math.pi*(pdot_x - r)**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['NW_NO2', pdot_size_resize]]\n # Case 3 : the circle is out of the zone - East side\n elif not in_interval(pdot_x + r, zone_NW[1]) and in_interval(pdot_x - r, zone_NW[1]):\n surface_diff = np.math.pi*(int(x_size/3) - (pdot_x + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['NW_NO2', pdot_size_resize], ['N_NO2', surface_diff]]\n \n # Case 4 : the circle is out of the zone - North side\n elif in_interval(pdot_y + r, zone_NW[0]) and not in_interval(pdot_y - r, zone_NW[0]):\n surface_diff = np.math.pi*(pdot_y - r)**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['NW_NO2', pdot_size_resize]]\n # Case 5 : the circle is out of the zone - South side\n elif not in_interval(pdot_y + r, zone_NW[0]) and in_interval(pdot_y - r, zone_NW[0]):\n surface_diff = np.math.pi*(int(y_size/2) - (pdot_y + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['NW_NO2', pdot_size_resize], ['SW_NO2', surface_diff]]\n \n # zone SW\n elif in_interval(pdot_y, zone_SW[0]) and in_interval(pdot_x, zone_SW[1]):\n # y axis\n if in_interval(pdot_y + r, zone_SW[0]) and in_interval(pdot_y - r, zone_SW[0]):\n # x axis\n # Case 1 : the entire circle is in the zone\n if in_interval(pdot_x + r, zone_SW[1]) and in_interval(pdot_x - r, zone_SW[1]):\n zones_and_percent = [['SW_NO2', pdot_size]]\n # Case 2 : the circle is out of the zone - West side\n elif in_interval(pdot_x + r, zone_SW[1]) and not in_interval(pdot_x - r, zone_SW[1]):\n surface_diff = np.math.pi*(pdot_x - r)**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['SW_NO2', pdot_size_resize]]\n # Case 3 : the circle is out of the zone - East side\n elif not in_interval(pdot_x + r, zone_SW[1]) and in_interval(pdot_x - r, zone_SW[1]):\n surface_diff = np.math.pi*(int(x_size/3) - (pdot_x + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['SW_NO2', pdot_size_resize], ['S_NO2', surface_diff]]\n # Case 4 : the circle is out of the zone - North side\n elif in_interval(pdot_y + r, zone_SW[0]) and not in_interval(pdot_y - r, zone_SW[0]):\n surface_diff = np.math.pi*(int(y_size/2) - (pdot_y - r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['NW_NO2', surface_diff], ['SW_NO2', pdot_size_resize]]\n # Case 5 : the circle is out of the zone - South side\n elif not in_interval(pdot_y + r, zone_SW[0]) and in_interval(pdot_y - r, zone_SW[0]):\n surface_diff = np.math.pi*(y_size - (pdot_y + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['SW_NO2', pdot_size_resize]]\n\n # zone N\n elif in_interval(pdot_y, zone_N[0]) and in_interval(pdot_x, zone_N[1]):\n # y axis\n if in_interval(pdot_y + r, zone_N[0]) and in_interval(pdot_y - r, zone_N[0]):\n # x axis\n # Case 1 : the entire circle is in the zone\n if in_interval(pdot_x + r, zone_N[1]) and in_interval(pdot_x - r, zone_N[1]):\n zones_and_percent = [['N_NO2', pdot_size]]\n # Case 2 : the circle is out of the zone - West side\n elif in_interval(pdot_x + r, zone_N[1]) and not in_interval(pdot_x - r, zone_N[1]):\n surface_diff = np.math.pi*(int(x_size/3) - (pdot_x - r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['N_NO2', pdot_size_resize], ['NW_NO2', surface_diff]]\n # Case 3 : the circle is out of the zone - East side\n elif not in_interval(pdot_x + r, zone_N[1]) and in_interval(pdot_x - r, zone_N[1]):\n surface_diff = np.math.pi*(int(2*x_size/3) - (pdot_x + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['N_NO2', pdot_size_resize], ['NE_NO2', surface_diff]]\n # Case 4 : the circle is out of the zone - North side\n elif in_interval(pdot_y + r, zone_N[0]) and not in_interval(pdot_y - r, zone_N[0]):\n surface_diff = np.math.pi*(pdot_y - r)**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['N_NO2', pdot_size_resize]]\n # Case 5 : the circle is out of the zone - South side\n elif not in_interval(pdot_y + r, zone_N[0]) and in_interval(pdot_y - r, zone_N[0]):\n surface_diff = np.math.pi*(int(y_size/2) - (pdot_y + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['N_NO2', pdot_size_resize], ['S_NO2', surface_diff]]\n \n # zone S\n elif in_interval(pdot_y, zone_S[0]) and in_interval(pdot_x, zone_S[1]):\n # y axis\n if in_interval(pdot_y + r, zone_S[0]) and in_interval(pdot_y - r, zone_S[0]):\n # x axis\n # Case 1 : the entire circle is in the zone\n if in_interval(pdot_x + r, zone_S[1]) and in_interval(pdot_x - r, zone_S[1]):\n zones_and_percent = [['S_NO2', pdot_size]]\n # Case 2 : the circle is out of the zone - West side\n elif in_interval(pdot_x + r, zone_S[1]) and not in_interval(pdot_x - r, zone_S[1]):\n surface_diff = np.math.pi*(int(x_size/3) - (pdot_x - r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['S_NO2', pdot_size_resize], ['SW_NO2', surface_diff]]\n # Case 3 : the circle is out of the zone - East side\n elif not in_interval(pdot_x + r, zone_S[1]) and in_interval(pdot_x - r, zone_S[1]):\n surface_diff = np.math.pi*(int(2*x_size/3) - (pdot_x + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['S_NO2', pdot_size_resize], ['SE_NO2', surface_diff]]\n # Case 4 : the circle is out of the zone - North side\n elif in_interval(pdot_y + r, zone_S[0]) and not in_interval(pdot_y - r, zone_S[0]):\n surface_diff = np.math.pi*(int(y_size/2) - (pdot_y - r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['N_NO2', surface_diff], ['S_NO2', pdot_size_resize]]\n # Case 5 : the circle is out of the zone - South side\n elif not in_interval(pdot_y + r, zone_S[0]) and in_interval(pdot_y - r, zone_S[0]):\n surface_diff = np.math.pi*(y_size - (pdot_y + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['S_NO2', pdot_size_resize]]\n \n # zone NE\n elif in_interval(pdot_y, zone_NE[0]) and in_interval(pdot_x, zone_NE[1]):\n # y axis\n if in_interval(pdot_y + r, zone_NE[0]) and in_interval(pdot_y - r, zone_NE[0]):\n # x axis\n # Case 1 : the entire circle is in the zone\n if in_interval(pdot_x + r, zone_NE[1]) and in_interval(pdot_x - r, zone_NE[1]):\n zones_and_percent = [['NE_NO2', pdot_size]]\n # Case 2 : the circle is out of the zone - West side\n elif in_interval(pdot_x + r, zone_NE[1]) and not in_interval(pdot_x - r, zone_NE[1]):\n surface_diff = np.math.pi*(int(2*x_size/3) - (pdot_x - r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['NE_NO2', pdot_size_resize], ['N_NO2', surface_diff]]\n # Case 3 : the circle is out of the zone - East side\n elif not in_interval(pdot_x + r, zone_NE[1]) and in_interval(pdot_x - r, zone_NE[1]):\n surface_diff = np.math.pi*(x_size - (pdot_x + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['NE_NO2', pdot_size_resize]]\n # Case 4 : the circle is out of the zone - North side\n elif in_interval(pdot_y + r, zone_NE[0]) and not in_interval(pdot_y - r, zone_NE[0]):\n surface_diff = np.math.pi*(pdot_y - r)**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['N_NO2', pdot_size_resize]]\n # Case 5 : the circle is out of the zone - South side\n elif not in_interval(pdot_y + r, zone_NE[0]) and in_interval(pdot_y - r, zone_NE[0]):\n surface_diff = np.math.pi*(int(y_size/2) - (pdot_y + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['NE_NO2', pdot_size_resize], ['SE_NO2', surface_diff]]\n \n \n # zone SE\n elif in_interval(pdot_y, zone_SE[0]) and in_interval(pdot_x, zone_SE[1]):\n # y axis\n if in_interval(pdot_y + r, zone_SE[0]) and in_interval(pdot_y - r, zone_SE[0]):\n # x axis\n # Case 1 : the entire circle is in the zone\n if in_interval(pdot_x + r, zone_SE[1]) and in_interval(pdot_x - r, zone_SE[1]):\n zones_and_percent = [['SE_NO2', pdot_size]]\n # Case 2 : the circle is out of the zone - West side\n elif in_interval(pdot_x + r, zone_SE[1]) and not in_interval(pdot_x - r, zone_SE[1]):\n surface_diff = np.math.pi*(int(2*x_size/3) - (pdot_x - r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['SE_NO2', pdot_size_resize], ['S_NO2', surface_diff]]\n # Case 3 : the circle is out of the zone - East side\n elif not in_interval(pdot_x + r, zone_SE[1]) and in_interval(pdot_x - r, zone_SE[1]):\n surface_diff = np.math.pi*(x_size - (pdot_x + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['SE_NO2', pdot_size_resize]]\n # Case 4 : the circle is out of the zone - North side\n elif in_interval(pdot_y + r, zone_SE[0]) and not in_interval(pdot_y - r, zone_SE[0]):\n surface_diff = np.math.pi*(int(y_size/2) - (pdot_y - r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['NE_NO2', surface_diff], ['SE_NO2', pdot_size_resize]]\n # Case 5 : the circle is out of the zone - South side\n elif not in_interval(pdot_y + r, zone_SE[0]) and in_interval(pdot_y - r, zone_SE[0]):\n surface_diff = np.math.pi*(y_size - (pdot_y + r))**2\n pdot_size_resize = pdot_size - surface_diff\n zones_and_percent = [['SE_NO2', pdot_size_resize]]\n \n else:\n zones_and_percent = [['HS', pdot_size]]\n \n return zones_and_percent","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def simulate_emissions(EFs, sim_df, emission_sources, all_data):\n x_size = 475; y_size=148\n for dt, day in sim_df.iterrows():\n pollution_dots = []\n for indx, e in emission_sources.iterrows():\n pollution_size = EFs.get(e['primary_fuel'])*e['daily_gwh']\n pollution_location = (e['source_px'][0] + all_data.loc[dt]['DAILY_DISP_X_PX'],\n e['source_px'][1] + all_data.loc[dt]['DAILY_DISP_Y_PX'])\n pollution_dots.append([pollution_size,pollution_location])\n\n # Let's now assign each pollution dot to the corresponding sector of the island, or ignore \n # them if they are outside the image boundaries \n day_pollution = pd.Series([0.0,0.0,0.0,0.0,0.0,0.0],index=['NW_NO2', 'SW_NO2','N_NO2','S_NO2', 'NE_NO2','SE_NO2'])\n \n for pdot in pollution_dots: \n pdot_size = pdot[0]; pdot_x = pdot[1][0];pdot_y = pdot[1][1]\n\n assigned_zones = assign_to_zone(pdot_x, pdot_y, pdot_size)\n\n for assigned_zone in assigned_zones:\n if assigned_zone[0] != 'HS':\n day_pollution[assigned_zone[0]] += assigned_zone[1]\n\n sim_df.loc[dt] = day_pollution \n \n return sim_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# First Emission Factors Guess, every type roughly equal\nEFs = {'Oil':.3334,'Coal':0.3333,'Gas':0.3333}\nsim_df = simulate_emissions(EFs, sim_df, emission_sources, all_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sim_df.sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def evaluate_simulation(sim_df, all_data):\n Fit = 0.0\n for region in sim_df.columns:\n Fit += (sim_df[region].astype('float').corr(all_data[region].astype('float')))**2\n Fit = np.sqrt(Fit/6.0)\n return 1 - Fit","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"###############################\n# lower indicates a better fit\n###############################\nevaluate_simulation(sim_df, all_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# let's try giving a higher value to coal\nEFs = {'Oil':.40,'Coal':0.40,'Gas':0.20}\nsimulate_emissions(EFs, sim_df, emission_sources, all_data)\nevaluate_simulation(sim_df, all_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"EFs = {'Oil':.50,'Coal':0.4,'Gas':0.10}\nsimulate_emissions(EFs, sim_df, emission_sources, all_data)\nevaluate_simulation(sim_df, all_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def random_EFs():\n r_1 = random.randint(0,1000)\n r_2 = random.randint(0,1000)\n r_3 = random.randint(0,1000)\n r_sum = r_1+r_2+r_3\n return {'Oil':r_1/r_sum,'Coal':r_2/r_sum,'Gas':r_3/r_sum}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"####################################################################\n# Let's test a hundred random EFs to get an idea of a reasonable fit\n####################################################################\n\nbest_fit_EFs = {}\nbest_fit = 1.0\nfor i in tqdm(range(100)):\n EFs = random_EFs()\n simulate_emissions(EFs, sim_df, emission_sources, all_data)\n fit = evaluate_simulation(sim_df, all_data)\n if fit < best_fit:\n best_fit = fit\n best_fit_EFs = EFs.copy()\n #print(\"Best fit sor far:\", best_fit_EFs)\nprint(\"##############\")\nprint(\"Best fit was: \", best_fit_EFs)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### It seems clear that gas is relatively clean when compared to coal and oil, but still is responsible for some emissions.\n\n## For the final EF's we'll use bayesian optimisation to get more accurate results."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Finally we will just adjust the best fit emission factor to the proper units\nEF_oil = all_data['NO2'].sum()*best_fit_EFs['Oil']/generation['estimated_generation_gwh'].sum()*1000\n # Mult. by 1000 to convert to kWH\nEF_coal = all_data['NO2'].sum()*best_fit_EFs['Coal']/generation['estimated_generation_gwh'].sum()*1000 \nEF_gas = all_data['NO2'].sum()*best_fit_EFs['Gas']/generation['estimated_generation_gwh'].sum()*1000 ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display(Math(r'EF_{Oil} = ' + str(round(EF_oil,2)) + '\\;NO2\\;units\\;per\\;kWh'))\ndisplay(Math(r'EF_{Coal} = ' + str(round(EF_coal,2)) + '\\;NO2\\;units\\;per\\;kWh'))\ndisplay(Math(r'EF_{Gas} = ' + str(round(EF_gas,2)) + '\\;NO2\\;units\\;per\\;kWh'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Hyperopt - for a better search strategy using Bayesian optimization"},{"metadata":{"trusted":true},"cell_type":"code","source":"from hyperopt import hp, fmin, tpe, STATUS_OK, Trials\nfrom sklearn.metrics import mean_squared_error","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### HyperOpt instead of random testing"},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.stats.stats import pearsonr","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def acc_model_corr(EFs, all_data, emission_sources, sim_df):\n \n if 'params' in EFs:\n EFs = EFs['params']\n \n sim_df = simulate_emissions(EFs, sim_df, emission_sources, all_data)\n \n # Correlation\n NW_corr_sq = pearsonr(sim_df['NW_NO2'], all_data['NW_NO2'])[0]**2\n SW_corr_sq = pearsonr(sim_df['SW_NO2'], all_data['SW_NO2'])[0]**2\n N_corr_sq = pearsonr(sim_df['N_NO2'], all_data['N_NO2'])[0]**2\n S_corr_sq = pearsonr(sim_df['S_NO2'], all_data['S_NO2'])[0]**2\n NE_corr_sq = pearsonr(sim_df['NE_NO2'], all_data['NE_NO2'])[0]**2\n SE_corr_sq = pearsonr(sim_df['SE_NO2'], all_data['SE_NO2'])[0]**2\n \n metric = 1 - np.sqrt((NW_corr_sq + SW_corr_sq + N_corr_sq + S_corr_sq + NE_corr_sq + SE_corr_sq)/6)\n\n return metric","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def f(EFs):\n best = 1.0\n acc = acc_model_corr(EFs, all_data, emission_sources, sim_df)\n if acc < best:\n best = acc\n return {'loss': best, 'status': STATUS_OK}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def random_EFs():\n r_1 = random.randint(0,1000)\n r_2 = random.randint(0,1000)\n r_3 = random.randint(0,1000)\n r_sum = r_1+r_2+r_3\n return {'Oil':r_1/r_sum,'Coal':r_2/r_sum,'Gas':r_3/r_sum}\n\nEFs_list = [random_EFs() for y in range(1000)]\nparams = {'params': hp.choice('params', EFs_list)}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"trials = Trials()\nparams_best = fmin(f, params, algo=tpe.suggest, max_evals=100, trials=trials)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"best_fit_EFs = EFs_list[params_best['params']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"best_fit_EFs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Finally we will just adjust the best fit emission factor to the proper units\nEF_oil = all_data['NO2'].sum()*best_fit_EFs['Oil']/generation['estimated_generation_gwh'].sum()*1000\n # Mult. by 1000 to convert to kWH\nEF_coal = all_data['NO2'].sum()*best_fit_EFs['Coal']/generation['estimated_generation_gwh'].sum()*1000 \nEF_gas = all_data['NO2'].sum()*best_fit_EFs['Gas']/generation['estimated_generation_gwh'].sum()*1000 ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display(Math(r'EF_{Oil} = ' + str(round(EF_oil,2)) + '\\;NO2\\;units\\;per\\;kWh'))\ndisplay(Math(r'EF_{Coal} = ' + str(round(EF_coal,2)) + '\\;NO2\\;units\\;per\\;kWh'))\ndisplay(Math(r'EF_{Gas} = ' + str(round(EF_gas,2)) + '\\;NO2\\;units\\;per\\;kWh'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}