{"cells":[{"metadata":{"colab_type":"text","id":"ZqQ0WNGv608i"},"cell_type":"markdown","source":"# **DS4G: Modelling of emission of power plants **\n\n# Project overview and problem statement"},{"metadata":{},"cell_type":"markdown","source":"PROJECT OVERVIEW\nDevelop a methodology to calculate an average historical emissions factor of electricity generated for a sub-national region, using remote sensing data and techniques.\n\nThe Environmental Insights Explorer team at Google is keen to gather insights on ways to improve calculations of global emissions factors for sub-national regions. The ultimate goal of this challenge is to test if calculations of emissions factors using remote sensing techniques are possible and on par with calculations of emissions factors from current methodologies.\n\nPROBLEM STATEMENT\nCurrent emissions factors methodologies are based on time-consuming data collection and may include errors derived from a lack of access to granular datasets, inability to refresh data on a frequent basis, overly general modeling assumptions, and inaccurate reporting of emissions sources like fuel consumption. This begs the question: What if there was a different way to calculate or measure emissions factors? We’re challenging the Kaggle community to see if it’s possible to use remote sensing techniques to better model emissions factors. You will develop a methodology to calculate an average historical emissions factor for electricity generation in a sub-national region.\n\nWe’ve provided an initial list of datasets covering the geographic boundary of Puerto Rico to serve as the foundation for this analysis. As an island, there are fewer confounding factors from nearby areas. Puerto Rico also offers a unique fuel mix and distinctive energy system layout that should make it easier to isolate pollution attributable to power generation in the remote sensing data.\n\nParticipants will be tasked with developing a methodology to calculate an average annual historical emissions factor for the sub-national region. Participants will also be asked to provide an explanation of the conditions that would result in a higher/lower emissions factor, as well as a recommendation for how the methodology could be applied to calculate the emissions factor of electricity for another geospatial area using similar techniques. Bonus points will be awarded for smaller time slices of the average historical emissions factors, such as one per month for the 12-month period, and additional bonus points will be awarded for participants that develop methodologies for calculating marginal emissions factors for the sub-national region."},{"metadata":{},"cell_type":"markdown","source":"This notebook uses the notebook of Paul Mooney as a starting point for analyzing geographical information (see link below)\n\nhttps://www.kaggle.com/paultimothymooney/how-to-get-started-with-the-earth-engine-data\n\nThe notebook has the following contents:  \n1. Exploration of properties of power plants on Puerto Rico  \n2. Exploration of regional absorption aerosol index (aai)  \n3. Exploration of local absorption aerosol index (aai)\n4. Exploration of regional climate variables (temperature, humidity, wind speed)\n5. Calculation of local climate variables (temperature, humidity, wind speed)  \n6. Modelling the local aai (y) with features (X) from climate and from emission factors   \n7. Historical emission factor"},{"metadata":{"colab_type":"text","id":"6DZldYs2608k"},"cell_type":"markdown","source":"# 1: Properties of power plants on Puerto Rico #"},{"metadata":{"_kg_hide-input":true,"colab":{},"colab_type":"code","id":"weDvKvFh608l","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport os\nimport rasterio as rio\nimport folium\nimport tifffile as tiff\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    color={ 'Hydro' : 'lightblue', 'Solar' : 'orange', 'Oil' : 'darkblue', 'Coal' : 'black', 'Gas' : 'lightgray', 'Wind' : 'green' }\n    for i in range(0,len(df)):\n        popup = folium.Popup(str(df.plant[i]))   #[i:i+1]))\n        folium.Marker([df[latitude_column].iloc[i],df[longitude_column].iloc[i]],popup=popup, \n                      icon=folium.Icon(color=color[df.primary_fuel.iloc[i]])).add_to(plot)\n    return(plot)\n\ndef overlay_image_on_puerto_rico(file_name,band_layer):\n    band = rio.open(file_name).read(band_layer)\n    m = folium.Map([lat, lon], zoom_start=8, width=500, height=400)\n    folium.raster_layers.ImageOverlay(\n        image=band,\n        bounds = [[18.6,-67.3,],[17.9,-65.2]],\n        colormap=lambda x: (1, 0, 0, x),\n    ).add_to(m)\n    return m\n\ndef plot_scaled(file_name):\n    vmin, vmax = np.nanpercentile(file_name, (5,95))  # 5-95% stretch\n    img_plt = plt.imshow(file_name, cmap='gray', vmin=vmin, vmax=vmax)\n    plt.show()\n\ndef split_column_into_new_columns(dataframe,column_to_split,new_column_one,begin_column_one,end_column_one):\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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"power_plants = pd.read_csv('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\npower_plants = split_column_into_new_columns(power_plants,'.geo','latitude',50,66)\npower_plants = split_column_into_new_columns(power_plants,'.geo','longitude',31,48)\npower_plants['latitude'] = power_plants['latitude'].astype(float)\na = np.array(power_plants['latitude'].values.tolist()) # 18 instead of 8\npower_plants['latitude'] = np.where(a < 10, a+10, a).tolist() ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Below is a list of power plants on Puerto Rico in descending order of capacity.  \n\nThere are 34 power plants on the island. 16 plants use fossil fuel. 18 plants use a renewable source of energy. All fossil fuel powered plants are located near the shore. Of the power plants that use a renewable source of energy some are located inland and some are located near the shore.  "},{"metadata":{"_kg_hide-input":true,"colab":{},"colab_type":"code","id":"0eDVlNlC608r","outputId":"72f1995e-1183-4ad9-8659-02ed8806a31b","trusted":true},"cell_type":"code","source":"\npower_plants_df = power_plants.sort_values('capacity_mw',ascending=False).reset_index()\n#power_plants_df['prod_hrs_p_day']=power_plants_df.estimated_generation_gwh*1000/(power_plants_df.capacity_mw*365)\n\n\nbounds = [[18.6,-67.3,],[17.9,-65.2]]\n\npower_plants_df['img_idx_lt']=(((18.6-power_plants_df.latitude)*148/(18.6-17.9))).astype(int)\npower_plants_df['img_idx_lg']=((67.3+power_plants_df.longitude.astype(float))*475/(67.3-65.2)).astype(int)\npower_plants_df['plant']=power_plants_df.name.str[:3]+power_plants_df.name.str[-1]+'_'+power_plants_df.primary_fuel\n#power_plants_df['offset_img_aai']=(3-(power_plants_df.longitude.astype(float)+67)*2).astype(int)\n\npower_plants=power_plants_df[['name','latitude','longitude','primary_fuel','capacity_mw','img_idx_lt','img_idx_lg','plant']]\n\npower_plants","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"colab":{},"colab_type":"code","id":"K1DlwAqB608o","outputId":"91b580a9-17db-4f63-ec3a-784f3334aaf5","trusted":true},"cell_type":"code","source":"lat=18.200178; lon=-66.3 #-66.664513\nplot_points_on_map(power_plants,0,425,'latitude',lat,'longitude',lon,9)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print('Total green (solar, wind, hydro) energy capacity in MW :', power_plants_df.loc[((power_plants_df['primary_fuel']=='Hydro') | (power_plants_df['primary_fuel']=='Solar') | (power_plants_df['primary_fuel']=='Wind'))\n                    ,'capacity_mw'].sum())\n\nprint('Total gray (oil, gas, coal) energy capacity in MW :',power_plants_df.loc[((power_plants_df['primary_fuel']=='Coal') | (power_plants_df['primary_fuel']=='Oil') | (power_plants_df['primary_fuel']=='Gas'))\n                    ,'capacity_mw'].sum())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import matplotlib.patches as mpatches\n\nfig1 = plt.figure(figsize=(10, 5))\n\ncolor={ 'Hydro' : 'lightblue', 'Solar' : 'orange', 'Oil' : 'darkblue', 'Coal' : 'black', 'Gas' : 'lightgray', 'Wind' : 'green' }\nbarcolor=[]\nfor fuel in power_plants_df.primary_fuel : barcolor.append(color[fuel]) \n\nfig1 = fig1.add_subplot(111)\nfig1.bar(x=power_plants_df.index, height=power_plants_df.capacity_mw, width=0.6, color=barcolor)     \n    \nplt.yscale('log')\nplt.title('Power plants in Puerto Rico by primary fuel and in descending order of capacity ')\nplt.ylabel('Capacity (MW, log-scale)')\nplt.xlabel('Powerplants in Puerto Rico')\n\npatches=[]\nfor key, value in color.items(): patches.append(mpatches.Patch(color=value, label=key))\nfig1.legend(handles=patches)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The total capacity of 'gray' fossil fuel powered energy plants is 5760 MW.  \nThe total capacity of 'green' renewable source energy plants is 387 MW.  \n\nThe largest power unit uses oil as primary fuel and has a capacity over 1000 MW.  \nThe largest wind unit and the largest solar unit approximate the capacity of one of the smaller oil/gas units (approx. 50 to 100 MW).    \nThe capacity of all hydro units together approximate the capacity of one of the smaller oil/gas units (approx. 100 MW).  "},{"metadata":{},"cell_type":"markdown","source":"Additional information is taken from the internet to calculate the activity factor for the powerplants. \n\nInformation from eia.gov on electricity consumption of Puerto Rico gives a power consumption of 19.48 billion kWh (=19.480.000 MWh)for the year 2019.  \n\nInformation from index.mundi.com on the fuel consumption of power generation in Puerto Rico gives a fuel distribution oil/gas/coal/renewables of 40%/40%/18%/2%.\n\nThe calculation of the emission factor is based on the following (values used for oil) :  \nThe emission factor per day for the oil powerplants is calculated as 40% of total power consumption on Puerto Rico divided by the number of days in a year.  \nThe emission factor per oil powerplant is based on the fraction of the capacity of the plant vs the total capacity of the oil powerplants.  \n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# add information on capacity, type of fuel and activity factor\n\n# Information from eia.gov on electricity consumption of Puerto Rico gives a power consumption of 19.48 billion kWh (=19.480.000 MWh)for the year 2019\n# Information from index.mundi.com on the fuel consumption of power generation in Puerte Rico gives a distribution of 40%/40%/18%/2% for oil/gas/coal/renewables\n\nProd_day=int(19480000/365)\nprint('Average emission factor per day (production in MWh/day) : ',Prod_day)\n\n\n# With above information the drivers for the emission factor is calculated on a daily basis:\n\nEF_oil=19480000*0.4/365   # MWh/day\nEF_gas=19480000*0.4/365   # MWh/day\nEF_coal=19480000*0.18/365  # MWh/day\n\n# With the available capacity for oil, gas and coal plants the daily activity factor A is calculated\n\nprint('Emission factor (production in MWh/day) per day for oil: ',int(EF_oil),' gas: ',int(EF_gas),' and coal: ',int(EF_coal))\n\n#print(gray.groupby(by='primary_fuel').capacity_mw.sum())\n\nA_oil=EF_oil/power_plants_df.loc[power_plants_df.primary_fuel=='Oil','capacity_mw'].sum() \nA_gas=EF_gas/power_plants_df.loc[power_plants_df.primary_fuel=='Gas','capacity_mw'].sum() \nA_coal=EF_coal/power_plants_df.loc[power_plants_df.primary_fuel=='Coal','capacity_mw'].sum() \n\nprint('Activity factor of power plants (average hrs/day) for oil, gas and coal ',A_oil/24,' gas: ',A_gas/24,' and coal: ',A_coal/24)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The conclusion from the calculated activity factor is that the gas and coal plants are running at high activity rates (80% for gas and 88% for coal).   \nThe oil plants are running at a low activity factor (21%). This could also mean that some oil plants are stand-by and not running at all, where other oil plants are running at a higher activity factor.   "},{"metadata":{"colab_type":"text","id":"G1ca13u6608t"},"cell_type":"markdown","source":"# 2. Regional absorption aerosol index (aai)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# inspection of image information\nimage = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180708T172237_20180714T190743.tif'\nimg=rio.open(image)\n\n# print('Shape of array with data points :',tiff.imread(image).shape)\n# img.descriptions\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"from datetime import datetime\n\nfiles=[]\nfor dirname, _, filenames in os.walk('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2'):\n    for filename in filenames:\n        files.append(os.path.join(dirname, filename))\n\n# read all the absorbing aerosol index data into one list of arrays\naai_first_day=[]\naai_first_key=[]\naai_last_day=[]\naai_arr=[]\nband=5 #aerosol index\n#band=6 # cloud fraction\nfor i in range(0,len(files)):\n    aai_first_day.append(datetime.strptime(files[i][76:91], '%Y%m%dT%H%M%S').date())\n    aai_first_key.append(datetime.strptime(files[i][76:91], '%Y%m%dT%H%M%S').toordinal()+1) # correction of + 1 day in order to sync on climate data\n    aai_last_day.append(datetime.strptime(files[i][92:107], '%Y%m%dT%H%M%S').date())\n    #aai_first_day.append(pd.Timestamp(files[i][76:91]))\n    #aai_last_day.append(pd.Timestamp(files[i][92:107]))\n    aai_arr.append(rio.open(files[i]).read(band+1))\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"**Data cleaning** "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from pandas.plotting import register_matplotlib_converters\nregister_matplotlib_converters()\n\n\n# a=[]\n# for i,arr in enumerate(aai_arr): a.append(np.nanmean(arr))\na=[]\na_pos=[]\nnll=[]\nfor i in range(0,len(aai_arr)): \n    a.append(np.nanmean(aai_arr[i]))\n    a_pos.append(np.nanmean(np.clip(aai_arr[i],0,10000)))\n    nll.append(pd.isnull(aai_arr[i]).sum().sum())\n\n#aai_rgn=pd.DataFrame({'first': aai_first_day,'last':aai_last_day,'aai_rgn' : a, 'nll' : nll })\naai_rgn=pd.DataFrame({ 'first': aai_first_day,'last':aai_last_day,'aai_rgn' : a_pos, 'nll' : nll, 'aai_raw' : a,'key_date' : aai_first_key })\naai_rgn=aai_rgn.sort_values('first')\naai_rgn=aai_rgn.reset_index()\n\nfig1 = plt.figure(figsize=(20, 10))\nfig1.suptitle(\"data cleaning - mean of aai and number of nan per observation for 1) raw data (upper graphs) 2) cleaned data, #nan < 5% (middle graphs) 3) clipped data (lower graphs)\")\nax1 = fig1.add_subplot(321)\nax1.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,5], label='average aai per day - raw data', color='b')\nax1.legend()\nax2 = fig1.add_subplot(322)\nax2.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,4], label='# nan per observation', color='b')\nax2.legend()\n\naai_rgn=aai_rgn.loc[aai_rgn.nll <3515,:] # only select observations with # nan < 5%\n\nax3 = fig1.add_subplot(323)\nax3.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,5], label='average aai per day - cleaned for data with nan > 5%', color='b')\nax3.legend()\nax4 = fig1.add_subplot(324)\nax4.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,4], label='# nan per observation', color='b')\nax4.legend()\n\naai_rgn=aai_rgn.loc[aai_rgn.nll <3515,:] # only select observations with # nan < 5%\n\nax5 = fig1.add_subplot(325)\nax5.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='average aai per day - data <0 clipped to 0', color='b') \nax5.legend()\nax6 = fig1.add_subplot(326)\nax6.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,4], label='# nan per observation', color='b')\nax6.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The graphs to the left show the average aerosol absorption index (aai) for the region as a function of time.  \nThe graphs to the right show the number of nan (=not a number) -values that are in the data. Some daily data has a very high number of nan-values. This data cannot be used. Data with a nan-value percentage higher than 5% is discarded. \n\nFrom top to bottom the effect of the data cleaning is presented:  \nupper graphs: raw data   \nmiddle graphs: cleaned data (#nan < 5%)   \nlower graphs: clipped data (negative values clipped to zero) - to visualize the effect of an absorption index with only positive values, see explanation below  "},{"metadata":{},"cell_type":"markdown","source":"https://disc.gsfc.nasa.gov/information/glossary?title=Aerosol%20Index  \n\nText from above reference:  \n\"Aerosol Index. It is an index that detects the presence of uv-absorbing aerosols such as dust and soot. Positive values of Aerosol Index generally represent absorbing aerosols (dust and smoke) while small or negative values represent nonabsorbing aerosols and clouds.\"\n\nAccording to above reference only positive values of aai are an indication of absorbation by dust and smoke. An aai with only positive values is visualized in the lower graphs. In the modelling stage it may be interesting to investigate the relationship between emission factor and clipped aai-values instead of \"raw\" aai-values. "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# read only the absorbing aerosol index arrays with a nan-percentage <5% into one list of arrays for calculation of local aai data\nfiles=[]\nfor dirname, _, filenames in os.walk('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2'):\n    for filename in filenames:\n        files.append(os.path.join(dirname, filename))\n\naai_first_day=[]\naai_last_day=[]\naai_first_key=[]\naai_arr=[]\n#aai_arr_pos=[]\nband=5 #aerosol index\n#band=6 # cloud fraction\nfor i in range(0,len(files)):\n    a=rio.open(files[i]).read(band+1)\n    if pd.isnull(a).sum().sum() < 3515:\n        aai_first_day.append(datetime.strptime(files[i][76:91], '%Y%m%dT%H%M%S').date())\n        aai_first_key.append(datetime.strptime(files[i][76:91], '%Y%m%dT%H%M%S').toordinal()+1) # correction of + 1 day in order to sync on climate data\n        aai_last_day.append(datetime.strptime(files[i][92:107], '%Y%m%dT%H%M%S').date())        \n        aai_arr.append(np.clip(a,0,10000))  # clip negative values to zero\n#        aai_arr.append(a)  # raw values\n        \n# aai_arr is list of arrays with cleaned aai values (negative values clipped to zero, # nan in array < 5%)\n# aai_rgn is list of mean values of aai for the whole region (based on cleaned data)\n#print(len(aai_rgn))\n#print(len(aai_arr))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Inspection of measurements with highest regional averaged aai**"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# inspection of measurements with highest aai\n\naai_rgn.loc[aai_rgn.aai_rgn>.8,['key_date','first','last','aai_rgn']]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"colab":{},"colab_type":"code","id":"BLzBQe_N608u","outputId":"088d622c-ae5e-4c73-e5e6-9fb783f2e7e5","trusted":true},"cell_type":"code","source":"image1 = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180711T162527_20180718T185658.tif'\nimage2 = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180816T164847_20180822T182145.tif'\n\n# image1 : aai_rgn>0.8\n# image2 : aai_rgn approx 0.01\n\n#inspection of aai\n\nimage=image1\nimg=rio.open(image)\n\nband=5\nprint('inspection of image with high absorbing_aerosol_index (average aai for the region > 0.8)') #, img.descriptions[band])\n\nimage_band = rio.open(image).read(band+1)\n\n#plot_scaled(image_band)\n\nf2 = folium.Figure(width=500, height=400, title=img.descriptions[band])\nm = folium.Map([lat, lon], min_zoom=8, max_zoom=8, width='100%', height='100%').add_to(f2)\nfolium.raster_layers.ImageOverlay(\n    image=image_band,\n    bounds = [[18.6,-67.3,],[17.9,-65.2]],\n    colormap=lambda x: (1, 0, 0, x),\n).add_to(m)\nf2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3. Local absorption aerosol index (aai)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gray=power_plants.loc[((power_plants['primary_fuel']=='Coal') | (power_plants['primary_fuel']=='Oil') | (power_plants['primary_fuel']=='Gas')),\n                         ['name','primary_fuel','capacity_mw','img_idx_lt','img_idx_lg','plant']]\ngray.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# defining mask for locations with gray powerplants\nlocations=np.zeros((148,475))\n\n# setting the geographical area (+/- n pixels) around a plant location for averaging local data\n#n=10\nn=11 #maximum value to stay within bounds of image\n\nfor j in range(0,len(gray)):\n    locations[gray.iloc[j,3]-n:gray.iloc[j,3]+n,gray.iloc[j,4]-n:gray.iloc[j,4]+n]=np.ones((2*n,2*n))\n\n#plot_scaled(locations)\nprint('Overview of areas on Puerto Rico with power plants running on oil, gas or coal. These areas are selected for modelling aai vs emission factors')\n\nf1 = folium.Figure(width=500, height=400)\nm = folium.Map([lat, lon], min_zoom=8, max_zoom=8, width='100%', height='100%').add_to(f1)  #zoom_start=8\nfolium.raster_layers.ImageOverlay(\n    image=locations,\n    bounds = [[18.6,-67.3,],[17.9,-65.2]],\n    colormap=lambda x: (1, 0, 0, x),\n).add_to(m)\nf1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# aai value in proximity of all plants with all locations in location mask - proximity is +/- n points from location of plant\n\naai=[]\nfor j in range(0,len(gray)):\n    idx_lt=gray.iloc[j,3]\n    idx_lg=gray.iloc[j,4]\n    \n    aai_j=[]\n    for i in range(0,len(aai_arr)):\n        aai_j.append(np.nanmean(aai_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\n    \n    aai.append(aai_j)\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"\naa=pd.DataFrame({'key_date':np.array(aai_first_key), 'first': aai_first_day,'last':aai_last_day}) \n\nfor j in range(0,len(gray)):\n    aa[gray.iloc[j,5]]=aai[j]  #add average of aai for location of plant to dataframe with column name from df gray.plant\n\nprint('size of dataframe with aai data for gray-energy power-plant locations: ',aa.shape)\n\n# sorting dataframe on date to produce ordered time series\naa=aa.sort_values('key_date')\naa=aa.reset_index()\naa=aa.drop(columns=['index'])\naa=aa.fillna(0)\naa.head()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"\nfig3 = plt.figure(figsize=(20, 10))\nfig3.suptitle(\"absorbing aerosol index as function of time for different power plants\")\n#fig3 = fig3.add_subplot(111)\noffset=0\n#offset=8\n\nax=[]\nfor i in range(0,8):\n    ax.append(fig3.add_subplot(421+i))\n    ax[i].plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='mean_region', color='r')\n    ax[i].plot(aa.iloc[:,0], aa.iloc[:,3+i+offset], label=aa.columns[3+i+offset], color='b')\n    ax[i].set(ylim=(0, 2))                       #xlim=(-3, 3), ylim=(-3, 3))\n    ax[i].set_xlabel('time')\n    ax[i].set_ylabel('absorbing aerosol index')\n#    ax[i].set_title(\"aerosol index as function of time\")\n    ax[i].legend()\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig4 = plt.figure(figsize=(20, 10))\nfig4.suptitle(\"absorbing aerosol index as function of time for different power plants\")\n#fig3 = fig3.add_subplot(111)\n#offset=0\noffset=8\n\nax2=[]\nfor i in range(0,8):\n    ax2.append(fig4.add_subplot(421+i))\n    ax2[i].plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='mean_region', color='r')\n    ax2[i].plot(aa.iloc[:,1], aa.iloc[:,3+i+offset], label=aa.columns[3+i+offset], color='b')\n    ax2[i].set(ylim=(0, 2))                       #xlim=(-3, 3), ylim=(-3, 3))\n    ax2[i].set_xlabel('time')\n    ax2[i].set_ylabel('absorbing aerosol index')\n#    ax2[i].set_title(\"aerosol index as function of time\")\n    ax2[i].legend()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#simplified emissions-factor as the average aai of the plant location divided by the capacity of the plant  \n\naai_plant=aa.drop(columns=['key_date','first','last']).mean()\n\nprint('yearly average aai for the whole region : ',aai_rgn.aai_rgn.mean())\n\nfig5 = plt.figure(figsize=(20, 5))\nfig5.suptitle(\"yearly average emission per plant location\")\nax5 = fig5.add_subplot(111)\nax5.plot(aai_plant.index, aai_plant.values, label='average per plant', color='b')\nax5.plot(aai_plant.index, np.ones((len(aai_plant)))*aai_rgn.aai_rgn.mean(), label='average for the region', color='r')\nax5.legend()\n# simplified_emissions_factor = float(average_no2_emission/quantity_of_electricity_generated)\n# print('Simplified emissions factor (S.E.F.) for a single power plant on the island of Vieques =  \\n\\n', simplified_emissions_factor, 'S.E.F. units')","execution_count":null,"outputs":[]},{"metadata":{"colab_type":"text","id":"1cWq-KgX608w"},"cell_type":"markdown","source":"# 4. Regional climate variables (temperature, humidity, wind speed)"},{"metadata":{"_kg_hide-input":true,"colab":{},"colab_type":"code","id":"L7M5yI_9608x","outputId":"df4fb3d2-9cac-4b10-c77a-399576e2cf94","trusted":true},"cell_type":"code","source":"# image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gldas/gldas_20180702_1500.tif'\n# img1=rio.open(image)\n# image_band = rio.open(image).read(3)\n# print(img1.descriptions)\n# plot_scaled(image_band)\n\n# # for i in range(1,13):\n# #     image_band = rio.open(image).read(i)\n# #     print(img1.descriptions[i-1])\n# #     plot_scaled(image_band)\n\n# # image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/gfs_2018072118.tif'\n# # img2=rio.open(image)\n# # image_band = rio.open(image).read(3)\n# # print(img2.descriptions)\n# # plot_scaled(image_band)\n\n# #overlay_image_on_puerto_rico(image,band_layer=3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/gfs_2018072118.tif'\nimg2=rio.open(image)\nprint('Available information on climate factors')\n\nfor i in range(1,7):\n    image_band = rio.open(image).read(i)\n    print(img2.descriptions[i-1])\n    plot_scaled(image_band)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"files=[]\nfor dirname, _, filenames in os.walk('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs'):\n    for filename in filenames:\n        files.append(os.path.join(dirname, filename))\n\n# read all the climate data into lists of arrays\ngfs_day=[]\ngfs_key=[]\ntemp_arr=[]\nspec_hum_arr=[]\nrel_hum_arr=[]\nu_wind_arr=[]\nv_wind_arr=[]\npr_water_arr=[]\n#band=0 #temperature_2m_above_ground\n#band=6 # cloud fraction\nfor i in range(0,len(files)):\n    gfs_day.append(datetime.strptime(files[i][68:78], '%Y%m%d%H').date())\n    gfs_key.append(datetime.strptime(files[i][68:78], '%Y%m%d%H').toordinal())\n    temp_arr.append(rio.open(files[i]).read(1)) #temperature_2m_above_ground\n    spec_hum_arr.append(rio.open(files[i]).read(2)) #specific_humidity_2m_above_ground\n    rel_hum_arr.append(rio.open(files[i]).read(3)) # relative_humidity_2m_above_ground\n    u_wind_arr.append(rio.open(files[i]).read(4)) # u_component_of_wind_10m_above_ground\n    v_wind_arr.append(rio.open(files[i]).read(5)) # v_component_of_wind_10m_above_ground\n    pr_water_arr.append(rio.open(files[i]).read(6)) # precipitable_water_entire_atmosphere\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#gfs data is clean - geen nan in data!\n\nt=[] ; s=[] ; r=[] ; u=[] ; v=[] ; p=[]\n\nfor i in range(0,len(temp_arr)): \n    t.append(np.nanmean(temp_arr[i]))\n    s.append(np.nanmean(spec_hum_arr[i]))\n    r.append(np.nanmean(rel_hum_arr[i]))\n    u.append(np.nanmean(u_wind_arr[i]))\n    v.append(np.nanmean(v_wind_arr[i]))\n    p.append(np.nanmean(pr_water_arr[i]))\n\ngfs_rgn=pd.DataFrame({'day': gfs_day,'temp' : t, 'spec_hum' : s, 'rel_hum' : r, 'u_wind' : u, 'v_wind' : v, 'pr_water' : p, 'gfs_key': gfs_key })\ngfs_rgn=gfs_rgn.sort_values('day')\ngfs_rgn=gfs_rgn.reset_index()\n\nfig10 = plt.figure(figsize=(20, 10))\nfig10.suptitle(\"Average values for the region: temperature, specific_humidity, relative_humidity, u_comp_of_wind, v_comp_of_wind, precipitable_water (all left axis), aai on right axis\")\nax1 = fig10.add_subplot(321)\nax1.plot(gfs_rgn.iloc[:,1], gfs_rgn.iloc[:,2], label='average temperature', color='b')\nax12 = ax1.twinx()\nax12.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='average aai per day ', color='r') # aai_rgn.iloc[:,5] : waarde 3 geeft alleen positieve waardes van aai\nax1.legend() ; ax12.legend()\nax2 = fig10.add_subplot(322)\nax2.plot(gfs_rgn.iloc[:,1], gfs_rgn.iloc[:,3], label='average specific_humidity', color='b')\nax22 = ax2.twinx()\nax22.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='average aai per day ', color='r')\nax2.legend() ; ax22.legend()\nax3 = fig10.add_subplot(323)\nax3.plot(gfs_rgn.iloc[:,1], gfs_rgn.iloc[:,4], label='average relative_humidity', color='b')\nax32 = ax3.twinx()\nax32.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='average aai per day ', color='r')\nax3.legend() ; ax32.legend()\nax4 = fig10.add_subplot(324)\nax4.plot(gfs_rgn.iloc[:,1], gfs_rgn.iloc[:,5], label='average u_comp_of_wind', color='b')\nax42 = ax4.twinx()\nax42.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='average aai per day ', color='r')\nax4.legend() ; ax42.legend()\nax5 = fig10.add_subplot(325)\nax5.plot(gfs_rgn.iloc[:,1], gfs_rgn.iloc[:,6], label='average v_comp_of_wind', color='b')\nax52 = ax5.twinx()\nax52.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='average aai per day ', color='r')\nax5.legend() ; ax52.legend()\nax6 = fig10.add_subplot(326)\nax6.plot(gfs_rgn.iloc[:,1], gfs_rgn.iloc[:,7], label='average precipitable_water', color='b')\nax62 = ax6.twinx()\nax62.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='average aai per day ', color='r')\nax6.legend() ; ax62.legend()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# 5. Local climate variables (temperature, humidity, wind speed)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# climate values in proximity of all plants with all locations in location mask - proximity is +/- n points from location of plant\n\ntemp=[]\nspec_hum=[]\nrel_hum=[]\nu_wind=[]\nv_wind=[]\npr_water=[]\nfor j in range(0,len(gray)):\n    idx_lt=gray.iloc[j,3]\n    idx_lg=gray.iloc[j,4]\n    \n    temp_j=[] ; rel_hum_j=[] ; spec_hum_j=[] ; u_wind_j=[] ; v_wind_j=[] ; pr_water_j=[]\n    for i in range(0,len(temp_arr)):\n        temp_j.append(np.nanmean(temp_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of temp for location of plant\n        spec_hum_j.append(np.nanmean(spec_hum_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n]))\n        rel_hum_j.append(np.nanmean(rel_hum_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n]))\n        u_wind_j.append(np.nanmean(u_wind_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n]))\n        v_wind_j.append(np.nanmean(v_wind_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n]))\n        pr_water_j.append(np.nanmean(pr_water_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n]))\n    temp.append(temp_j)\n    spec_hum.append(spec_hum_j)\n    rel_hum.append(rel_hum_j)\n    u_wind.append(u_wind_j)\n    v_wind.append(v_wind_j)\n    pr_water.append(pr_water_j)\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gray=gray.iloc[:,:6]\n\n# weight of each powerplant as input to the emission model is the same. The model will calculate the relative weights of each plant.\ngray.loc[:,'EF_wght']=1\n\n# aggregation of climate data per plant location into one dataframe, addition of aai data per plant location \n# only use data for the dates that coincide for aai-data and for climate data\n\nww=pd.DataFrame({'key_date':gfs_key})\n\nXX=pd.DataFrame({})\nfor j in range(0,len(gray)):\n    #ww[gray.iloc[j,5]]=temp[j]  #add average of aai for location of plant to dataframe with column name from df gray.plant\n    ww['temp']=temp[j]\n    ww['spec_hum']=spec_hum[j]\n    ww['rel_hum']=rel_hum[j]\n    ww['u_wind']=u_wind[j]\n    ww['v_wind']=v_wind[j]\n    ww['pr_water']=pr_water[j]\n        \n    x=ww.groupby(by='key_date').agg(['max','min','mean'])\n    \n    X=pd.merge(aa.loc[:,['key_date',gray.iloc[j,5]]], x, how='inner', on='key_date')\n    X=X.rename(columns = {gray.iloc[j,5]:'aai'})\n    \n    c=gray.iloc[j,5]   #'EF_'+gray.iloc[j,5]\n    X[c]=np.ones((len(X)))*gray.iloc[j,6] # addition of EF_wght for each plant to the dataframe\n    \n    XX=pd.concat([XX,X], axis=0, sort=False) # aggregation of dataframe per plantlocation into one dataframe\n\nXX=XX.fillna(0) \nXX=XX.reset_index()\nXX    \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# i=13 #5 #13 #0-13\n# j=13 #5 #13 #0-13\n# offset=0\n# fig12 = plt.figure(figsize=(20, 5))\n# fig12.suptitle(\"Average values for the region: temperature, specific_humidity, relative_humidity, u_comp_of_wind, v_comp_of_wind, precipitable_water\")\n# ax1 = fig12.add_subplot(211)\n# ax1.plot(XX.iloc[i*324:(i+1)*324,1], XX.iloc[i*324:(i+1)*324,18], label=XX.columns[20], color='b')\n# ax12 = ax1.twinx()\n# ax12.plot(XX.iloc[i*324:(i+1)*324,1], XX.iloc[i*324:(i+1)*324,2], label=XX.columns[23+i], color='r') # aai_rgn.iloc[:,5] : waarde 3 geeft alleen positieve waardes van aai\n# ax1.legend() ; ax12.legend()\n# ax2 = fig12.add_subplot(212)\n# ax2.plot(XX.iloc[j*324:(j+1)*324,1], XX.iloc[j*324:(j+1)*324,18], label=XX.columns[20], color='b')\n# ax22 = ax2.twinx()\n# ax22.plot(XX.iloc[j*324+offset:(j+1)*324+offset,1], XX.iloc[j*324:(j+1)*324,2], label=XX.columns[23+i], color='r')\n# ax2.legend() ; ax22.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 6. Modelling aai (y) with features (X) from climate and from emission factors"},{"metadata":{},"cell_type":"markdown","source":"The general model for the actual emission is proposed as:\n\nE=EF x A x (1-R)\n\nE = actual emission (measured by aai)  \nEF = emission factor (production capacity and type of powerplant)  \nA = activity factor (fraction of the day that power plant is active or fraction of capacity that power plant is running at)  \nR = external factors that reduce the actual emission (climate variabeles)  \n\nIn the model that is calculated EF and A are combined into one variable (EF_powerplant for each plant) reflecting the emission factor and activity of the power plant.  \n\nAlso in the model the factor R is expanded to 6 climate variabeles (temperature, relative and specific humidity, horizontal and vertical wind-speed and precipitable water) that may have an effect on the actual emissions.  \n\nBoth factors are considered features (X) of the model in order to predict the aai (y) as a reference for the total emission.   \n\nThe feature importance as an output of the model will give an indication as to the importance of reduction factors versus the emission factors for predicting the total emission.  \n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"\ny=XX['aai']\n\nX=XX.drop(columns=['index','key_date','aai'])\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_squared_error\nimport xgboost as xgb\n\nmax_depth = 3\nmin_child_weight = 10\nsubsample = 0.5\ncolsample_bytree = 0.6\nobjective = \"reg:squarederror\" #'reg:linear',#\"reg:squarederror\"\nnum_estimators = 500 #2000 #1000 #3000  #200\nlearning_rate =  0.01 #0.01  #0.05 #0.003 # 0.3\n\nxgb_reg = xgb.XGBRegressor(max_depth=max_depth,\n            min_child_weight=min_child_weight,\n            subsample=subsample,\n            colsample_bytree=colsample_bytree,\n            objective=objective,\n            n_estimators=num_estimators,\n            learning_rate=learning_rate,\n            early_stopping_rounds=100,\n            num_boost_round = 2000)\n\nkf = KFold(n_splits=5, random_state=42, shuffle=True) # n_splits was 5\n\n\ni=0\ntestscore=[]\nfeature_imp=pd.DataFrame({'feature': X.columns})\nfor train_index, test_index in kf.split(X, y):\n#for train_index, test_index in gkf.split(X, y, groups):\n    X_train, X_test = X.iloc[train_index,:], X.iloc[test_index,:]\n    y_train, y_test = y[train_index], y[test_index]\n                      \n    xgb_reg.fit(X_train, np.ravel(y_train)) \n\n    y_pred = xgb_reg.predict(X_test)\n    test_score1 = mean_squared_error(y_test, y_pred)  \n    \n    testscore.append(test_score1)\n    i=i+1\n    feature_imp['importance'+str(i)]=xgb_reg.feature_importances_\n\nfeature_imp['mean']=feature_imp.iloc[:,1:i].mean(axis=1)\nfeature_imp['std']=feature_imp.iloc[:,1:i].std(axis=1)    \nprint('mean_squared_error on test_set:', testscore, np.mean(testscore))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig12 = plt.figure(figsize=(20, 5))\nfig12.suptitle(\"visual inspection of prediction of aai (y_pred) vs measured value of aai (y_test)\")\nax1 = fig12.add_subplot(111)\nax1.plot(range(0,len(y_test)), y_test, label='y_test', color='b')\n#ax12 = ax1.twinx()\nax1.plot(range(0,len(y_test)), y_pred, label='y_pred', color='r') # aai_rgn.iloc[:,5] : waarde 3 geeft alleen positieve waardes van aai\nax1.legend() #; ax12.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 7. Historical emission factor"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"feature_imp=feature_imp.sort_values('mean', ascending=False)\n#feature_imp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"feature_imp=feature_imp.sort_values('mean', ascending=True)\n\nplt.figure(figsize=(16, 12))\nplt.title(\"Feature importances in emission model of power plants in Puerto Rico\")\nplt.barh(range(X.shape[1]), feature_imp['mean'],\n       color=\"r\", xerr=feature_imp['std'], align=\"center\")\n# If you want to define your own labels,\n# change indices to a list of labels on the following line.\nplt.yticks(range(X.shape[1]), feature_imp['feature'])\nplt.ylim([-1, X.shape[1]])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Interpretation of feature importance**\n\nThe feature importance represents the weights in the calculated model that explain the relation between the features (climate factors and production factors) and the aerosol absorption index.\n\nThe feature importance for the production factors can be interpreted as the contribution of electricity production to the total emission for the specific location. E.g. the contribution to the total emission at Mayz_Gas is approx. two to three times higher than the contribution at Vega_Gas. In this case the higher emission is primarily explained by a higher production level (capacity of Mayz_Gas is 5 times higher than capacity of Vega_Gas).\n\nThe feature importance for the climate factors can be interpreted as the reduction factors in the modelling of emissions. \n\nThe feature importance shows that a large part (80% to 90%) of the variability in the model for predicting emission (=aai) is explained by climate factors and 10% to 20% of the variability is explained by production factors.\n\nThe standard deviation in feature importances for the production factors is slightly higher than the standard deviation in feature importances for the climate factors. It appears that the uncertainty in the production factors is higher, possible due to the uncertainty in production levels over time which are assumed constant but may vary in the actual situation. \n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gray=gray.iloc[:16,:7]\n\nprod_features=feature_imp.rename(columns= {'feature':'plant'})\ngray=pd.merge(gray, prod_features.loc[:,['plant','mean']], how='left', on='plant')\ngray=gray.rename(columns= {'mean':'emission_contrib'}) # emission contribution as calculated from the feature importances of the model\n#gray","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print('Contribution of production factors in the model to the measured total emissions (%) :', gray.emission_contrib.sum()*100)\n\ngray=gray.sort_values('emission_contrib', ascending=False)\n\n# calculation of maximum daily production (MWh) based on max. capacity of power plants\ngray['EF_max_MWh_day']=(gray['capacity_mw']*24).astype(int)\n\n# daily energy production for the region (MWh) distributed to plants according to emission distribution from the model\ngray['hist_emission_MWh_day']=(gray['emission_contrib']*Prod_day/gray.emission_contrib.sum()).astype(int)\n\n# activity factor calculated from model emission distribution and maximum daily production\ngray['Activity_%']=(gray['hist_emission_MWh_day']*100/gray['EF_max_MWh_day']).clip(upper=100).astype(int)\n\n# pollution factor: if activity > 100 (%) then emissions cannot be explained by power production because the plant is running beyond maximum capacity.\n# the pollution factor can explain these additional emissions attributable to choice of primary_fuel and generation of technology (age of plant).\ngray['Pollution_factor']=(gray['hist_emission_MWh_day']/gray['EF_max_MWh_day']).astype(int)\ngray","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The table is based on the model and gives an overview of powerplants in descending order of contribution to the total emission on Puerto Rico.  \n\nThe column 'emission_contrib' shows the weight of the feature importance as calculated by the model for each power plant.  \n\nThe column 'EF_max_MWh_day' gives the maximum production per day in MWh (capacity * 24 hrs) as a measure of maximum emission per day.  \n\nThe column 'hist_emission_MWh_day' gives the calculated production per day in MWh. For this calculation the yearly production of 19.480 GWh from eia.gov is divided by 365 days to reach an average daily production of 53.369 MWh for the region. The daily production is distributed to the locations by using the fraction of the emission_contrib for the location divided by the sum of the emission_contrib for all locations. \n\nThe column 'Activity (%)' is 'hist_emission_MWh_day' divided by 'EF_max_MWh_day'. The maximum is 100% since a plant cannot produce more power than its maximum capacity. \n\nWhen comparing the values in the columns 'hist_emission_MWh_day' and 'EF_max_MWh_day' one notices that values in the first column are sometimes higher, meaning that a higher production level is required to explain the emissions. Since Activity is already at 100% other mechanisms are needed to explain the additional emissions.  \n\nFor this explanation an additional pollution factor is calculated (column 'Pollution_factor') which can be thought of as a multiplier of emissions attributable to e.g. type of primary fuel or generation of power plant technology. Additional emissions may also be attributable to power generation if the data on power plant capacity is incorrect and too low. A check on the capacity specifications might be useful. Finally the additional emissions can also come from additional sources not part of this model.  \n\nWith the information in the table the distribution of used energy sources can be calculated. Earlier we saw a reported fuel usage distribution oil/gas/coal/renewables of 40%/40%/18%/2%. Let's make the calculation as a cross-check of the predictions of the model. \n\nBased on emissions and with ignoring capacity restrictions, column 'hist_emission_MWh_day' gives a distribution oil/gas/coal of 26%/67%/7%.  \n\nTaking into account capacity restrictions, a distribution oil/gas/coal of 23%/62%/15% is reached. Both distributions are somewhat in line with the reported fuel distribution, but there is still room for improvement. It seems that the model overstates the use of gas in the fuel usage distribution."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# fuel distribution purely based on historical emissions\nprint('fuel distribution based on historical emissions')\nprint(gray.groupby(by='primary_fuel').hist_emission_MWh_day.sum()/gray.hist_emission_MWh_day.sum())\nprint('  ')\n# fuel distribution based on historical emissions with capacity restrictions used where applicable. \ngray['est_prod']=gray['Activity_%']*gray['EF_max_MWh_day']\n\nprint('fuel distribution based on historical emissions with capacity restrictions')\nprint(gray.groupby(by='primary_fuel').est_prod.sum()/gray.est_prod.sum())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Investigations that were done but did not lead to additional information: \n* introduction of a reference point located in the middle of Puerto Rico, with no emission source but with the assumption that measured emission is partly attributable to all power plants in the region. The idea is that this might influence the distribution of the weights in feature importance compared to the model above. The calculated effects appeared to be small.  \n\n\nFurther investigations and improvements may include:\n* check on robustness and accuracy of prediction by using a different ML-model instead of XG-Boost\n* investigating why average regional emission is higher than average local emission for 14 out of 16 locations (see graph on \"yearly average emission per plant location\") \n* simplification in number of climate factors used\n\n* adding additional data from other sources\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"colab_type":"text","id":"-13q1iOg608z"},"cell_type":"markdown","source":"**Step 4: Connect to the Google Earth Engine API**"},{"metadata":{"_kg_hide-input":true,"colab":{},"colab_type":"code","id":"TFXRp_Hy6082","trusted":true},"cell_type":"code","source":"# from kaggle_secrets import UserSecretsClient\n# from google.oauth2.credentials import Credentials\n# import ee\n# import folium\n\n# def add_ee_layer(self, ee_image_object, vis_params, name):\n#   # https://github.com/google/earthengine-api/blob/master/python/examples/ipynb/ee-api-colab-setup.ipynb\n#   map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)\n#   folium.raster_layers.TileLayer(\n#     tiles = map_id_dict['tile_fetcher'].url_format,\n#     attr = 'Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n#     name = name,\n#     overlay = True,\n#     control = True\n#   ).add_to(self)\n\n# def plot_ee_data_on_map(dataset,column,begin_date,end_date,minimum_value,maximum_value,latitude,longitude,zoom):\n#     # https://github.com/google/earthengine-api/blob/master/python/examples/ipynb/ee-api-colab-setup.ipynb\n#     folium.Map.add_ee_layer = add_ee_layer\n#     vis_params = {\n#       'min': minimum_value,\n#       'max': maximum_value,\n#       'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5']}\n#     my_map = folium.Map(location=[latitude,longitude], zoom_start=zoom, height=500)\n#     s5p = ee.ImageCollection(dataset).filterDate(\n#         begin_date, end_date)\n#     my_map.add_ee_layer(s5p.first().select(column), vis_params, 'Color')\n#     my_map.add_child(folium.LayerControl())\n#     display(my_map)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#!cat ~/.config/earthengine/credentials","execution_count":null,"outputs":[]},{"metadata":{"colab":{},"colab_type":"code","id":"aqCw05AR6086","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# user_secret = \"earth_engine\" # Your user secret, defined in the add-on menu of the notebook editor\n# refresh_token = UserSecretsClient().get_secret(user_secret)\n# credentials = Credentials(\n#         None,\n#         refresh_token=refresh_token,\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)\n# ee.Initialize(credentials=credentials)","execution_count":null,"outputs":[]},{"metadata":{"colab_type":"text","id":"4Mrhh2896088"},"cell_type":"markdown","source":"Congratulations! You are connected!"},{"metadata":{"colab_type":"text","id":"5XUi-MLC6089"},"cell_type":"markdown","source":"** Step 5: Explore the emissions data using the Google Earth Engine API**\n* Plots will exist for 5 days or until they are regenerated."},{"metadata":{"_kg_hide-input":true,"colab":{},"colab_type":"code","id":"EILEV9yp608-","outputId":"7c157a29-ff4d-4824-85db-4e9648faa3ce","trusted":true},"cell_type":"code","source":"# dataset = \"COPERNICUS/S5P/NRTI/L3_NO2\"\n# column = 'absorbing_aerosol_index'\n# begin_date = '2018-07-08'\n# end_date = '2018-07-14'\n# minimum_value = 0.1 # 0.0000000001\n# maximum_value = 0.4 # 1\n# latitude = 60.17\n# longitude = 25.94\n# zoom = 5\n# plot_ee_data_on_map(dataset,column,begin_date,end_date,minimum_value,maximum_value,latitude,longitude,zoom)\n","execution_count":null,"outputs":[]},{"metadata":{"colab_type":"text","id":"YeLxgsiI609A"},"cell_type":"markdown","source":"** Step 6: Explore the weather data data using the Google Earth Engine API**\n* Plots will exist for 5 days or until they are regenerated."},{"metadata":{"_kg_hide-input":true,"colab":{},"colab_type":"code","id":"ddLBiG9D609B","outputId":"5c537320-b109-4120-aac5-c8d136f2b74a","trusted":true},"cell_type":"code","source":"# dataset = \"NOAA/GFS0P25\"\n# column = 'temperature_2m_above_ground'\n# begin_date = '2018-07-08'\n# end_date = '2018-07-14'\n# minimum_value = 0\n# maximum_value = 50\n# latitude = 18.20\n# longitude = -66.66\n# zoom = 8\n# plot_ee_data_on_map(dataset,column,begin_date,end_date,minimum_value,maximum_value,latitude,longitude,zoom)\n\n# dataset = \"NASA/GLDAS/V021/NOAH/G025/T3H\"\n# column = 'Tair_f_inst'\n# begin_date = '2018-07-08'\n# end_date = '2018-07-14'\n# minimum_value = 270\n# maximum_value = 310\n# latitude = 18.20\n# longitude = -66.66\n# zoom = 8\n# plot_ee_data_on_map(dataset,column,begin_date,end_date,minimum_value,maximum_value,latitude,longitude,zoom)","execution_count":null,"outputs":[]}],"metadata":{"colab":{"name":"how-to-get-started-with-the-earth-engine-data.ipynb","provenance":[]},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.6"}},"nbformat":4,"nbformat_minor":1}