{"cells":[{"metadata":{"colab_type":"text","id":"ZqQ0WNGv608i"},"cell_type":"markdown","source":"# **DS4G: Predicting emissions 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)\nhttps://www.kaggle.com/paultimothymooney/how-to-get-started-with-the-earth-engine-data\n\nThe suggestion of Rasyid Ridha is acknowledged and the y-axis is reverted in the processing of images to bring the images in line with images of Earth Engine\nhttps://www.kaggle.com/rasyidstat/fixing-the-reversed-y-axis\n\n\nThe notebook has the following contents:\n1. Summary\n2. Properties of power plants on Puerto Rico  \n3. Emission data on NO2 and aerosols  \n4. A prediction model for emission with local emission data\n5. Emission data for locations with and without a fossil fuel power plant\n6. Climate data on temperature, humidity and wind speed \n7. A prediction model for emission with local emission and weather data   \n"},{"metadata":{},"cell_type":"markdown","source":"# 1: Summary #\n\n\nThis notebook investigates the usability of emission data (from Sentinel 5P satellite) and of weather data (from GFS satellite) for the calculation of a historical emissions factor for power plants on Puerto Rico.\n\nA historical emissions factor is succesfully calculated for each fossil fuel powered plant on Puerto Rico. The calculation only needs a few data sources such as NO in the troposphere from Sentinel 5P, a list with locations and capacity of powerplants in the region and the yearly electrical power generation/consumption for the region. The calculation distributes the total electrical power for the region to the individual plants according to the geographical distribution of NO2 emissions. This distribution is further optimized considering constraints such as the maximum capacity of plants. With this distribution the historical emission level is determined. \n\nFor 2019 the generated power for Puerto Rico is reported as 19480 GWh. With this reference the model calculates the emission from power generation for the region as 37477 kg NO2/year.\n\nSeveral investigations were done to check if additional information would lead to better results. A calculation is made for all power plants, including reusable energy powered plants. It is expected that there is a difference in emissions between the locations of a fossil fuel plant and the location of a plant not running on fossil fuel. However this difference could not be observed. For Puerto Rico local NO2 emissions cannot distinguish between these types of locations.\n\nAlso it was investigated if a weekly pattern (weekday vs weekend) can be identified. Such a pattern is common for power plants in an industrial area. In the NO2 emissions data for Puerto Rico such a pattern is not observed. Possible explanations are the absence of industrial activities on Puerto Rico or the granularity of emissions data measurements where NO2 in the higher atmosphere and/or from other sources is blurring the weekly variations in the lower atmosphere which are expected from industrial activities such as power generation . \n\nFinally an investigation was done if data on the weather could improve the model described. Weather data is available on temperature, humidity, wind speed and precipitable water. The data for local weather and the data on local emissions are used to calculate a model that predicts local emissions. The relative importance of the features in the model is evaluated which of the features contribute to a better prediction. This investigation did not lead to improved results in calculation of the emission factor compared to the first more simple model. Therefor the more simple model is proposed to derive historical emissions for power plants in a region.\n\nThe proposed model is based on analysis for Puerto Rico. Further exploration for other types of areas needs to be investigated to see if the model is also valid for other types of environments. \n"},{"metadata":{"colab_type":"text","id":"6DZldYs2608k"},"cell_type":"markdown","source":"# 2: 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\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    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":"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":{"_kg_hide-input":true,"colab":{},"colab_type":"code","id":"0eDVlNlC608r","outputId":"72f1995e-1183-4ad9-8659-02ed8806a31b","trusted":true},"cell_type":"code","source":"power_plants_df = power_plants.sort_values('capacity_mw',ascending=False).reset_index()\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\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":{},"cell_type":"markdown","source":"Above 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 (see map). Of the power plants that use a renewable source of energy some are located inland and some are located near the shore. \n"},{"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 a simple 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\nThis emission factor is expressed in MWh/day production per plant. The idea behind it is that production is the driving force behind emission of the plant and it is assumed that there is a linear relationship between the two. A production figure can easily be compared to capacity constraints of a plant."},{"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) # MWh/day\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":"# 3. Emission data on NO2 and aerosols  #"},{"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)\nimg.descriptions\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The satellite images from Sentinel 5P contain information presented in the list above.  \n\nA first quick analysis to the relevance of this information is done for the first 7 bands. From this quick scan the following bands are analysed in further depth: \n* NO2_column_number_density   \n* tropospheric_NO2_column_number_density  \n* absorbing_aerosol_index  \n* cloud_fraction"},{"metadata":{"trusted":true,"_kg_hide-input":true},"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\n#date_sync=1     # MSE = 0.0187 - startwaarde en waarde ingeleverd in baseline versie    \n#date_sync=2 # MSE 0.0211\n#date_sync=3 # MSE 0.0213, slecht resultaat features (variatie > dan feature importance) \n#date_sync=4 # MSE 0.0211, lijkt niet zo kritisch voor uitkomst features ???\n#date_sync=0 # MSE = 0.01836, beste resultaat ook met features\n#date_sync=-2 # MSE = 0.01857, y_pred/y_test minder goed, features minder goed\ndate_sync=-1 # MSE = 0.01630, y_pred/y_test heel goed, features iets minder dan sync=0 - beste MSE dus default, date_sync=0 kan ook goed zijn\n\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=[]\n#band=0 # no2 column - redelijke locale informatie\nband=1 # no2 troposfeer - goede locale differentiatie zichtbaar\n#band=2 # no2 stratosfeer - geen locale informatie\n#band=3 # no2 slant column - redelijke locale informatie, maar minder dan band 1\n\n#band=5 #aerosol index - goede locale differentiatie zichtbaar\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()+date_sync) # 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(rio.open(files[i]).read(band+1))\n    \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\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 - regional mean of emission data and number of nan per observation for 1) raw data (upper graphs) 2) cleaned data, #nan < 5% (middle graphs) 3) #nan < 1% (lower graphs)\")\nax1 = fig1.add_subplot(321)\nax1.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,5], label='average NO2 in troposphere 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 NO2 in troposphere 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/5,:] # only select observations with # nan < 1%\n\nax5 = fig1.add_subplot(325)\n#ax5.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,3], label='average aai per day - data <0 clipped to 0', color='b') \nax5.plot(aai_rgn.iloc[:,1], aai_rgn.iloc[:,5], label='average NO2 in troposphere per day - cleaned for data with nan > 1%', 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 NO2 in troposphere 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: cleaned data (#nan < 1%)  "},{"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 emission 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\nef_day=[]\nef_dow=[]\nef_last_day=[]\nef_first_key=[]\nno2_col_arr=[] ; trop_arr=[] ; aai_arr=[] ; clfr_arr=[]\nno2_rgn=[] ; trop_rgn=[] ; aai_rgn=[] ; clfr_rgn=[]\n#band=0 # no2 column - redelijke locale informatie\n#band=1 # no2 troposfeer - goede locale differentiatie zichtbaar\n#band=2 # no2 stratosfeer - geen locale informatie\n#band=3 # no2 slant column - redelijke locale informatie, maar minder dan band 1\n#band=5 #aerosol index - goede locale differentiatie zichtbaar\n#band=6 # cloud fraction\n#print(band,img.descriptions[band])\nmax_nan=3515 # max 5% (3515) of data points nan\nfor i in range(0,len(files)):\n#     a=rio.open(files[i]).read(0+1) # no2 column\n#     b=rio.open(files[i]).read(1+1) # no2 troposfeer\n#     c=rio.open(files[i]).read(5+1) #aerosol index\n#     d=rio.open(files[i]).read(6+1) # cloud fraction\n    a=np.flip(rio.open(files[i]).read(0+1),0) # no2 column - image flipped to revert y-axis as suggested by https://www.kaggle.com/rasyidstat/fixing-the-reversed-y-axis\n    b=np.flip(rio.open(files[i]).read(1+1),0) # no2 troposfeer - see remark above\n    c=np.flip(rio.open(files[i]).read(5+1),0) #aerosol index - see remark above\n    d=np.flip(rio.open(files[i]).read(6+1),0) # cloud fraction - see remark above\n    if ((pd.isnull(a).sum().sum() < max_nan) & (pd.isnull(b).sum().sum() < max_nan) & (pd.isnull(c).sum().sum() < max_nan) & (pd.isnull(d).sum().sum() < max_nan)):\n        ef_day.append(datetime.strptime(files[i][76:91], '%Y%m%dT%H%M%S').date()) \n        ef_dow.append(datetime.strptime(files[i][76:91], '%Y%m%dT%H%M%S').isoweekday())\n        ef_first_key.append(datetime.strptime(files[i][76:91], '%Y%m%dT%H%M%S').toordinal()+date_sync) # correction of + 1 day in order to sync on climate data\n        ef_last_day.append(datetime.strptime(files[i][92:107], '%Y%m%dT%H%M%S').date())        \n        no2_col_arr.append(a)\n        trop_arr.append(b)\n        aai_arr.append(np.clip(c,0,10000))  # clip negative values to zero\n#        aai_arr.append(c)  # clip negative values to zero\n        clfr_arr.append(d)  \n        no2_rgn.append(np.nanmean(a))\n        trop_rgn.append(np.nanmean(b))\n        aai_rgn.append(np.nanmean(np.clip(c,0,10000)))        \n#        aai_rgn.append(np.nanmean(c))\n        clfr_rgn.append(np.nanmean(d))\n\n        \nef_rgn=pd.DataFrame({'day': ef_day,'no2' : no2_rgn, 'trop' : trop_rgn, 'aai' : aai_rgn, 'clfr' : clfr_rgn, 'ef_key': ef_first_key, 'dow' : ef_dow})\nef_rgn=ef_rgn.sort_values('day')\nef_rgn=ef_rgn.reset_index()        \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"** Visual inspection of measurements with emission**"},{"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'\nimage3 = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180708T172237_20180714T190743.tif'\n\nimage=image1\nimg=rio.open(image)\n\nprint('inspection of an image with some emissions  - orginal images') #, img.descriptions[band])\n\nimage_band = rio.open(image).read(band+1)\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":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"#with y-axis flipped\n\nimage=image1\nimg=rio.open(image)\n\nprint('inspection of an image with some emissions - images flipped along y-axis') #, img.descriptions[band])\n\nimage_band = np.flip(rio.open(image).read(band+1),0)\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":"# 4. A prediction model for emission with local emission data"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_locations_model(n):\n\n    # defining mask for locations with gray powerplants\n    locations=np.zeros((148,475))\n\n    for 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)\n    print('Overview of areas on Puerto Rico that are selected for modelling local emission factors')\n\n    f1 = folium.Figure(width=500, height=400)\n    m = folium.Map([lat, lon], min_zoom=8, max_zoom=8, width='100%', height='100%').add_to(f1)  #zoom_start=8\n    folium.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)\n    f1\n\n    return f1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_timeview_local_emission(df,mn,days):   # df is array with emission factor per location (in the columns) as function of time (in the rows)\n                                            # mn is list of mean for the region as function of time\n    fig3 = plt.figure(figsize=(20, 10))\n    fig3.suptitle(\"emission as function of time for different power plants\")\n\n    wkday=[]\n    for i,x in enumerate(ef_rgn.dow) : \n        if x>5: #<6: \n            wkday.append(ef_rgn.day[i]) \n\n    offset=0\n\n    ax=[]\n    for i in range(0,np.min((8,df.shape[1]-3))):\n        ax.append(fig3.add_subplot(421+i))\n        ax[i].plot(df.iloc[:days,1], mn[:days], label='mean_region', color='r')\n#        ax[i].plot(df.iloc[:days,1], wktrop[:days], label='weekdays', color='g')\n        ax[i].plot(df.iloc[:days,1], df.iloc[:days,3+i+offset], label=df.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('emission')\n    #    ax[i].set_title(\"aerosol index as function of time\")\n        ax[i].legend()\n        for x in wkday[:np.int(days*5/7)]:\n            ax[i].axvline(x, color='g',linewidth=4,alpha=0.15)\n    \n    if df.shape[1]-3 > 8 :\n        fig4 = plt.figure(figsize=(20, 10))\n        \n        offset=8\n\n        ax2=[]\n        for i in range(0,min((8,df.shape[1]-3-8-1))):\n            ax2.append(fig4.add_subplot(421+i))\n            ax2[i].plot(df.iloc[:,1], mn, label='mean_region', color='r')\n            ax2[i].plot(df.iloc[:,1], df.iloc[:,3+i+offset], label=df.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('emission')\n        #    ax2[i].set_title(\"aerosol index as function of time\")\n            ax2[i].legend()\n            for x in wkday[:np.int(days*5/7)]:\n                ax2[i].axvline(x, color='g',linewidth=4,alpha=0.15)\n    \n    return\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_locationsview_local_data(df1,rgn1,titel1,df2,rgn2,titel2 ):\n\n    #simplified emissions-factor as the average aai of the plant location divided by the capacity of the plant  \n\n    df_plant=df1.drop(columns=['key_date','first','dow']).mean()\n\n    #print('yearly average emission for the whole region : ',np.mean(rgn1))\n    labels_x=df_plant.index\n    x = np.arange(len(labels_x))  # the label location\n    width = 0.35 \n    fig5 = plt.figure(figsize=(20, 4))\n    fig5.suptitle(\"yearly average data per plant location :  \"+titel1+\" (left)  -   \"+titel2+\"  (right)\")\n    ax5 = fig5.add_subplot(121)\n    ax5.bar(x- width/2, df_plant.values, width, label='local average per cluster', color='b')\n    #ax5.bar(x+ width/2, np.ones((len(aai_plant)))*aai_rgn.aai_rgn.mean(), width, label='average for the region', color='r')\n    ax5.plot(x, np.ones((len(df_plant)))*np.mean(rgn1), label='average for the whole region', color='r')\n    ax5.plot(x, np.ones((len(df_plant)))*np.mean(df_plant.values), label='average for the locations', color='g')\n    ax5.set_xticks(x)\n    ax5.set_xticklabels(labels_x)\n    plt.setp(ax5.get_xticklabels(), rotation=45, ha=\"right\", rotation_mode=\"anchor\")\n    ax5.legend()\n    \n    df_plant=df2.drop(columns=['key_date','first','dow']).mean()\n\n    #print('yearly average emission for the whole region : ',np.mean(rgn2))\n    labels_x=df_plant.index\n    x = np.arange(len(labels_x))  # the label location\n    width = 0.35 \n    #fig5.suptitle(\"yearly average data per plant location - \"+titel2)\n    ax5 = fig5.add_subplot(122)\n    ax5.bar(x- width/2, df_plant.values, width, label='local average per cluster', color='b')\n    #ax5.bar(x+ width/2, np.ones((len(aai_plant)))*aai_rgn.aai_rgn.mean(), width, label='average for the region', color='r')\n    ax5.plot(x, np.ones((len(df_plant)))*np.mean(rgn2), label='average for the whole region', color='r')\n    ax5.plot(x, np.ones((len(df_plant)))*np.mean(df_plant.values), label='average for the locations', color='g')\n    ax5.set_xticks(x)\n    ax5.set_xticklabels(labels_x)\n    plt.setp(ax5.get_xticklabels(), rotation=45, ha=\"right\", rotation_mode=\"anchor\")\n    ax5.legend()\n        \n    return","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gray_plant=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']]\n#gray_plant=power_plants.loc[:,['name','primary_fuel','capacity_mw','img_idx_lt','img_idx_lg','plant']]\n\n\ngray_plant=gray_plant.sort_values('primary_fuel')\n\n# calculation of maximum daily production (MWh) based on max. capacity of power plants\ngray_plant['Max_Prod']=(gray_plant['capacity_mw']*24).astype(int)\n\nprint('Dataframe for local calculation of emission factors : ')\n#print(gray_plant.head())\n\nn=11 # size of area around a plant where emission is averaged for further calculation\n#n=4\n#n=3\n\ngray=pd.DataFrame({})\nfor j in range(0,len(gray_plant)):\n    idx_lt=gray_plant.iloc[j,3]\n    idx_lg=gray_plant.iloc[j,4]\n    temp=gray_plant.loc[(((gray_plant.img_idx_lt>idx_lt-n) & (gray_plant.img_idx_lt<idx_lt+n)) & \n                         ((gray_plant.img_idx_lg>idx_lg-n) & (gray_plant.img_idx_lg<idx_lg+n))),gray_plant.columns[:7]]\n    if len(temp)>1 : # if there are more plants in 1 location, group them together into 1 source\n        temp1=pd.DataFrame([[str(len(temp))+' locations','comb',np.int(temp.capacity_mw.sum()),np.int(temp.img_idx_lt.mean()),np.int(temp.img_idx_lg.mean()),\n                           temp.plant.sum(),temp.Max_Prod.sum()]], columns=temp.columns)\n        gray_plant.loc[(((gray_plant.img_idx_lt>idx_lt-n) & (gray_plant.img_idx_lt<idx_lt+n)) & \n                         ((gray_plant.img_idx_lg>idx_lg-n) & (gray_plant.img_idx_lg<idx_lg+n))),'cluster']=temp.plant.sum()\n        gray=pd.concat([gray,temp1])\n    else :    # if there is 1 plant in a location, keep this one as the source\n        gray=pd.concat([gray,temp])\n        gray_plant.loc[(((gray_plant.img_idx_lt>idx_lt-n) & (gray_plant.img_idx_lt<idx_lt+n)) & \n                         ((gray_plant.img_idx_lg>idx_lg-n) & (gray_plant.img_idx_lg<idx_lg+n))),'cluster']=temp.plant.sum()\n        \ngray=gray.drop_duplicates()\n\ngray=gray.rename(columns = { 'plant': 'cluster', 'Max_Prod': 'Max_Prod_clst'})\n\nprint(gray)\n\nplot_locations_model(n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"See above a picture of the locations of fossil fuel powered plants that are used in further calculations. Overlapping areas of power plants that are in close proximity of eachother are merged into one area (cluster) with a production capacity equal to the sum of the capacity of the individual plants.  \n\nThe table above gives the clusters that are generated for the calculation. Max_Prod_clst is the maximum daily production of the cluster in MWh/day."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# emission value in proximity of all plants with all locations in location mask - proximity is +/- n points from location of plant\nno2=[]\ntrop=[]\naai=[]\nclfr=[]\nfor j in range(0,len(gray)):\n    idx_lt=gray.iloc[j,3]\n    idx_lg=gray.iloc[j,4]\n        \n    no2_j=[]\n    trop_j=[]\n    aai_j=[]\n    clfr_j=[]\n    for i in range(0,len(trop_arr)):\n        no2_j.append(np.nanmean(no2_col_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\n        trop_j.append(np.nanmean(trop_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\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        clfr_j.append(np.nanmean(clfr_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\n        \n    no2.append(no2_j)\n    trop.append(trop_j)\n    aai.append(aai_j)\n    clfr.append(clfr_j)\n    \n    \n# list to DataFrame with dates   \nno=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow})\ntro=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow}) \naa=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow})\nclf=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow})\n\nfor j in range(0,len(gray)): #add average of emission for each location (plant) to dataframe with in the column name taken from df gray.cluster\n    no[gray.iloc[j,5]]=no2[j]\n    tro[gray.iloc[j,5]]=trop[j]\n    aa[gray.iloc[j,5]]=aai[j]  #add average of aai for each location (plant) to dataframe with in the column name taken from df gray.cluster\n    clf[gray.iloc[j,5]]=clfr[j]\n    \n#print('size of dataframe with emission data for gray-energy power-plant locations: ',no.shape, tro.shape, aa.shape, clf.shape)\n\n# sorting dataframe on date to produce ordered time series\nno=no.sort_values('key_date') ; tro=tro.sort_values('key_date') ; aa=aa.sort_values('key_date') ; clf=clf.sort_values('key_date')\n\nno=no.reset_index() ; tro=tro.reset_index() ; aa=aa.reset_index() ; clf=clf.reset_index()\nno=no.drop(columns=['index']) ; tro=tro.drop(columns=['index']) ; aa=aa.drop(columns=['index']) ; clf=clf.drop(columns=['index'])\nno=no.fillna(0) ; tro=tro.fillna(0) ; aa=aa.fillna(0) ; clf=clf.fillna(0)\n    \n#emission_correlation() # show locations with high correlation\n\n#plot_timeview_local_emission() # show timeview of emission per location\n\nplot_locationsview_local_data(no,no2_rgn,'NO2 column',tro,trop_rgn, 'NO2 in troposphere') # show emissions per location vs average for the region\nplot_locationsview_local_data(aa,aai_rgn, 'Aerosol absorption index (aai)',clf,clfr_rgn, 'Cloud fraction') # show emissions per location vs average for the region\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Above 4 graphs that present the local measurements of emission data averaged over the year. A first glance shows that NO2 column and NO2 in troposphere have a similar pattern but NO2 in troposphere has a better distinction between locations. In the lower graphs the high values for Mayz_Gas are curious since Mayz_Gas does not have the highest production capacity.\n\nFor aerosol absorption index some pre-processing of the data is done. According to below reference only positive values of aerosol absorption index are an indication of absorption by dust and smoke. The pre-processing is done on the negative values in the data which are clipped to zero. \n\nhttps://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.\""},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot_timeview_local_emission(tro,ef_rgn.trop,324) # show timeview of emission per location\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The graphs shows the measured emission for each location over the year. The green background in the graphs represent the weekends (saturday and sunday). The idea behind this is that for man-made emissions a weekly pattern can be observed (emissions higher on weekdays and lower in the weekend). From visual inspection it is unclear if such a pattern exists in the emission data.  "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"\ndata_dow=ef_rgn.pivot(columns='dow',values='trop')\n\ndays=['mon','tue', 'wed','thu','fri','sat','sun']\n\nf11, ax11 = plt.subplots(1, 2, figsize=(20, 4))\nchart=sns.boxplot(data=data_dow, ax=ax11[0])\nchart.set_xticklabels(days, rotation=45, horizontalalignment='right')\nchart.set_title('Regional NO2 in troposphere')\nchart.set_ylabel('emission (NO2 troposphere)') \n#data_dow\n\n\np=4 # number selects location to show\ndata_dow=tro.pivot(columns='dow',values=tro.columns[p])\n\nchart=sns.boxplot(data=data_dow, ax=ax11[1])\nchart.set_xticklabels(days, rotation=45, horizontalalignment='right')\nchart.set_title(tro.columns[p])\nchart.set_ylabel('emission (NO2 troposphere)') ;\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The box-plots show the measured emission per weekday for the region (left) and for one location (right).  \nNo difference in emissions between weekdays can be observed. A weekly pattern is common for power plants in industrialized areas. Since it cannot be observed in emissions on Puerto Rico either:  \n* there is no weekly pattern because there is little industry   \n* the weekly pattern is blurred because of background emissions  \n* aggregation of NO2 in the higher atmosphere conceils weekly variations in the lower parts.     "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# calculation of monthly emission factors\nimport seaborn as sns\n#sns.set(style=\"ticks\", palette=\"pastel\")\n\n\ndef monthly_data_location(df1,df2,titel1,titel2) :\n\n    maand=df1\n    t=[]\n    for i, x in enumerate (maand.iloc[:,0]) : t.append(datetime.fromordinal(x+1).month)\n    maand['month']=t\n    maand=maand.drop(columns=['key_date','first','dow'])\n    maand=maand.groupby(by='month').agg('mean')\n\n    emis_maand=(maand.T*Prod_day/maand.T.sum(axis=0)).astype(int)\n\n    f1, ax1 = plt.subplots(1, 2, figsize=(20, 4))\n    chart=sns.boxplot(data=emis_maand.T, ax=ax1[0])\n    chart.set_xticklabels(chart.get_xticklabels(), rotation=45, horizontalalignment='right')\n    chart.set_title(titel1)\n    chart.set_ylabel('production (MWh/day)') \n    #chart.set_xlabel('location')\n        \n    maand=df2\n    t=[]\n    for i, x in enumerate (maand.iloc[:,0]) : t.append(datetime.fromordinal(x+1).month) # correction +1 to sync first date on 1st day of July\n    maand['month']=t\n    maand=maand.drop(columns=['key_date','first','dow'])\n    maand=maand.groupby(by='month').agg('mean')\n\n    emis_maand=(maand.T*Prod_day/maand.T.sum(axis=0)).astype(int)\n\n    chart=sns.boxplot(data=emis_maand.T, ax=ax1[1])\n    chart.set_xticklabels(chart.get_xticklabels(), rotation=45, horizontalalignment='right')    \n    chart.set_title(titel2)\n    chart.set_ylabel('production (MWh/day)') \n    #chart.set_xlabel('location')\n    return\n\nmonthly_data_location(no,tro,'monthly data NO2 column','monthly data troposphere' )\nmonthly_data_location(aa,clf,'monthly data aerosol absorption index','monthly data cloud fraction' )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The box-plots show the variability in measured emission per location per month. Emission is calculated as a corresponding production figure in order to compare it with available production capacity (see graph below)   \n\nThe upper plots give the best distinction between locations. Variation of emissions levels over the months is more or less the same for all locations. From this it is concluded that there is little monthly variation in emissions between locations. An emission factor calculated for the year or for the month will not differ very much in this case. If larger variations in monthly emissions per location can be observed that would be an indication that emission factors per location change over time.   \n\nSince the emission factor NO2 in troposphere gives the best distcintion between location this emission measurement is used to determine the emission factor in further calculations. "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig6 = plt.figure(figsize=(20, 4))\n#fig5.suptitle(\"yearly average data per plant location :  \"+titel1+\" (left)  -   \"+titel2+\"  (right)\")\nax6 = fig6.add_subplot(121)\nlabels_x=gray.cluster\nx = np.arange(len(labels_x))  # the label location\nwidth = 0.35 \nax6.bar(x-width/2,gray.Max_Prod_clst.values, width, label='Max_Prod (daily production capacity for the cluster)', color='b')\n#ax5.bar(x+ width/2, np.ones((len(aai_plant)))*aai_rgn.aai_rgn.mean(), width, label='average for the region', color='r')\nax6.plot(np.ones((len(gray)))*np.mean(gray.Max_Prod_clst), label='average daily capacity per location for the region', color='r')\nax6.plot(np.ones((len(gray)))*Prod_day/len(gray), label='daily production level for the region averaged per location', color='g')\nax6.set_xticks(x)\nax6.set_xticklabels(labels_x)\nax6.set_ylabel('Production (MWh/day)')\nplt.setp(ax6.get_xticklabels(), rotation=45, ha=\"right\", rotation_mode=\"anchor\")\nax6.legend();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A graph which shows the available production capacity (= driver of emissions) per cluster.\n\nThe table below gives the measured emission (EF_contrib) per location. The total daily production of electricity is distributed to the locations by the weight of the emission at the location divided by the sum of emissions for all locations. This leads to a daily production per location (hist_prd in MWh per day) and an activity factor (Activity %) (see next table).    "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gray=gray.iloc[:,:7]\n\nlocal_plant=tro.drop(columns=['key_date','first','dow']).mean()\n\navg_location=pd.DataFrame(local_plant, columns=['EF_contrib']).reset_index().rename(columns = {'index': 'cluster'})\n\ngray=pd.merge(gray, avg_location.loc[:,['cluster','EF_contrib']], how='left', on='cluster')\n\ngray=gray.drop(columns=['img_idx_lt','img_idx_lg'])\n\ngray=gray.fillna(0)\ngray\n\n","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.EF_contrib.sum()*100)\n\ngray=gray.sort_values('EF_contrib', ascending=False)\n\n# # calculation of maximum daily production (MWh) based on max. capacity of power plants\n# gray['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_prod_day']=(gray['EF_contrib']*Prod_day/gray.EF_contrib.sum()).astype(int)\n\n# activity factor calculated from model emission distribution and maximum daily production\ngray['Activity_%']=(gray['hist_prod_day']*100/gray['Max_Prod_clst']).clip(upper=100).astype(int)\n\n# multiplier: if activity > 100 (%) then emissions cannot be explained by power production because the plant is running beyond maximum capacity.\n# the multiplier can explain these additional emissions attributable to choice of primary_fuel and generation of technology (age of plant).\ngray['multiplier']=(gray['hist_prod_day']/gray['Max_Prod_clst']).astype(int)\n\ngray['EF_MaxCap']=gray['EF_contrib']\ngray.loc[gray.hist_prod_day>gray.Max_Prod_clst,'EF_MaxCap']=0\ngray['hist2_prod_day']=(gray['EF_MaxCap']*(Prod_day-gray.loc[gray.hist_prod_day>gray.Max_Prod_clst,'Max_Prod_clst'].sum())/gray.EF_MaxCap.sum()).astype(int)\ngray.loc[gray.hist_prod_day>gray.Max_Prod_clst,'hist2_prod_day']=gray.Max_Prod_clst\ngray['Activity2_%']=(gray['hist2_prod_day']*100/gray['Max_Prod_clst']).astype(int)     #.clip(upper=100).astype(int)\n\ngray['EF_hist2']=gray['EF_MaxCap']\ngray.loc[gray.EF_MaxCap==0,'EF_hist2']=gray['EF_contrib']*gray['hist2_prod_day']/gray['hist_prod_day']\n\ngray_cluster=gray\ngray_cluster","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The column 'hist_prod_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_prod_day' divided by 'Max_prod_clst'. The maximum is 100% since a plant cannot produce more power than its maximum capacity.\n\nWhen comparing the values in the columns 'hist_prod_day' and 'Max_prod_clst' 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 multiplier is calculated (column 'multiplier') which can be thought of as a multiplier to the emission from the power plant running at maximum capacity to reach the measured level. The assumption is that the measured level can only in part be explained by power generation since it exceeds the maximum emission capability of the plant.\n\nTo get a better distribution of power generation across plants, with the constraint of maximum production capacity per plant, a second distribution is calculated ('hist2_prod_day'). In this calculation the production of power plants with a multiplier is set at the maximum daily production. At the same time the contribution of these plants to the emission is neglected (EF_max_cap set to zero) since it is lower then can be explained by its capacity. The production for the other plants is calculated according to the remaining weights in the measured emission.\n\nThis leads to a production distribution in column 'hist2_prod_day' and an associated activity factor Activity2_%, not capped to 100% to check if the calculated plant production level remains within the available capacity.\n\nFinally a historical emission factor is calculated 'EF_hist2' based on historical measurements and scaled back from the measurement if the maximum production of the plant cannot explain (is lower then) the measured emission level.  \n\nEF-hist2 has the same dimension as the measurement NO2 in troposphere (mol/m^2) from which an absolute emission value needs to be derived. \n\nA different approach to reach an absolute emission value is based on NO2 emission factors known for electricity generation from coal, oil and gas e.g from http://data.ec.gc.ca/data/substances/monitor/canada-s-official-greenhouse-gas-inventory/E-Tables-Electricity-Canada-Provinces-Territories/?lang=en referenced on EIE page https://insights.sustainability.google/methodology#some-notes   \n\nThe reference gives values for coal: 0.01 gr NO2/kWh (Tab A13-11), oil: 0.001 gr NO2/kWh (Tab A13-2) and natural gas: 0.0007 gr NO2/kWh (Tab A13-11). With columns 'primary_fuel' and 'hist_prod_day' in the table below the absolute NO2 emission per plant can be calculated.\n\n*It is interesting to note that the calculated absolute emission for plant 'A.E._Coal' is 10x the calculated absolute emission for plant 'Palo_Oil'. Both run at approx. the same production level (6000 MWh/day +/- 10%). The emission factor from the troposphere measurements give values that are more or less the same: 6.5*e-6 vs 7.2*e-6. The higher absolute emission level of a coal plant as reported in the reference is not observed in the local emission measurements.*"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# post-processing - for each plant\ngray_plant=gray_plant.iloc[:,:8]\n#gray=gray.rename(columns = {'plant':'cluster'})\n\ngray_plant=pd.merge(gray_plant, gray[['hist2_prod_day','Max_Prod_clst','cluster','EF_hist2']], how='left', on='cluster')\n\ngray_plant['hist_prod_day']=(gray_plant['Max_Prod']*gray_plant['hist2_prod_day']/gray_plant['Max_Prod_clst']).astype(int)\n\ngray_plant['Activity_%']=(gray_plant['hist_prod_day']*100/gray_plant['Max_Prod']).astype(int)  \n\ngray_plant['EF_hist']=(gray_plant['Max_Prod']*gray_plant['EF_hist2']/gray_plant['Max_Prod_clst'])\n\nNO_factors={ 'Coal': 0.01 , 'Oil': 0.001, 'Gas': 0.0007 } # gr/kWh is same as kg/MWh\n\nfor i,x in enumerate(gray_plant.primary_fuel) : gray_plant.loc[gray_plant.primary_fuel==x,'EF_abs']=NO_factors[x]*gray_plant.loc[gray_plant.primary_fuel==x,'hist_prod_day']\n\ngray_plant=gray_plant.drop(columns = ['hist2_prod_day', 'Max_Prod_clst','EF_hist2'])\n  \ngray_location=gray_plant\n\ngray_location","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# fuel distribution purely based on historical emissions\nprint('fuel distribution based on historical emissions')\nprint(gray_plant.groupby(by='primary_fuel').hist_prod_day.sum()/gray_plant.hist_prod_day.sum())\nprint('  ')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The baseline for emission of NO2 from electricity generation in Puerto Rico is based on the daily power generation * the emission factor for coal/oil/gas."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('Baseline for emission of NO2 from electricity generation in Puerto Rico is ',int(gray_location.EF_abs.sum()*365),'kg NO2/year' )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To set this value into perspective a simple approximation of emission is calculated as the product of electricity generation in the year times the fuel distribution times the emission per kWh for the fuel. "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('Simple approximation of emission of NO2 from electricity generation in Puerto Rico : ',int(19480000*(.2*0.01+.4*0.001+.4*0.0007)),'kg NO2/year' )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This value is higher then the value from the model due to a different fuel distributon calculated from the model.  \n\nEstimation of emission of NO2 from satellite measurements is based on the average NO2 density (mol/m^2) on plant locations. Molecular weight of NO2 is 46,0055 g/mol.   \nFurther we need to estimate the area for which the emission is generated and the lifetime of the generated emissions. For simplicity we use 1 square kilometer for the area and 1 day for the lifetime of the emissions"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print('Estimated NO2 emission based on measurements at plant locations ', int(tro.iloc[:,3:13].mean().mean()*365*10*1000*1000/1000), 'kg NO2/year')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This does not give us a very useful value considering the baseline. If we would want to bring this value in line with the baseline value by, for instance, reducing the lifetime of NO2 emissions at the plant location (1 sq km) this would lead to a lifetime of NO2 emissions at the plant location of a couple of minutes.    \n\nReading about conversion from NO2 density measurements to absolute values of NO2 does not give straightforward conversion methods. Indications are that density measurements are useful for comparing emissions as it is used in estimating the distribution of electricity generation between plants.\nAlso it may be used to confirm/question the absolute emission factors between plants. For this purpose we calculate the emission factors per MWH for coal/gas and oil as seen in the measurements and compare these with the reported values (see table below).\n\nThere is no differentation in the measured values for emission per MWh between coal, gas and oil plants."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"kengetal=gray_plant.loc[:,['primary_fuel','EF_hist','EF_abs','hist_prod_day']].groupby(by='primary_fuel').agg('sum')\nkengetal['EF_measured_MWh']=kengetal['EF_hist']/kengetal['hist_prod_day']\nkengetal['EF_abs_MWh']=kengetal['EF_abs']/kengetal['hist_prod_day']\nkengetal","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"raw","source":"My interpretation from these results is that the measurements do not have sufficient resolution to show such a distinction possibly due to the fact that emissions spread quite fast and/or dissappear in the conversion to other substances or due to weather conditions. A similar effect is observed at small plants such as the Vieques power plant. The emissions of this small plant seem to 'drown' in the emissions from the surroundings leading to a prediction of power generation which is much higher then the available capacity."},{"metadata":{},"cell_type":"markdown","source":"With the knowledge from the reported distribution of oil/gas/coal/renewables of 40%/40%/18%/2%, a better calculated distribution may be obtained by a different choice for calculating individual plant contributions within a cluster.  \nInstead of using a attribution to plants pro rato the capacity within a cluster a preference for gas is used. That is: all gas plants in a cluster are set at maximum production (only if this does not exceed calculated cluster production). The remaining production in the cluster is attributed to oil and coal.  \n\nThis gives the following table and a fuel usage distribution coal/gas/oil of 10%/39%/51% : which is closer to the reported distribution for gas but still needs tweaking to get the distribution for coal and oil in line with reported values. It shows that the model has the capability of bringing calculated historical production and emission figures more in line with reported figures if additional and relevant information is added. \n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# post-processing - for each plant\ngray_plant=gray_plant.iloc[:,:8]\n\ngray_plant=pd.merge(gray_plant, gray[['hist2_prod_day','Max_Prod_clst','cluster','EF_hist2']], how='left', on='cluster')\n\nfor x in gray_plant.cluster :\n    gray_plant.loc[((gray_plant.primary_fuel=='Gas')&(gray_plant.cluster==x)),'scen_gas_max']=gray_plant['Max_Prod']\n    if len(gray_plant.loc[(gray_plant.primary_fuel!='Gas')&(gray_plant.cluster==x),:])>0 :\n        gray_plant.loc[((gray_plant.primary_fuel!='Gas')&(gray_plant.cluster==x)),'scen_gas_max']=((gray.loc[gray.cluster==x,'hist2_prod_day'].values-gray_plant.loc[((gray_plant.primary_fuel=='Gas')&(gray_plant.cluster==x)),'scen_gas_max'].sum())*gray_plant['Max_Prod']/gray_plant.loc[((gray_plant.primary_fuel!='Gas')&(gray_plant.cluster==x)),'Max_Prod'].sum()).astype(int)\n    \n    if gray_plant.loc[((gray_plant.primary_fuel=='Gas')&(gray_plant.cluster==x)),'scen_gas_max'].sum() > gray.loc[gray.cluster==x,'hist2_prod_day'].values:\n        gray_plant.loc[((gray_plant.primary_fuel=='Gas')&(gray_plant.cluster==x)),'scen_gas_max']=(gray_plant['Max_Prod']*gray.loc[gray.cluster==x,'hist2_prod_day'].values/gray_plant.loc[((gray_plant.primary_fuel=='Gas')&(gray_plant.cluster==x)),'Max_Prod'].sum()).astype(int)     \n        gray_plant.loc[((gray_plant.primary_fuel!='Gas')&(gray_plant.cluster==x)),'scen_gas_max']=0\n        \n\ngray_plant['Activity_%']=(gray_plant['scen_gas_max']*100/gray_plant['Max_Prod']).astype(int)  \n\ngray_plant['EF_gas_max']=(gray_plant['Max_Prod']*gray_plant['EF_hist2']/gray_plant['Max_Prod_clst'])\n\n#gray_plant['EF_p_MWh']=(gray_plant['EF_hist']/gray_plant['hist_prod_day'])\n \n#gray_plant","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print('fuel distribution based on historical emissions')\nprint(gray_plant.groupby(by='primary_fuel').scen_gas_max.sum()/gray_plant.scen_gas_max.sum())\nprint('  ')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 5. Emission data for locations with and without a fossil fuel power plant #"},{"metadata":{},"cell_type":"markdown","source":"Further exploration into emission factor in other areas then fossil fuel power plants is done to find out if emissions at fossil fuel power plant locations stand out from other areas. To this end the full table of power plants is analyzed including hydro, wind and solar locations from which no emissions are expected. "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"gray_plant=power_plants.loc[:,['name','primary_fuel','capacity_mw','img_idx_lt','img_idx_lg','plant']]\n\n\ngray_plant=gray_plant.sort_values('primary_fuel')\n\n# calculation of maximum daily production (MWh) based on max. capacity of power plants\ngray_plant['Max_Prod']=(gray_plant['capacity_mw']*24).astype(int)\n\nprint('Dataframe for local calculation of emission factors : ')\n#print(gray_plant.head())\n\nn=11 # size of area around a plant where emission is averaged for further calculation\n#n=4\n#n=3\n\ngray=pd.DataFrame({})\nfor j in range(0,len(gray_plant)):\n    idx_lt=gray_plant.iloc[j,3]\n    idx_lg=gray_plant.iloc[j,4]\n    temp=gray_plant.loc[(((gray_plant.img_idx_lt>idx_lt-n) & (gray_plant.img_idx_lt<idx_lt+n)) & \n                         ((gray_plant.img_idx_lg>idx_lg-n) & (gray_plant.img_idx_lg<idx_lg+n))),gray_plant.columns[:7]]\n    if len(temp)>1 : # if there are more plants in 1 location, group them together into 1 source\n        temp1=pd.DataFrame([[str(len(temp))+' locations','comb',np.int(temp.capacity_mw.sum()),np.int(temp.img_idx_lt.mean()),np.int(temp.img_idx_lg.mean()),\n                           temp.plant.sum(),temp.Max_Prod.sum()]], columns=temp.columns)\n        gray_plant.loc[(((gray_plant.img_idx_lt>idx_lt-n) & (gray_plant.img_idx_lt<idx_lt+n)) & \n                         ((gray_plant.img_idx_lg>idx_lg-n) & (gray_plant.img_idx_lg<idx_lg+n))),'cluster']=temp.plant.sum()\n        gray=pd.concat([gray,temp1])\n    else :    # if there is 1 plant in a location, keep this one as the source\n        gray=pd.concat([gray,temp])\n        gray_plant.loc[(((gray_plant.img_idx_lt>idx_lt-n) & (gray_plant.img_idx_lt<idx_lt+n)) & \n                         ((gray_plant.img_idx_lg>idx_lg-n) & (gray_plant.img_idx_lg<idx_lg+n))),'cluster']=temp.plant.sum()\n        \ngray=gray.drop_duplicates()\n\ngray=gray.rename(columns = { 'plant': 'cluster', 'Max_Prod': 'Max_Prod_clst'})\n\nprint(gray)\n\nplot_locations_model(n)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# emission value in proximity of all plants with all locations in location mask - proximity is +/- n points from location of plant\nno2=[]\ntrop=[]\naai=[]\nclfr=[]\nfor j in range(0,len(gray)):\n    idx_lt=gray.iloc[j,3]\n    idx_lg=gray.iloc[j,4]\n        \n    no2_j=[]\n    trop_j=[]\n    aai_j=[]\n    clfr_j=[]\n    for i in range(0,len(trop_arr)):\n        no2_j.append(np.nanmean(no2_col_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\n        trop_j.append(np.nanmean(trop_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\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        clfr_j.append(np.nanmean(clfr_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\n        \n    no2.append(no2_j)\n    trop.append(trop_j)\n    aai.append(aai_j)\n    clfr.append(clfr_j)\n    \n    \n# list to DataFrame with dates   \nno=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow})\ntro=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow}) \naa=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow})\nclf=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow})\n\nfor j in range(0,len(gray)): #add average of emission for each location (plant) to dataframe with in the column name taken from df gray.cluster\n    no[gray.iloc[j,5]]=no2[j]\n    tro[gray.iloc[j,5]]=trop[j]\n    aa[gray.iloc[j,5]]=aai[j]  #add average of aai for each location (plant) to dataframe with in the column name taken from df gray.cluster\n    clf[gray.iloc[j,5]]=clfr[j]\n    \nprint('size of dataframe with emission data for gray-energy power-plant locations: ',no.shape, tro.shape, aa.shape, clf.shape)\n\n# sorting dataframe on date to produce ordered time series\nno=no.sort_values('key_date') ; tro=tro.sort_values('key_date') ; aa=aa.sort_values('key_date') ; clf=clf.sort_values('key_date')\n\nno=no.reset_index() ; tro=tro.reset_index() ; aa=aa.reset_index() ; clf=clf.reset_index()\nno=no.drop(columns=['index']) ; tro=tro.drop(columns=['index']) ; aa=aa.drop(columns=['index']) ; clf=clf.drop(columns=['index'])\nno=no.fillna(0) ; tro=tro.fillna(0) ; aa=aa.fillna(0) ; clf=clf.fillna(0)\n    \n#emission_correlation() # show locations with high correlation\n\n#plot_timeview_local_emission() # show timeview of emission per location\n\nplot_locationsview_local_data(no,no2_rgn,'NO2 column',tro,trop_rgn, 'NO2 in troposphere') # show emissions per location vs average for the region\nplot_locationsview_local_data(aa,aai_rgn, 'Aerosol absorption index (aai)',clf,clfr_rgn, 'Cloud fraction') # show emissions per location vs average for the region\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Measured values of emission on fossil fuel power plant locations do not stand out from locations without fossil fuel power plants. The same comment holds for the box-plots below.   \n\nMain conclusion is that locations of fossil fuel power plants cannot be distinguished solely from emission measurements of NO2 in troposphere. Additional information on locations of fossil fuel power plants is necessary to construct a model."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\ndata_dow=ef_rgn.pivot(columns='dow',values='trop')\n\ndays=['mon','tue', 'wed','thu','fri','sat','sun']\n\nf11, ax11 = plt.subplots(1, 2, figsize=(20, 4))\nchart=sns.boxplot(data=data_dow, ax=ax11[0])\nchart.set_xticklabels(days, rotation=45, horizontalalignment='right')\nchart.set_title('Regional NO2 in troposphere')\nchart.set_ylabel('emission (NO2 troposphere)') \n\np=13 # number selects location to show\ndata_dow=tro.pivot(columns='dow',values=tro.columns[p])\n\nchart=sns.boxplot(data=data_dow, ax=ax11[1])\nchart.set_xticklabels(days, rotation=45, horizontalalignment='right')\nchart.set_title(tro.columns[p])\nchart.set_ylabel('emission (NO2 troposphere)') ;\n\nmonthly_data_location(no,tro,'monthly data NO2 column','monthly data troposphere' )\nmonthly_data_location(aa,clf,'monthly data aerosol absorption index','monthly data cloud fraction' )","execution_count":null,"outputs":[]},{"metadata":{"colab_type":"text","id":"1cWq-KgX608w"},"cell_type":"markdown","source":"# 6. Climate data on temperature, humidity and wind speed #"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gray_plant=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']]\n\n\ngray_plant=gray_plant.sort_values('primary_fuel')\n\n# calculation of maximum daily production (MWh) based on max. capacity of power plants\ngray_plant['Max_Prod']=(gray_plant['capacity_mw']*24).astype(int)\n\nprint('Dataframe for local calculation of emission factors : ')\n#print(gray_plant.head())\n\nn=11 # size of area around a plant where emission is averaged for further calculation\n#n=4\n#n=3\n\ngray=pd.DataFrame({})\nfor j in range(0,len(gray_plant)):\n    idx_lt=gray_plant.iloc[j,3]\n    idx_lg=gray_plant.iloc[j,4]\n    temp=gray_plant.loc[(((gray_plant.img_idx_lt>idx_lt-n) & (gray_plant.img_idx_lt<idx_lt+n)) & \n                         ((gray_plant.img_idx_lg>idx_lg-n) & (gray_plant.img_idx_lg<idx_lg+n))),gray_plant.columns[:7]]\n    if len(temp)>1 : # if there are more plants in 1 location, group them together into 1 source\n        temp1=pd.DataFrame([[str(len(temp))+' locations','comb',np.int(temp.capacity_mw.sum()),np.int(temp.img_idx_lt.mean()),np.int(temp.img_idx_lg.mean()),\n                           temp.plant.sum(),temp.Max_Prod.sum()]], columns=temp.columns)\n        gray_plant.loc[(((gray_plant.img_idx_lt>idx_lt-n) & (gray_plant.img_idx_lt<idx_lt+n)) & \n                         ((gray_plant.img_idx_lg>idx_lg-n) & (gray_plant.img_idx_lg<idx_lg+n))),'cluster']=temp.plant.sum()\n        gray=pd.concat([gray,temp1])\n    else :    # if there is 1 plant in a location, keep this one as the source\n        gray=pd.concat([gray,temp])\n        gray_plant.loc[(((gray_plant.img_idx_lt>idx_lt-n) & (gray_plant.img_idx_lt<idx_lt+n)) & \n                         ((gray_plant.img_idx_lg>idx_lg-n) & (gray_plant.img_idx_lg<idx_lg+n))),'cluster']=temp.plant.sum()\n        \ngray=gray.drop_duplicates()\n\ngray=gray.rename(columns = { 'plant': 'cluster', 'Max_Prod': 'Max_Prod_clst'})\n\nprint(gray)\n\nplot_locations_model(n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# emission value in proximity of all plants with all locations in location mask - proximity is +/- n points from location of plant\nno2=[]\ntrop=[]\naai=[]\nclfr=[]\nfor j in range(0,len(gray)):\n    idx_lt=gray.iloc[j,3]\n    idx_lg=gray.iloc[j,4]\n        \n    no2_j=[]\n    trop_j=[]\n    aai_j=[]\n    clfr_j=[]\n    for i in range(0,len(trop_arr)):\n        no2_j.append(np.nanmean(no2_col_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\n        trop_j.append(np.nanmean(trop_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\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        clfr_j.append(np.nanmean(clfr_arr[i][idx_lt-n:idx_lt+n,idx_lg-n:idx_lg+n])) # calculate average of aai for location of plant\n        \n    no2.append(no2_j)\n    trop.append(trop_j)\n    aai.append(aai_j)\n    clfr.append(clfr_j)\n    \n    \n# list to DataFrame with dates   \nno=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow})\ntro=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow}) \naa=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow})\nclf=pd.DataFrame({'key_date':np.array(ef_first_key), 'first': ef_day,'dow':ef_dow})\n\nfor j in range(0,len(gray)): #add average of emission for each location (plant) to dataframe with in the column name taken from df gray.cluster\n    no[gray.iloc[j,5]]=no2[j]\n    tro[gray.iloc[j,5]]=trop[j]\n    aa[gray.iloc[j,5]]=aai[j]  #add average of aai for each location (plant) to dataframe with in the column name taken from df gray.cluster\n    clf[gray.iloc[j,5]]=clfr[j]\n    \nprint('size of dataframe with emission data for gray-energy power-plant locations: ',no.shape, tro.shape, aa.shape, clf.shape)\n\n# sorting dataframe on date to produce ordered time series\nno=no.sort_values('key_date') ; tro=tro.sort_values('key_date') ; aa=aa.sort_values('key_date') ; clf=clf.sort_values('key_date')\n\n#aa=pd.merge(aa,aai_rgn.loc[:,['key_date','aai_rgn']], how='inner',on='key_date')\n\nno=no.reset_index() ; tro=tro.reset_index() ; aa=aa.reset_index() ; clf=clf.reset_index()\nno=no.drop(columns=['index']) ; tro=tro.drop(columns=['index']) ; aa=aa.drop(columns=['index']) ; clf=clf.drop(columns=['index'])\nno=no.fillna(0) ; tro=tro.fillna(0) ; aa=aa.fillna(0) ; clf=clf.fillna(0)\n#aa.head()\n    ","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['wind']=np.sqrt(np.multiply(gfs_rgn['u_wind'],gfs_rgn['u_wind'])+np.multiply(gfs_rgn['v_wind'],gfs_rgn['v_wind']))\ngfs_rgn['log_wind']=np.log(gfs_rgn['wind'])\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, wind_speed, precipitable_water (all left axis), NO2 in troposhere 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(ef_rgn.iloc[:,6], ef_rgn.iloc[:,3], label='average '+ef_rgn.columns[3] , 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(ef_rgn.iloc[:,6], ef_rgn.iloc[:,3], label='average '+ef_rgn.columns[3] , 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(ef_rgn.iloc[:,6], ef_rgn.iloc[:,3], label='average '+ef_rgn.columns[3] , color='r')\nax3.legend() ; ax32.legend()\nax4 = fig10.add_subplot(324)\nax4.plot(gfs_rgn.iloc[:,1], gfs_rgn.iloc[:,9], label='average wind', color='b')\nax42 = ax4.twinx()\nax42.plot(ef_rgn.iloc[:,6], ef_rgn.iloc[:,3], label='average '+ef_rgn.columns[3] , color='r')\nax4.legend() ; ax42.legend()\nax5 = fig10.add_subplot(325)\nax5.plot(gfs_rgn.iloc[:,1], gfs_rgn.iloc[:,7], label='average precipitable_water', color='b')\nax52 = ax5.twinx()\nax52.plot(ef_rgn.iloc[:,6], ef_rgn.iloc[:,3], label='average '+ef_rgn.columns[3] , color='r')\nax5.legend() ; ax52.legend();\n# ax6 = fig10.add_subplot(326)\n# ax6.plot(gfs_rgn.iloc[:,1], gfs_rgn.iloc[:,7], label='average precipitable_water', color='b')\n# ax62 = ax6.twinx()\n# ax62.plot(ef_rgn.iloc[:,6], ef_rgn.iloc[:,3], label='average '+ef_rgn.columns[3] , color='r')\n# ax6.legend() ; ax62.legend()","execution_count":null,"outputs":[]},{"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[:,:7]\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#gray.loc[:,'EF_wght']=gray_location.hist2_EF_MWh_day/gray_location.hist2_EF_MWh_day.sum()\n\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['temp_2d']=ww['temp']+ww['temp'].shift(1)\n#    ww['temp_3d']=ww['temp']+ww['temp'].shift(1)+ww['temp'].shift(2)\n#    ww['temp_5d']=ww['temp']+ww['temp'].shift(1)+ww['temp'].shift(2)+ww['temp'].shift(3)+ww['temp'].shift(4)   \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    ww['wind']=np.sqrt(np.multiply(ww['u_wind'],ww['u_wind'])+np.multiply(ww['v_wind'],ww['v_wind']))\n#    ww['wind_2d']=ww['wind']+ww['wind'].shift(1)\n#    ww=ww.drop(columns=['u_wind','v_wind'])\n#    ww['wind_^3']=np.multiply(np.multiply(ww['wind'],ww['wind']),ww['wind'])\n#    ww['wind_3d']=ww['wind']+ww['wind'].shift(-1)+ww['wind'].shift(-2)\n#    ww['pr_water_2d']=ww['pr_water']+ww['pr_water'].shift(1)\n#    ww['pr_water^3']=np.multiply(np.multiply(ww['pr_water'],ww['pr_water']),ww['pr_water'])\n#    ww['pr_water_3d']=ww['pr_water']+ww['pr_water'].shift(1)+ww['pr_water'].shift(2)\n    \n    x=ww.groupby(by='key_date').agg(['max','min','mean','std'])\n    \n    X=pd.merge(clf.loc[:,['key_date',gray.iloc[j,5]]], x, how='inner', on='key_date')  #addition of cloud fraction\n    X=X.rename(columns = {gray.iloc[j,5]:'clf'})\n    \n#     c=gray.iloc[j,5]   #'EF_'+gray.iloc[j,5]\n#     X[c]=np.ones((len(X)))*gray.iloc[j,7] # addition of EF_wght for each plant to the dataframe\n# #    X['EF_distribution']=np.ones((len(X)))*gray.iloc[j,7] # addition of EF_wght for each plant to the dataframe\n\n\n#   XX=pd.concat([XX,X], axis=0, sort=False) # aggregation of dataframe per plantlocation into one dataframe\n\n    X=pd.merge(tro.loc[:,['key_date',gray.iloc[j,5]]], X, how='inner', on='key_date')  #addition of no2 troposphere\n    X=X.rename(columns = {gray.iloc[j,5]:'trop'})\n    \n    c=gray.iloc[j,5]   #'EF_'+gray.iloc[j,5]\n    X[c]=np.ones((len(X)))*gray.iloc[j,7] # addition of EF_wght for each plant to the dataframe\n#    X['EF_distribution']=np.ones((len(X)))*gray.iloc[j,7] # addition of EF_wght for each plant to the dataframe\n    \n  \n    XX=pd.concat([XX,X], axis=0, sort=False) # aggregation of dataframe per plantlocation into one dataframe\n    \n        \nXX=XX.fillna(0) \nXX=XX.reset_index()\n\nprint('DataFrame with local features from climate factors and presence of specific power plants one-hot encoded' )\nXX    \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 7. A prediction model for emission with local emission and weather data  "},{"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 (for this calculation NO2 in troposphere is taken as input)  \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 the factor R is expanded to 5 climate variabeles (temperature, relative and specific humidity, 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 measured emission (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\n\nFor features with a high correlation one is kept and the other removed to have a set of independent features.\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"y=XX['trop']\n\nX=XX.drop(columns=['index','key_date','trop'])\n\n# removing features with high correlation\nfeatures=X.columns\n\ncounter = 0\nto_remove = []\nfor feat_a in features:\n    for feat_b in features:\n        if feat_a != feat_b and feat_a not in to_remove and feat_b not in to_remove:\n            c = np.corrcoef(X[feat_a], X[feat_b])[0][1]\n            if c > 0.90:  #c > 0.995\n                counter += 1\n                to_remove.append(feat_b)\n                print('{}: FEAT_A: {} FEAT_B: {} - Correlation: {}'.format(counter, feat_a, feat_b, c))\n                \nprint('Features that are removed : ', to_remove)\n\nXT=X.drop(columns=to_remove)\nfeature_imp=pd.DataFrame({'feature': XT.columns})\n\nXT.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#effect of scaling: no visible difference in y_test vs y_pred, MSE is approx 0.8, not much difference in features\n\nfrom sklearn.preprocessing import StandardScaler\n\nscaler = StandardScaler()\n\nX=scaler.fit_transform(XT)\n\ny=scaler.fit_transform(np.ravel(y).reshape(-1, 1))\n\nfeature_imp=pd.DataFrame({'feature': XT.columns})\n\n#print(XT.shape,y.shape,X.shape)\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 = 4\nmin_child_weight = 5\nsubsample = 0.7\ncolsample_bytree = 0.7\nobjective = \"reg:squarederror\" #'reg:linear',#\"reg:squarederror\"\nnum_estimators = 500 #2000 #1000 #3000  #200\nlearning_rate =  0.1 #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=4, random_state=42, shuffle=True) # n_splits was 5\n\n\ni=0\ntestscore=[]\n#feature_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,:] # for unscaled dataframe\n    X_train, X_test = X[train_index,:], X[test_index,:] # for scaled version\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 emission (y_pred) vs measured value of emission (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":{"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 measured emission.\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. \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 of the variability in the model for predicting emission is explained by emission factors. Climate factors are not dominant in this model.\n\nWe use the feature importances of the model to calculate emission factors for the locations: \n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gray=gray.iloc[:,:7]\n\nprod_features=feature_imp.rename(columns= {'feature':'cluster'})\n\n#uit kfold\ngray=pd.merge(gray, prod_features.loc[:,['cluster','mean']], how='left', on='cluster')\ngray=gray.rename(columns= {'mean':'EF_contrib'}) # emission contribution as calculated from the feature importances of the model\ngray=gray.fillna(0)\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.EF_contrib.sum()*100)\n\ngray=gray.sort_values('EF_contrib', ascending=False)\n\ngray=gray.drop(columns=['img_idx_lt','img_idx_lg'])\n\n# # calculation of maximum daily production (MWh) based on max. capacity of power plants\n# gray['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_prod_day']=(gray['EF_contrib']*Prod_day/gray.EF_contrib.sum()).astype(int)\n\n# activity factor calculated from model emission distribution and maximum daily production\ngray['Activity_%']=(gray['hist_prod_day']*100/gray['Max_Prod_clst']).clip(upper=100).astype(int)\n\n# multiplier: if activity > 100 (%) then emissions cannot be explained by power production because the plant is running beyond maximum capacity.\n# the multiplier can explain these additional emissions attributable to choice of primary_fuel and generation of technology (age of plant).\ngray['multiplier']=(gray['hist_prod_day']/gray['Max_Prod_clst']).astype(int)\n\ngray['hist_prod_mdl1']=gray_cluster['hist_prod_day']\ngray['activity_mdl1']=gray_cluster['Activity_%']\n#gray_plant['multiplier_mdl1']=gray_cluster['multiplier']\ngray['multiplier_mdl1']=gray_cluster['multiplier']\n\ngray","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The table above is based on the weight of the feature importances. \n\nTo the right of the table three columns are added with the values from the earlier model based only on local emissions. \n\nIf the columns are compared one notices that the values are different but that a multiplier is calculated for the same locations as in the first model. A second step is still necessary to take into account capacity constraints of power plants. The model with climate features does not lead to a better distribution of emissions to the power plants. The observation is therefor that the model with the climate factors does not add much to the simpler model only using location. Because of this the advice is to use the more simple model.\n"}],"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":4}