{"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\">Sub-region and State wise E.F and Evaluation</font></h1> "},{"metadata":{},"cell_type":"markdown","source":"# Notebook overview\nThis notebook implements the methodology to calculate the emission factor from power generation in sub-region and states level. In the first notebook of my solution, I describe the basic theory part of the methodology. In this notebook, I will describe the theory in more detail and implement the theory. After that, I will evaluate the results and see how our model perform in a real-world example. "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"##### \n# importing necessary libraries\n####\n\nimport 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\n\nimport matplotlib.pyplot as plt\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Abbreviation\n* GEE - Google Earth Engine\n* TVCD - Troposporic vertical column density\n* E.F - Emission Factor \n* AOI - Area Of Interest"},{"metadata":{},"cell_type":"markdown","source":"# Earth engine overview and why?\nWe will extensively use Google Earth Engine(GEE) for our research and implementation. Because Google Earth Engine is the most advanced cloud-based geospatial processing platform in the world! The purpose of the Earth Engine is to:\n\n* Perform highly-interactive algorithm development at a global scale\n* Push the edge of the envelope for big data in remote sensing\n* Enable high-impact, data-driven science\n* Make substantive progress on global challenges that involve large geospatial datasets\n\nBecause in our kaggle kernel we don't have enough compute resources, we will rely on GEE backend for our analysis. To understand the technical details and implementation, having some background in GEE will help. But I will try my best to explain important parts of the implementation. First, we have to authenticate and initialize our GEE python API. "},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"# Trigger the authentication flow. \n# ee.Authenticate()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# !cat ~/.config/earthengine/credentials","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# initializing earth engine\n\nuser_secret = \"earth_engine_3\" # 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":{"_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\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":"# Methodology for calculating E.F for sub-region and states\nOverview: For every power plant we draw a circle centring individual power plant with a radius that relative to power plant yearly emission(GWH). After that exclude all other pixel value of satellite image, as a result, we are left with individual power plants area with the circle. Then we sum all NO2 using a weighted approach where weight is defined how far a pixel is located from the centre( power plant ). This is just an overview. Below we break down each step with a clear explanation of our method. We are going to calculate the emission and emission factor for Peurto Rico and from that, we will learn how this methodology work.\n\n1. [Importing and filtering satellite images](#importing_satellite_images)\n2. [Importing power plant database, filter them and calculate AOI](#calculate_aoi)\n3. [Convert daily scattered image tiles to daily mosaic](#daily_mosaic)\n4. [Unit conversion](#unit_convert)\n5. [Building a yearly composite of AOI](#composite_aoi)\n6. [Finally, calculate total NO2 for a region using weighted reductions](#finally_weighted_reduction)"},{"metadata":{},"cell_type":"markdown","source":"<a id='importing_satellite_images'></a>\n### 1. Importing and filtering satellite images\nOur first step is to import our Sentinel-5p dataset and select TVCD band and filter bounds to Peurto Rico. \n\nWhy TVCD(Trophosphoric Vertical Column Density)? Because our satellite covers daily global coverage, we will be using tropospheric vertical column, because it is the nearest stage of the atmosphere where all weather takes place. Because no2 has very shot lifespan it is the ideal column for our research."},{"metadata":{"trusted":true},"cell_type":"code","source":"peurto_rico_state_code = ['72'] # peurto rico code\n\n# import tiger states boudary data for getting states boundary\npeurto_rico_geometry = ee.FeatureCollection(\"TIGER/2018/States\")\\\n.filter(ee.Filter.inList('STATEFP', peurto_rico_state_code))\n\n# import sentinel-5p dataset, filter band to TVCDs and filter bounds (peurto rico)\nno2_satellite_data_for_peurto_rico = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2')\\\n    .select('tropospheric_NO2_column_number_density')\\\n    .filterBounds(peurto_rico_geometry.geometry())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='calculate_aoi'></a>\n### 2. Importing power plant database, filter them and calculate AOI\nWe want to exclude all the pixels that don't fall in our AOI. But how do we calculate our AOI? We will take each power plant and calculate a circle with a relative radius where the radius is a factor of annual total emissions(GWH) of that power plant. We will do that for all power plants in desired states and exclude pixels that don't fall into the AOI.\n\nFirst, we will import power plant database from GEE. It stores as FeatureCollection. After that, we filter the power plant situated in the USA. In this data, power plant data stored as a point. Each point contains latitude and longitude of the power plant. We use that point to draw a circle of relative radius around the power plant. We will do that for each plant in the USA. Then we exclude other pixels that don't fall into the AOI using masking and store the data as a feature collection for future when we mask this image with no2 image for excluding pixels out of our AOI. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# import power plant database from GEE\npower_plant_feature_collection = ee.FeatureCollection(\"WRI/GPPD/power_plants\")\n\n# filter USA power plant only\ncountry_filter = ee.Filter.eq('country', 'USA') \npower_plant_feature_collection = power_plant_feature_collection.filter(country_filter)\n\ndef calc_buffer(feature):\n    \"\"\"\n    this function calculates the circle around each power plant using relative \n    radius and update the geometry\n    \"\"\"\n    keepProperties = ['country','country_lg','name','gppd_idnr','capacitymw',\n                      'latitude','longitude','fuel1','fuel2','fuel3','fuel4','comm_year',\n                      'owner','source','url','src_latlon','cap_year','gwh_2013',\n                      'gwh_2014','gwh_2015','gwh_2016','gwh_estimt']\n    buffer_amount = ee.Number(feature.get('gwh_2016')).multiply(1.1) # circle radius\n     # create that circle using buffer\n    buffer = feature.geometry().buffer(buffer_amount, maxError=200)\n    return ee.Feature(buffer).copyProperties(feature, keepProperties)\n\n# apply the above function to create circle around each power plant of USA\npower_plant_feature_collection_buffer = power_plant_feature_collection.map(\n    lambda feature: calc_buffer(feature))\n\n# create the mask for exlcuding pixels that don't fall in the circle\npp_masks = ee.Image().toByte().paint(power_plant_feature_collection_buffer, 1)\npp_masks = pp_masks.updateMask(pp_masks)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# let's plot the mask\nminimum_value = 0\nmaximum_value = 0.0002\nlatitude = 37.5010\nlongitude = -122.1899\nzoom =4.3\nplot_ee_data_on_map(pp_masks, minimum_value, maximum_value, latitude, longitude, zoom)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are seeing on the map that our AOI are showing in red, in where we will calculate no2. "},{"metadata":{},"cell_type":"markdown","source":"<a id='daily_mosaic'></a>\n### 3. Convert daily scattered image tiles to daily mosaic\n\nThese Sentinel-5p images come with different tiles. Sometimes a single location has multiple tiles. We will use earth engine mosaic feature to turn all tiles into a single image for each day for Peurto Rico as an example. And we will do that for all days between 2018-06-28 to 2019-06-28, generating a total of 365 images. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# yearly first date and last date \nstart = ee.Date('2018-06-28')\nfinish = ee.Date('2019-06-28')\n\ndiff = finish.difference(start, 'day') # 365 days \nrange_date = ee.List.sequence(0, diff.subtract(1)).map(lambda day: start.advance(day,'day'))\n\ndef day_mosaics(date, newlist):\n    \"\"\"\n    this function convert daily tiles of image of a region to a full mosaic image\n    \"\"\"\n    date = ee.Date(date)\n    newlist = ee.List(newlist)\n    filtered = no2_satellite_data_for_peurto_rico.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# generating 365 images of daily in a year\nno2_satellite_data_for_peurto_rico = ee.ImageCollection(ee.List(range_date.iterate(day_mosaics,\n                                                                                   ee.List([]))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='unit_convert'></a>\n### 4. Unit conversion\nSentinel-5p data comes in mol/m^2 unit. We want to convert that unit to locally based unit like a ton or lb. For now, let's convert mol/m^2 to a ton. But how? Let's understand the unit first. \n\nLet's talk about a single pixel. Sentinel-5p single-pixel size (area) is 7x3.5 km^2(24,500,000 m^2). Sentinel-5p data gives value for each pixel (e.g a pixel value is 0.89 mol/m^2). If we multiply that with each pixel area then we negate m^2 and we get unit as mol. Next, we convert mol to gram by multiplying 46 and gram to tons by multiplying gram_to_ton factor. We will apply these process to each pixel in sentinel-5p images and in that way we will get no2 emission in tons.\n\nI inspired this method conversion from [a. this paper](http://www.scielo.org.mx/scielo.php?script=sci_arttext&pid=S0187-62362018000300189), [b. question 1](https://gis.stackexchange.com/questions/353368/correlating-mol-m2-with-ton-and-lb), [c. question 2](https://earthscience.stackexchange.com/questions/19444/how-to-convert-mol-m2-to-total-mass-e-g-gram-kg-etc)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# multipication factor for converting gram to ton\nus_ton_conve_mult_factor = (1.10231e-6) \nno2_satellite_data_for_peurto_rico = no2_satellite_data_for_peurto_rico.map(\n    lambda image: image.multiply(24500000)) # convert mol/m^2 to mol\nno2_satellite_data_for_peurto_rico = no2_satellite_data_for_peurto_rico.map(\n    lambda image: image.multiply(46)) # convert mol no2 to gram no2\nno2_satellite_data_for_peurto_rico = no2_satellite_data_for_peurto_rico.map(\n    lambda image: image.multiply(us_ton_conve_mult_factor)) # convert gram to ton","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='composite_aoi'></a>\n### 5. Building a yearly composite of AOI\nNow we want to all those 365 images into one image using composite features of earth engine. Where each pixel is the yearly sum of no2. Then from that image, we will exclude pixels that don't fall in our AOI. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# convert 365 images to one image composite using ee.Reducer.sum()\nno2_satellite_data_for_peurto_rico = \\\nno2_satellite_data_for_peurto_rico.reduce(ee.Reducer.sum())\n\n# exlclude pixels that don't fall in the AOIs\nno2_satellite_data_for_peurto_rico = no2_satellite_data_for_peurto_rico.updateMask(pp_masks)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='finally_weighted_reduction'></a>\n### 6. Finally, calculate total NO2 from power plant using weighted reduction\nFinally, we will now calculate the total no2 of Puerto Rico. But we are not doing a normal reduction(sum), we will be using a weighted reduction by distance from the power plant as a weight. We can easily calculate sum reduction over an AOI but how do we calculate the weights. After researching a lot I have found the solution. \n\nWe want the weight to be calculated in a way where each pixel get less weight as it's moving away from the source ( power plant ). For implementation this we first declare max distance and then use GEE FeatureCollection distance function calculate the distance for each pixel and then calculate weights from 0 to 1. After we use GEE reduceRegion with sum reducer to reduce our state no2 to a single sum value. \n\nThis method is very useful because it takes a lot of factors considered. For example, wind moves no2 from its source. Within our AOI, our satellite will take pictures and if no2 moving away from its source it will have a lot of noise around the circle. For this, weighted distance gives less weight as no2 move from its source. By this way, we can minimize a lot of noise. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculating weights \nmaxDist = 5000\ndistance = power_plant_feature_collection.distance(searchRadius=maxDist, maxError=100)\nweight = distance.subtract(maxDist).abs().divide(maxDist)\n# adding weights as a band to our no2 composite images \nno2_satellite_data_for_peurto_rico = no2_satellite_data_for_peurto_rico.addBands(weight)\n\n# we will calculate total no2 of a region, using reduceRegion function\n# but first we will calculate without weight\ntotal_no2_peurto_rico = no2_satellite_data_for_peurto_rico.reduceRegion(\n    reducer = ee.Reducer.sum(),\n    geometry= peurto_rico_geometry.geometry(),\n    scale= 1113).getInfo()\nprint('yearly total no2 emission from power generation (ton) without weight: ' \\\n      + str(total_no2_peurto_rico.get('tropospheric_NO2_column_number_density_sum'))) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"total_no2_peurto_rico = no2_satellite_data_for_peurto_rico.reduceRegion(reducer = ee.Reducer.sum().splitWeights(),\n                                                  geometry= peurto_rico_geometry.geometry(),\n                                                  scale= 1113).getInfo()\nprint('yearly total no2 emission from power generation (ton) with weight: ' \\\n      + str(total_no2_peurto_rico.get('sum')))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, we get annual total emission from power generation in Puerto Rico. As we already have data on yearly total energy generation we can easily calculate the emission factor. But how do we evaluate the result? For this reason, we are going to do another major step. "},{"metadata":{},"cell_type":"markdown","source":"# How do we evaluate our methodology?\n\nTo know how correct and efficient a model is it crucial to evaluate the model with a real-world scenario. For this reason, we have developed a very well defined evaluation method to test our model. As we know that Sentinel-5p has data from 2018-06-28, so we cant compare our calculated value against historical value, so how do we evaluate our model. For solving that, what we did is calculate emission for all states and compare them with 2018 and 2019 data. That way we can evaluate our results perfectly. \n\nOur searching and reading a lot, we have found only two datasets that useful for us. The first one is bottom-up emission data from EPA and the second individual power plant emission data from EPA. The problem is EPA up until right now doesn’t publish statewide emission data for 2019. the [latest data is for 2018](https://www.epa.gov/energy/emissions-generation-resource-integrated-database-egrid). And this dataset contains the yearly sum of no2 emission, emission factor for each state, plant on a total yearly average. But our sentinel-5p satellite data is available from July 2018. That has a six-month gap from the total yearly 2018 data. \n\nOn the other hand, we have [no2 data for 2019](https://www.epa.gov/airmarkets/power-plant-emission-trends) that contains emission from power plant wise and but it doesn't have all power plant emission. That makes it unsuitable for our evaluation. \n\nThen what to do? Don’t worry! After some researching, we have found a solution for that. \n    • 2018 EPA contains information in ozone season. That is good news for us. Because mostly ozone season happens after 5 -6 month. So we have calculated total emission from satellite imagery in that date bound and evaluate against ozone emission. This will give us an accurate understanding of our methodology. \nWe will use the 2019 power plant wise emission data to evaluate our power plant wise emission calculation ( marginal emission ) in the next notebook. "},{"metadata":{},"cell_type":"markdown","source":"# Calculate emission and E.F for all US states\nFor evaluation, we like to calculate emission and emission factor for us states. This will help us to more accurately measure our methodology performances.\n\n1. [Create a function for calculating emission for all US states](#function)\n2. [Create a batch of running for computation limit and calculate](#batch)\n3. [Process for evaluation](#process_evaluation)"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"<a id='function'></a>\n### 1. Create a function for calculating emission for all US states\nFor calculating emission for all us states we are going to use a function for easy to use. This function works the same way as we describe above. We just wrap it into a function and apply for all us states. "},{"metadata":{"trusted":true},"cell_type":"code","source":"us_ton_conve_mult_factor = (1.10231e-6)\n\ndef calc_each_states_total_no2(feature): \n    \n    \"\"\"\n    this function calculates emission for each states. \n    here feature represents a single states \n    all the functionality are same described above\n    \"\"\"\n    \n    collection_for_specific_states = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2')\\\n    .select('tropospheric_NO2_column_number_density')\\\n    .filterBounds(feature.geometry())\n    \n    start = ee.Date('2018-06-28')\n    finish = ee.Date('2018-11-28')\n\n    diff = finish.difference(start, 'day')\n    range_date = ee.List.sequence(0, diff.subtract(1)).map(\n        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_specific_states.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    s5p_mosaic_for_each_states = ee.ImageCollection(\n        ee.List(range_date.iterate(day_mosaics, ee.List([]))))\n    \n#     s5p_mosaic_for_each_states = s5p_mosaic_for_each_states.map(\n        #lambda image: image.convolve(gauss_kernel))\n    \n    s5p_mosaic_for_each_states = s5p_mosaic_for_each_states.map(\n        lambda image: image.multiply(24500000))\n    s5p_mosaic_for_each_states = s5p_mosaic_for_each_states.map(\n        lambda image: image.multiply(46))\n    s5p_mosaic_for_each_states = s5p_mosaic_for_each_states.map(\n        lambda image: image.multiply(us_ton_conve_mult_factor))\n    s5p_mosaic_for_each_states = s5p_mosaic_for_each_states.reduce(ee.Reducer.sum())\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    s5p_mosaic_for_each_states = s5p_mosaic_for_each_states.addBands(weight)\n    #########\n    \n    s5p_mosaic_for_each_states = s5p_mosaic_for_each_states.updateMask(pp_masks)\n        \n    total_no2_yearly_each_states = s5p_mosaic_for_each_states.reduceRegion(\n        reducer = ee.Reducer.sum().splitWeights(),\n        geometry= feature.geometry(),scale= 1113)\n    \n    return feature.set({'total_no2' : total_no2_yearly_each_states.get('sum')})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='batch'></a>\n### 2. Create a batch of running for computation limit and calculation\nBecause we have more than 50 states to calculate it take quite a few time to calculate. To avoid GEE time out error, we separate our calculation in two batches. The first batch will calculate emission for half states and the second batch will calculate for other half states. \n\nAfter that, we will calculate emission for all states. "},{"metadata":{"trusted":true},"cell_type":"code","source":"us_states_collection = ee.FeatureCollection(\"TIGER/2018/States\")\nall_us_state_code = us_states_collection.reduceColumns(\n    reducer=ee.Reducer.toList().repeat(2),\n    selectors=['STATEFP', 'NAME']).getInfo().get('list')[0]\n\nfirst_batch_us_state_code = all_us_state_code[:len(all_us_state_code)//2]\nsecond_batch_us_state_code = all_us_state_code[len(all_us_state_code)//2:]\n\n# create feature collection for half of states\nfirst_batch_states_collection = ee.FeatureCollection(\"TIGER/2018/States\")\\\n.filter(ee.Filter.inList('STATEFP', first_batch_us_state_code))\n\n# create feature collection with second half of states\nsecond_batch_states_collection = ee.FeatureCollection(\"TIGER/2018/States\")\\\n.filter(ee.Filter.inList('STATEFP', second_batch_us_state_code))\n\n\ntotal_no2_first_batch = first_batch_states_collection.map(\n    lambda feature: calc_each_states_total_no2(feature))\ntotal_no2_second_batch = second_batch_states_collection.map(\n    lambda feature: calc_each_states_total_no2(feature))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# start computation and get desired column as featurecollection\nfirst_batch_results = total_no2_first_batch.reduceColumns(\n    reducer=ee.Reducer.toList().repeat(4),\n    selectors=['NAME', 'STATEFP', 'STUSPS', 'total_no2']).getInfo()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"second_batch_results = total_no2_second_batch.reduceColumns(\n    reducer=ee.Reducer.toList().repeat(4),\n    selectors=['NAME', 'STATEFP', 'STUSPS', 'total_no2']).getInfo()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='process_evaluation'></a>\n### 3. Process for evaluation\nFinally, we have calculated emission for all states. But to evaluate against bottom-up emission we have to process and merge data with bottom-up data. First, we will convert satellite data into pandas dataframe and then we will merge it with bottom-up emission data. **After that we convert emission data from ton to lb and calculate emission factor by dividing emission by total electricity generation. **"},{"metadata":{"trusted":true},"cell_type":"code","source":"# create empty dataframe for storing featurecollection from computaion\ndf_first_batch_total_no2 = pd.DataFrame(columns=['NAME', 'STATEFP', 'STUSPS', 'total_no2'])\ndf_second_batch_total_no2 = pd.DataFrame(columns=['NAME', 'STATEFP', 'STUSPS', 'total_no2'])\n\ndef convert_feature_collection_results_to_df(results, dataframe):\n    dataframe['NAME'] = pd.Series(results.get('list')[0])\n    dataframe['STATEFP'] = pd.Series(results.get('list')[1])\n    dataframe['STUSPS'] = pd.Series(results.get('list')[2])\n    dataframe['total_no2'] = pd.Series(results.get('list')[3])\n    return dataframe","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# convert our calculated featurecollection to pandas dataframe\ndf_first_batch_total_no2 = convert_feature_collection_results_to_df(\n    first_batch_results, df_first_batch_total_no2)\ndf_second_batch_total_no2 = convert_feature_collection_results_to_df(\n    second_batch_results, df_second_batch_total_no2)\n\ndf_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":"# import epa 2018 yearly emission data and merge with our calculated emission data\n\ncolumns = ['FIPSST', 'STNGENAN', 'STNGENOZ', 'STNOXAN', 'STNOXOZ', 'STNOXRTA', 'STNOXRTO']\ndf_2018_epa_emission = pd.read_csv(\n    '../input/state-wise-2018-epa-data/state_wise_epa_data_2018_version_3.csv',\n    header=0, skiprows=1, usecols=columns)\ndf_2018_epa_emission.columns = ['STATEFP', 'annual net electricity generation',\n                                'ozone net electricity generation', 'annual nox emissions',\n                                'ozone nox emissions', 'annual nox output rate',\n                                'ozone nox output rate']\ndf_total_no2['STATEFP'] = df_total_no2['STATEFP'].astype(int)\n\ndf_total_no2_with_bottom_up = df_total_no2.merge(\n    df_2018_epa_emission, on='STATEFP', how='inner').sort_values('STATEFP')\n\n# convert ton to lb, because ground e.f are in lb/MWH\n\ndf_total_no2_with_bottom_up['total_no2_lb'] = df_total_no2_with_bottom_up['total_no2'] * 2000\n\n# calculate ozone season emission factor\ndf_total_no2_with_bottom_up['ozone_emission_factor'] = \\\ndf_total_no2_with_bottom_up['total_no2_lb'] / \\\ndf_total_no2_with_bottom_up['ozone net electricity generation']\n\ndf_total_no2_with_bottom_up.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Evaluation\nAfter a lot of work, we have finally come to the evaluation part. In this part, we are going to evaluate our methodology against bottom-up emission to see how our model performance. How we like to evaluate is described in the upper section. Now we are going to compare our calculated satellite image form 2018-06-28 to 2018-11-28 to bottom-up emission data from EPA of 2018. We have already discussed the issue in the upper section. So let's just see what happens. "},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = go.Figure()\nfig.add_trace(go.Scatter(\n    x=df_total_no2_with_bottom_up['total_no2'].values.tolist(),\n    y=df_total_no2_with_bottom_up['NAME'].values.tolist(),\n    marker=dict(color=\"crimson\", size=12),\n    mode=\"markers\",\n    name=\"Satellite (2018-06-28 to 2018-11-28)\",\n))\n\nfig.add_trace(go.Scatter(\n    x=df_total_no2_with_bottom_up['ozone nox emissions'].values.tolist(),\n    y=df_total_no2_with_bottom_up['NAME'].values.tolist(),\n    marker=dict(color=\"gold\", size=12),\n    mode=\"markers\",\n    name=\"Bottom up (Ozone season 2018)\"\n))\n\nfig.update_layout(title=\"Satellite emissions vs bottom up emissions for ozone season in 2018\",\n                  xaxis_title=\"Emissions (tons)\",autosize=False,\n                  width=700,height=1200,\n                  yaxis_title=\"State names\", yaxis={'nticks':60})\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = go.Figure()\nfig.add_trace(go.Scatter(\n    x=df_total_no2_with_bottom_up['ozone_emission_factor'].values.tolist(),\n    y=df_total_no2_with_bottom_up['NAME'].values.tolist(),\n    marker=dict(color=\"crimson\", size=12), \n    mode=\"markers\",\n    name=\"Satellite (2018-06-28 to 2018-11-28)\",\n))\n\nfig.add_trace(go.Scatter(\n    x=df_total_no2_with_bottom_up['ozone nox output rate'].values.tolist(),\n    y=df_total_no2_with_bottom_up['NAME'].values.tolist(),\n    marker=dict(color=\"gold\", size=12),\n    mode=\"markers\",\n    name=\"Bottom up (Ozone season 2018)\"\n))\n\nfig.update_layout(title=\"Satellite emissions factor vs bottom up emission factor\",\n                  xaxis_title=\"Emissions factor(lb/MWh)\",autosize=False,\n                  width=700,height=1200,\n                  yaxis_title=\"State names\", yaxis={'nticks':60})\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Wow! we are seeing a very astonishing result. The satellite emission data correlate heavily with bottom-up data. We can see for each state our satellite image data produces very similar data. Of course, there is some ambiguity of the date of calculation, for that reason we are seeing some error. But overall we are confident that our methodology has successfully calculated emission from satellite images. "},{"metadata":{},"cell_type":"markdown","source":"# Conclusion\nWe have finally implemented our methodology for calculating the emission factor for sub-national regions and individual states. Also, we successfully evaluate our results and found very promising results. In the next notebook, we are going to see how to use this methodology for individual power plants."}],"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}