{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook is a part of the solution for DSG: EIE competition. The solution splited into 4 parts. Here is the list of notebook in correct order. The part of solution you are currently reading is highlighted in bold.\n\n[1. Introduction to the solution of DSG: EIE](https://www.kaggle.com/niyamatalmass/1-introduction-to-the-solution-of-dsg-eie)\n\n[2. Sub-region and State wise E.F and Evaluation](https://www.kaggle.com/niyamatalmass/2-sub-region-and-state-wise-e-f-and-evaluation)\n\n[**3. Individual Power plant E.F and Evaluation**](https://www.kaggle.com/niyamatalmass/3-individual-power-plant-e-f-and-evaluation)\n\n[4. Final thoughts, recommendation](https://www.kaggle.com/niyamatalmass/4-final-thoughts-recommendation)\n***\n<br/>"},{"metadata":{},"cell_type":"markdown","source":"<h1 align=\"center\"><font color=\"#5831bc\" face=\"Comic Sans MS\">Individual Power plant E.F and Evaluation</font></h1> "},{"metadata":{},"cell_type":"markdown","source":"# Notebook Overview\nWe previously calculate emission and emission factor of Peurto Rico and other states, evaluate with 2018 EPA emissions dataset. We did that excluding pixel that doesn't fall on the power plant area. But we also discuss the trade-off we have to face for evaluation. In this notebook, we will calculate the emission and emission factor for an individual power plant for 2019. After that, we will evaluate with 2019 EPA power plantwise emission data. \n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np\nimport math \nimport pandas as pd\nfrom scipy.ndimage import gaussian_filter\nimport glob \nimport os\nimport time\nfrom tqdm import tqdm_notebook as tqdm\n\nimport geopandas \nimport rasterio as rio\nimport folium\nimport tifffile as tiff\n\nimport ee\nfrom kaggle_secrets import UserSecretsClient\nfrom google.oauth2.credentials import Credentials\n\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\n\nfrom plotly.offline import init_notebook_mode\ninit_notebook_mode(connected=True)\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\n\n################\n# below code is earth enginie initialization \n# paste it in different cell for initializaiton for yours personal use\n################\n\n# Trigger the authentication flow.\n# ee.Authenticate()\n\n\n# !cat ~/.config/earthengine/credentials\n\n\nuser_secret = \"earth_engine_2\" # Your user secret, defined in the add-on menu of the notebook editor\nrefresh_token = UserSecretsClient().get_secret(user_secret)\ncredentials = 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)\nee.Initialize(credentials=credentials)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Calculate emission for all coal power plant\nWe previously saw the methodology for calculating emission and emission factor for sub-region and individual states using satellite data. Also, we have tested our methodology against bottom-up emission and found a very good correlation. In this notebook, we are going to calculate the emission and emission factor for an individual power plant. In the above, we described why we should calculate emission and E.F for individual power plant. \n\nMethodology: Basically, all the method for calculating emission for an individual power plant is almost identical to our previous method for calculating emission factor for each state. Here we just select each power plant and draw a circle with a relative radius. The relative radius is calculated using the same method. After that, we excluded all the pixels that don't fall into our circle. And finally, reduce that AOI(circle) into total no2 using distance weighted approach. \n\n1. [Importing power plant database, filter them](#power_plant)\n2. [Create a function for calculating emission for all power plants](#function)\n3. [Create a batch of running for computation limit and calculate](#batch)\n4. [Calculate total no2 emission for each power plant](#calculate)\n5. [Process for evaluation](#process)"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"band_viz = {\n  'min': 0,\n  'max': 0.0002,\n  'palette': ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']\n}\n\nimport 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    \n# Define a method for displaying Earth Engine image tiles on a folium map.\ndef add_ee_layer(self, ee_object, vis_params, name):\n    \n    try:    \n        # display ee.Image()\n        if isinstance(ee_object, ee.image.Image):    \n            map_id_dict = ee.Image(ee_object).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n            ).add_to(self)\n        # display ee.ImageCollection()\n        elif isinstance(ee_object, ee.imagecollection.ImageCollection):    \n            ee_object_new = ee_object.mosaic()\n            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n            ).add_to(self)\n        # display ee.Geometry()\n        elif isinstance(ee_object, ee.geometry.Geometry):    \n            folium.GeoJson(\n            data = ee_object.getInfo(),\n            name = name,\n            overlay = True,\n            control = True\n        ).add_to(self)\n        # display ee.FeatureCollection()\n        elif isinstance(ee_object, ee.featurecollection.FeatureCollection):  \n            ee_object_new = ee.Image().paint(ee_object, 0, 2)\n            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n        )\n        # display ee.FeatureCollection()\n        elif isinstance(ee_object, ee.feature.Feature):  \n            ee_object_new = ee.Image().paint(ee_object, 0, 2)\n            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n        ).add_to(self)\n    \n    except:\n        print(\"Could not display {}\".format(name))\n    \ndef plot_ee_data_on_map(dataset,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': ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']}\n    my_map = folium.Map(location=[latitude,longitude], zoom_start=zoom, height=500)\n    my_map.add_ee_layer(dataset, vis_params, 'Color')\n    my_map.add_child(folium.LayerControl())\n    display(my_map)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='power_plant'></a>\n### 1. Importing power plant database, filter them\nAs like the previous notebook, we first start by importing our power plant database. Also, in here we imported 2019 EPA emission data. This dataset contains about 1400 power plants emission data for 2019. This makes this dataset perfect for us. Because we have full coverage Sentinel-5p satellite image for whole 2019. So if we calculate individual power plant emission from satellite images we can easily evaluate against this dataset, unlike the previous notebook where we evaluate our model data with slightly different dates. \n\nAfter importing power plant data, we filter them for containing power plants that also available in the EPA 2019 data and for computation easy we only select coal and hydropower plant in this notebook but this can be changed with simple code. "},{"metadata":{"trusted":true},"cell_type":"code","source":"df_2019_epa_emission_data = pd.read_csv('../input/2019-epa-nox-emission/epa_no2_emission_2019.csv')\npower_plant_feature_collection = ee.FeatureCollection(\"WRI/GPPD/power_plants\")\n\ncountry_filter = ee.Filter.eq('country', 'USA')\nname_filter = ee.Filter.inList('name', df_2019_epa_emission_data['Facility Name'].tolist())\nlatitude_filter = ee.Filter.inList('latitude', df_2019_epa_emission_data['Latitude (degrees)'].tolist())\nlongitude_filter = ee.Filter.inList('longitude', df_2019_epa_emission_data['Longitude (degrees)'].tolist())\ntype_filter = ee.Filter.inList('fuel1', ['Coal', 'Hydro'])\n\npower_plant_feature_collection = power_plant_feature_collection.filter(\n    ee.Filter.And(country_filter,type_filter,\n                  ee.Filter.Or([name_filter,\n                                latitude_filter,\n                                longitude_filter])))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='function'></a>\n### 2. Create a function for calculating emission for all power plants\nAs we are already familiar with the methodology we don't need much explanation. In the previous notebook, we explained the methodology by giving an individual example but in here we just calculate all power plant emission by one chance. For that, we just declare our function. The function is documented very well and explained every single part of the process. "},{"metadata":{"trusted":true},"cell_type":"code","source":"us_ton_conve_mult_factor = (1.10231e-6)\n\ndef calc_each_pp_total_no2(feature):\n    buffer_amount = ee.Number(feature.get('gwh_2016')).multiply(1.1)\n    buffer_point = feature.geometry().buffer(distance= buffer_amount, maxError=200)\n    collection_for_pp = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2')\\\n    .select('tropospheric_NO2_column_number_density')\\\n    .filterBounds(buffer_point)\n    \n    start = ee.Date('2019-01-01')\n    finish = ee.Date('2020-01-01')\n\n    diff = finish.difference(start, 'day')\n    range_date = ee.List.sequence(0, diff.subtract(1)).map(lambda day: start.advance(day,'day'))\n    \n    def day_mosaics(date, newlist):\n        date = ee.Date(date)\n        newlist = ee.List(newlist)\n        filtered = collection_for_pp.filterDate(date, date.advance(1,'day'))\n        image = ee.Image(filtered.mosaic())\n        return ee.List(ee.Algorithms.If(filtered.size(), newlist.add(image), newlist))\n    \n    collection_for_specific_pp = ee.ImageCollection(ee.List(range_date.iterate(day_mosaics, ee.List([]))))\n    \n    \n#     s5p_mosaic_for_each_states = s5p_mosaic_for_each_states.map(lambda image: image.reproject(\n#         crs=collection.first().projection().crs().getInfo(),\n#         scale = collection.first().projection().nominalScale().getInfo()))\n    \n    #######\n    maxDist = 20000\n    distance = power_plant_feature_collection.distance(searchRadius=maxDist, maxError=1000)\n    weight = distance.subtract(maxDist).abs().divide(maxDist)\n    \n    \n    #######\n    \n    \n    collection_for_specific_pp = collection_for_specific_pp.map(lambda image: image.multiply(24500000))\n    collection_for_specific_pp = collection_for_specific_pp.map(lambda image: image.multiply(46))\n    collection_for_specific_pp = collection_for_specific_pp.map(lambda image: image.multiply(us_ton_conve_mult_factor))\n    collection_for_specific_pp = collection_for_specific_pp.reduce(ee.Reducer.sum())\n    \n    collection_for_specific_pp = collection_for_specific_pp.addBands(weight)\n    \n    \n    ### applying gaussian\n    gauss_kernel = ee.Kernel.gaussian(\n        radius= 3, sigma= 1,\n        units= 'pixels')\n    collection_for_specific_pp = collection_for_specific_pp.convolve(gauss_kernel)\n    \n    \n    total_no2_yearly_specific_pp = collection_for_specific_pp.reduceRegion(reducer = ee.Reducer.sum().splitWeights(),\n                                                  geometry= buffer_point,\n                                                  scale= 1113)\n    \n    return feature.set({'total_no2' : total_no2_yearly_specific_pp.get('sum')})\n#     return feature.set({'total_no2' : total_no2_yearly_specific_pp.get('list').get('tropospheric_NO2_column_number_density_sum')})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='batch'></a>\n### 3. Create a batch of running for computation limit and calculate\nFor avoiding computation limit, we batched our computation in different parts. Though in this notebook, we don't need that but still we keep the functionality in the case in future we need it. "},{"metadata":{"trusted":true},"cell_type":"code","source":"all_pp_code = power_plant_feature_collection.reduceColumns(\n    reducer=ee.Reducer.toList(),\n    selectors=['gppd_idnr']).getInfo().get('list')\n\nfirst_batch_pp_code = all_pp_code[:len(all_pp_code)//1]\nsecond_batch_pp_code = all_pp_code[len(all_pp_code)//2:]\n\n\nfirst_batch_pp_collection = power_plant_feature_collection\\\n.filter(ee.Filter.inList('gppd_idnr', first_batch_pp_code))\n\nsecond_batch_pp_collection = power_plant_feature_collection\\\n.filter(ee.Filter.inList('gppd_idnr', second_batch_pp_code))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"first_batch_total_no2 = first_batch_pp_collection.map(lambda feature: calc_each_pp_total_no2(feature))\n# second_batch_total_no2 = second_batch_pp_collection.map(lambda feature: calc_each_pp_total_no2(feature))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='calculate'></a>\n### 4. Calculate total no2 emission for each power plant\nFinally, using the function we created earlier we calculate all power plant emission using a distance weighted sum approach. How the function works is describes thoroughly in function with comment. "},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nfirst_batch_results = first_batch_total_no2.reduceColumns(\n    reducer=ee.Reducer.toList().repeat(5),\n    selectors=['gppd_idnr', 'name', 'latitude', 'longitude', 'total_no2']).getInfo()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# second_batch_results = second_batch_total_no2.reduceColumns(\n#     reducer=ee.Reducer.toList().repeat(5),\n#     selectors=['gppd_idnr', 'name', 'latitude', 'longitude', 'total_no2']).getInfo()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='process'></a>\n### 5. Process for evaluation\nFinally, we have calculated emission for each power plant. But to evaluate, we have to process the output a little bit. First, we convert our function returns to dataframe and then merge with 2019 power plant wise data. "},{"metadata":{"trusted":true},"cell_type":"code","source":"df_first_batch_total_no2 = pd.DataFrame(columns=['gppd_idnr', 'name', 'latitude', 'longitude', 'total_no2'])\n# df_second_batch_total_no2 = pd.DataFrame(columns=['gppd_idnr', 'name', 'latitude', 'longitude', 'total_no2'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def convert_feature_collection_results_to_df(results, dataframe):\n    dataframe['gppd_idnr'] = pd.Series(results.get('list')[0])\n    dataframe['name'] = pd.Series(results.get('list')[1])\n    dataframe['latitude'] = pd.Series(results.get('list')[2])\n    dataframe['longitude'] = pd.Series(results.get('list')[3])\n    dataframe['total_no2'] = pd.Series(results.get('list')[4])\n    return dataframe","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_first_batch_total_no2 = convert_feature_collection_results_to_df(first_batch_results, df_first_batch_total_no2)\n# df_second_batch_total_no2 = convert_feature_collection_results_to_df(second_batch_results, df_second_batch_total_no2)\n\n# df_pp_total_no2 = pd.concat([df_first_batch_total_no2, df_second_batch_total_no2])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = df_first_batch_total_no2.merge(df_2019_epa_emission_data, left_on=['latitude', 'longitude'],\n                               right_on=['Latitude (degrees)', 'Longitude (degrees)'], how='inner')\n\n\n\n\ntemp.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Evaluation\nIn this step, we are going to evaluate our model output for individual power plant. In the previous notebook, we evaluated all states emission factor with bottom-up emission using 2018. But for evaluating power plant emission factor we have full power plant wise data available from EPA for 2019. \n\nWe first see how far our satellite calculated emission from local emission and then we will see some compare with the distribution. "},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = go.Figure()\nfig.add_trace(go.Scatter(\n    x=temp['total_no2'].values.tolist(),\n    y=temp['name'].values.tolist(),\n    marker=dict(color=\"crimson\", size=12),\n    mode=\"markers\",\n    name=\"Satellite Image(2019)\",\n))\n\nfig.add_trace(go.Scatter(\n    x=temp['2019 NOX Emissions (tons)'].values.tolist(),\n    y=temp['name'].values.tolist(),\n    marker=dict(color=\"gold\", size=12),\n    mode=\"markers\",\n    name=\"Bottom up(2019)\",\n))\n\nfig.update_layout(title=\"Satellite emission data vs bottom up emission data for individual power plants\",\n                  xaxis_title=\"NO2 amount (tons)\",\n                  yaxis_title=\"Power plant\" ,autosize=False,\n                  width=700,height=1500,yaxis= {'nticks': 100})\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are seeing very promising results. We are seeing that our model output is very correlated with bottom-up emission which means that our model performs very well. Now we can easily get the emission factor by dividing the emission of each power plant with its generation. "}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}